--- /dev/null
--- /dev/null
++cmake_minimum_required(VERSION 2.8.6)
++
++project( MathGL )
++
++if(NOT CMAKE_BUILD_TYPE)
++ set(CMAKE_BUILD_TYPE "Release" CACHE STRING "Choose the type of build, options are:
++ None(CMAKE_CXX_FLAGS or CMAKE_C_FLAGS used) Debug Release RelWithDebInfo MinSizeRel." FORCE)
++endif(NOT CMAKE_BUILD_TYPE)
++
++set(CMAKE_VERBOSE_MAKEFILE ON)
++set(MathGL_VERSION_MAJOR 2)
++set(MathGL_VERSION_MINOR 2.0)
++set(MathGL_SOVERSION 7.1.0)
++
++set(CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} -Wl,-z,relro")
++set(CMAKE_SHARED_LINKER_FLAGS "${CMAKE_SHARED_LINKER_FLAGS} -Wl,-z,relro")
++set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -Wl,-z,relro")
++
++MACRO(MGL_DEPENDENT_OPTION option doc default depends1 force1 depends2 force2)
++ IF(${option}_ISSET MATCHES "^${option}_ISSET$")
++ SET(${option}_AVAILABLE 1)
++ IF(${force1})
++ FOREACH(d ${depends1})
++ STRING(REGEX REPLACE " +" ";" CMAKE_DEPENDENT_OPTION_DEP "${d}")
++ IF(${CMAKE_DEPENDENT_OPTION_DEP})
++ ELSE(${CMAKE_DEPENDENT_OPTION_DEP})
++ SET(${option}_AVAILABLE 0)
++ SET(depends1_AVAILABLE 1)
++ ENDIF(${CMAKE_DEPENDENT_OPTION_DEP})
++ ENDFOREACH(d)
++ ENDIF(${force1})
++ IF(${force2})
++ FOREACH(d ${depends2})
++ STRING(REGEX REPLACE " +" ";" CMAKE_DEPENDENT_OPTION_DEP "${d}")
++ IF(${CMAKE_DEPENDENT_OPTION_DEP})
++ ELSE(${CMAKE_DEPENDENT_OPTION_DEP})
++ SET(${option}_AVAILABLE 0)
++ SET(depends2_AVAILABLE 1)
++ ENDIF(${CMAKE_DEPENDENT_OPTION_DEP})
++ ENDFOREACH(d)
++ ENDIF(${force2})
++ IF(${option}_AVAILABLE)
++ OPTION(${option} "${doc}" "${default}")
++ SET(${option} "${${option}}" CACHE BOOL "${doc}" FORCE)
++ ELSE(${option}_AVAILABLE)
++ IF(${option} MATCHES "^${option}$")
++ ELSE(${option} MATCHES "^${option}$")
++ SET(${option} "${${option}}" CACHE INTERNAL "${doc}")
++ ENDIF(${option} MATCHES "^${option}$")
++ IF(depends1_AVAILABLE)
++ SET(${option} OFF)
++ ELSEIF(depends2_AVAILABLE)
++ SET(${option} ON)
++ ENDIF(depends1_AVAILABLE)
++ ENDIF(${option}_AVAILABLE)
++ ELSE(${option}_ISSET MATCHES "^${option}_ISSET$")
++ SET(${option} "${${option}_ISSET}")
++ ENDIF(${option}_ISSET MATCHES "^${option}_ISSET$")
++ENDMACRO(MGL_DEPENDENT_OPTION)
++include(CMakeDependentOption)
++
++set(MGL_LIB_INSTALL_DIR "lib" CACHE STRING "Set library install directory")
++
++option(enable-double "Enable double precision in MathGL library" ON)
++option(enable-simple "Slightly increase drawing speed but disable mglDataA class")
++<<<<<<< HEAD
++option(enable-mpi "Enable mpi" ON)
++option(enable-doc "Enable documentation building" ON)
++option(enable-all "Enable all core features" ON)
++option(enable-all-widgets "Enable all Widgets" ON)
++option(enable-all-swig "Enable all SWIG based interfaces" ON)
++option(enable-pthread "Enable POSIX threads support" ON)
++option(enable-lgpl "Enable only LGPL part of MathGL" OFF)
++option(enable-mgl2 "Use names 'libmgl2-*' instead of 'libmgl-*'" OFF)
++=======
++option(enable-mpi "Enable mpi")
++option(enable-opengl "Enable OpenGL support" ON)
++option(enable-all-docs "Enable all documentation building")
++#option(enable-doc "Enable documentation building")
++option(enable-all "Enable all core features")
++option(enable-all-widgets "Enable all Widgets")
++option(enable-all-swig "Enable all SWIG based interfaces")
++option(enable-pthread "Enable POSIX threads support" OFF)
++option(enable-openmp "Enable POSIX threads support" ON)
++option(enable-lgpl "Enable only LGPL part of MathGL")
++option(enable-mgl2 "Use names 'libmgl2-*' instead of 'libmgl-*'")
++>>>>>>> upstream/2.2
++#option(enable-ltdl "Enable loading modules support")
++CMAKE_DEPENDENT_OPTION(enable-doc-site "Enable HTML documentation for website" OFF "NOT enable-all-docs" ON)
++CMAKE_DEPENDENT_OPTION(enable-doc-html "Enable HTML documentation" OFF "NOT enable-all-docs" ON)
++CMAKE_DEPENDENT_OPTION(enable-doc-info "Enable INFO documentation" OFF "NOT enable-all-docs" ON)
++CMAKE_DEPENDENT_OPTION(enable-doc-pdf-ru "Enable Russian PDF documentation" OFF "NOT enable-all-docs" ON)
++CMAKE_DEPENDENT_OPTION(enable-doc-pdf-en "Enable English PDF documentation" OFF "NOT enable-all-docs" ON)
++CMAKE_DEPENDENT_OPTION(enable-doc-prc "Enable PDF samples for HTML docs" OFF "NOT enable-all-docs" ON)
++CMAKE_DEPENDENT_OPTION(enable-doc-json "Enable JSON samples for HTML docs" OFF "NOT enable-all-docs" ON)
++option(enable-texi2html "Use texi2html (obsolete package) instead of texi2any" OFF)
++
++CMAKE_DEPENDENT_OPTION(enable-zlib "Enable zlib support" ON "NOT enable-all" ON)
++CMAKE_DEPENDENT_OPTION(enable-png "Enable png support" ON "NOT enable-all" ON)
++<<<<<<< HEAD
++CMAKE_DEPENDENT_OPTION(enable-jpeg "Enable jpeg support" ON "NOT enable-all" ON)
++MGL_DEPENDENT_OPTION(enable-gsl "Enable gsl support" ON "NOT enable-lgpl" ON "NOT enable-all" ON)
++MGL_DEPENDENT_OPTION(enable-hdf4 "Enable hdf4 support" ON "NOT enable-lgpl" ON "NOT enable-all" ON)
++MGL_DEPENDENT_OPTION(enable-hdf5 "Enable hdf5 support" ON "NOT enable-lgpl" ON "NOT enable-all" ON)
++CMAKE_DEPENDENT_OPTION(enable-pdf "Enable pdf support" ON "NOT enable-all" ON)
++CMAKE_DEPENDENT_OPTION(enable-gif "Enable gif support" ON "NOT enable-all" ON)
++CMAKE_DEPENDENT_OPTION(enable-opengl "Enable OpenGL support" ON "NOT enable-mpi" ON)
++MGL_DEPENDENT_OPTION(enable-glut "Enable glut support" ON "NOT enable-lgpl;NOT enable-mpi" OFF "NOT enable-all-widgets" OFF)
++MGL_DEPENDENT_OPTION(enable-fltk "Enable fltk widget" ON "NOT enable-lgpl;NOT enable-mpi" OFF "NOT enable-all-widgets" OFF)
++CMAKE_DEPENDENT_OPTION(enable-wx "Enable wxWidget widget" ON "NOT enable-lgpl;NOT enable-mpi" ON)
++MGL_DEPENDENT_OPTION(enable-qt "Enable Qt4 widget" ON "NOT enable-lgpl;NOT enable-mpi" OFF "NOT enable-all-widgets" OFF)
++MGL_DEPENDENT_OPTION(enable-python "Enable python interface" ON "NOT enable-lgpl" ON "NOT enable-all-swig" OFF)
++MGL_DEPENDENT_OPTION(enable-octave "Enable octave interface" ON "NOT enable-lgpl" ON "NOT enable-all-swig" OFF)
++MGL_DEPENDENT_OPTION(enable-octave-install "Octave interface will install for all users" ON "NOT enable-lgpl" ON "NOT enable-all-swig" OFF)
++=======
++CMAKE_DEPENDENT_OPTION(enable-jpeg "Enable jpeg support" OFF "NOT enable-all" ON)
++MGL_DEPENDENT_OPTION(enable-gsl "Enable gsl support" OFF "NOT enable-lgpl" ON "NOT enable-all" ON)
++MGL_DEPENDENT_OPTION(enable-hdf4 "Enable hdf4 support" OFF "NOT enable-lgpl" ON "NOT enable-all" ON)
++MGL_DEPENDENT_OPTION(enable-hdf5 "Enable hdf5 support" OFF "NOT enable-lgpl" ON "NOT enable-all" ON)
++CMAKE_DEPENDENT_OPTION(enable-pdf "Enable pdf support" OFF "NOT enable-all" ON)
++CMAKE_DEPENDENT_OPTION(enable-gif "Enable gif support" OFF "NOT enable-all" ON)
++MGL_DEPENDENT_OPTION(enable-glut "Enable glut support" OFF "NOT enable-lgpl" ON "NOT enable-all-widgets" ON)
++MGL_DEPENDENT_OPTION(enable-fltk "Enable fltk widget" OFF "NOT enable-lgpl" ON "NOT enable-all-widgets" ON)
++CMAKE_DEPENDENT_OPTION(enable-wx "Enable wxWidget widget" OFF "NOT enable-lgpl" OFF)
++MGL_DEPENDENT_OPTION(enable-qt "Enable Qt4 widget" OFF "NOT enable-lgpl" ON "NOT enable-all-widgets" ON)
++MGL_DEPENDENT_OPTION(enable-python "Enable python interface" OFF "NOT enable-lgpl" ON "NOT enable-all-swig" ON)
++MGL_DEPENDENT_OPTION(enable-lua "Enable Lua (v.5.1) interface" OFF "NOT enable-lgpl" ON "NOT enable-all-swig" ON)
++MGL_DEPENDENT_OPTION(enable-octave "Enable octave interface" OFF "NOT enable-lgpl" ON "NOT enable-all-swig" ON)
++MGL_DEPENDENT_OPTION(enable-octave-install "Octave interface will install for all users" ON "NOT enable-lgpl" ON "NOT enable-all-swig" ON)
++>>>>>>> upstream/2.2
++
++include_directories( ${MathGL_SOURCE_DIR}/include ${MathGL_BINARY_DIR}/include)
++set(MGL_INCLUDE_PATH "${CMAKE_INSTALL_PREFIX}/include/mgl2")
++set(MGL_CGI_PATH "${CMAKE_INSTALL_PREFIX}/share/mathgl" CACHE STRING "Set CGI install directory")
++
++if(NOT WIN32)
++# set(MGL_CGI_PATH "${CMAKE_INSTALL_PREFIX}/share/mathgl")
++ set(MGL_DATA_PATH "${CMAKE_INSTALL_PREFIX}/share/mathgl")
++ set(MGL_DOC_PATH "${CMAKE_INSTALL_PREFIX}/share/doc/mathgl")
++ set(MGL_MAN_PATH "${CMAKE_INSTALL_PREFIX}/share/man")
++ set(MGL_INFO_PATH "${CMAKE_INSTALL_PREFIX}/share/info")
++ set(MGL_FONT_PATH "${MGL_DATA_PATH}/fonts")
++ install(DIRECTORY fonts/ DESTINATION ${MGL_FONT_PATH}
++ PATTERN ".svn" EXCLUDE
++ PATTERN "*.vfm")
++endif(NOT WIN32)
++
++include(CheckFunctionExists)
++include(CMakePushCheckState)
++include(TestBigEndian)
++
++TEST_BIG_ENDIAN(WORDS_BIGENDIAN)
++if(WORDS_BIGENDIAN)
++ ADD_DEFINITIONS("-DWORDS_BIGENDIAN")
++endif(WORDS_BIGENDIAN)
++
++CHECK_FUNCTION_EXISTS(sin MGL_SIN)
++CHECK_FUNCTION_EXISTS(memrchr HAVE_MEMRCHR)
++if(NOT MGL_SIN)
++ cmake_push_check_state()
++ set(CMAKE_REQUIRED_LIBRARIES ${CMAKE_REQUIRED_LIBRARIES} m)
++ CHECK_FUNCTION_EXISTS(sin MGL_SIN_M)
++ cmake_pop_check_state()
++ if(MGL_SIN_M)
++ set(M_LIB m)
++ elseif(MGL_SIN_M)
++ message(SEND_ERROR "Math library not found")
++ endif(MGL_SIN_M)
++endif(NOT MGL_SIN)
++if(HAVE_MEMRCHR)
++ ADD_DEFINITIONS("-DHAVE_MEMRCHR")
++endif(HAVE_MEMRCHR)
++
++if(enable-double)
++ set(MGL_USE_DOUBLE 1)
++else(enable-double)
++ set(MGL_USE_DOUBLE 0)
++endif(enable-double)
++
++if(enable-simple)
++ set(MGL_NO_DATA_A 1)
++message(STATUS "Class mglDataA is switched off.")
++else(enable-simple)
++ set(MGL_NO_DATA_A 0)
++endif(enable-simple)
++
++if(enable-openmp)
++ find_package(OpenMP)
++ if(OPENMP_FOUND)
++ set(MGL_HAVE_OMP 1)
++ set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
++ set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
++ else(OPENMP_FOUND)
++ message(SEND_ERROR "Couldn't find OpenMP. You can enable POSIX threads instead.")
++ set(MGL_HAVE_OMP 0)
++ endif(OPENMP_FOUND)
++else(enable-openmp)
++ set(MGL_HAVE_OMP 0)
++endif(enable-openmp)
++
++if(enable-mpi)
++ set(MGL_HAVE_MPI 1)
++ find_package(MPI REQUIRED)
++ set(CMAKE_CXX_COMPILE_FLAGS ${CMAKE_CXX_COMPILE_FLAGS} ${MPI_COMPILE_FLAGS})
++ set(CMAKE_CXX_LINK_FLAGS ${CMAKE_CXX_LINK_FLAGS} ${MPI_LINK_FLAGS})
++ include_directories(${MPI_INCLUDE_PATH})
++else(enable-mpi)
++ set(MGL_HAVE_MPI 0)
++endif(enable-mpi)
++
++if(enable-pthread)
++ set(MGL_HAVE_PTHREAD 1)
++ include(FindThreads)
++ if(NOT CMAKE_USE_PTHREADS_INIT)
++ message(SEND_ERROR "Couldn't find POSIX threads library.")
++ endif(NOT CMAKE_USE_PTHREADS_INIT)
++else(enable-pthread)
++ set(MGL_HAVE_PTHREAD 0)
++endif(enable-pthread)
++
++if(enable-gsl)
++ set(MGL_HAVE_GSL 1)
++ find_library(GSL_LIB gsl)
++ find_library(GSL_CBLAS_LIB gslcblas)
++ find_path(GSL_INCLUDE_DIR gsl/gsl_fft_complex.h)
++ if(NOT GSL_LIB OR NOT GSL_CBLAS_LIB OR NOT GSL_INCLUDE_DIR)
++ message(SEND_ERROR "${GSL_LIB}")
++ message(SEND_ERROR "${GSL_CBLAS_LIB}")
++ message(SEND_ERROR "${GSL_INCLUDE_DIR}")
++ message(SEND_ERROR "Couldn't find GSL libraries.")
++ endif(NOT GSL_LIB OR NOT GSL_CBLAS_LIB OR NOT GSL_INCLUDE_DIR)
++else(enable-gsl)
++ set(MGL_HAVE_GSL 0)
++endif(enable-gsl)
++
++#if(enable-all OR enable-ltdl)
++# set(MGL_HAVE_LTDL 1)
++# find_library(LTDL_LIB ltdl)
++# find_path(LTDL_INCLUDE_DIR ltdl.h)
++# if(NOT LTDL_LIB OR NOT LTDL_INCLUDE_DIR)
++# message(SEND_ERROR "${LTDL_LIB}")
++# message(SEND_ERROR "${LTDL_INCLUDE_DIR}")
++# message(SEND_ERROR "Couldn't find LTDL library.")
++# endif(NOT LTDL_LIB OR NOT LTDL_INCLUDE_DIR)
++#else(enable-all OR enable-ltdl)
++# set(MGL_HAVE_LTDL 0)
++#endif(enable-all OR enable-ltdl)
++
++if(enable-hdf4)
++ set(MGL_HAVE_HDF4 1)
++ find_library(HDF4_LIB df)
++ find_library(HDF4MF_LIB mfhdf)
++ find_path(HDF4_INCLUDE_DIR hdf/mfhdf.h)
++ if(NOT HDF4_LIB OR NOT HDF4MF_LIB OR NOT HDF4_INCLUDE_DIR)
++ message(SEND_ERROR "${HDF4_LIB}")
++ message(SEND_ERROR "${HDF4MF_LIB}")
++ message(SEND_ERROR "${HDF4_INCLUDE_DIR}")
++ message(SEND_ERROR "Couldn't find HDF4 libraries.")
++ endif(NOT HDF4_LIB OR NOT HDF4MF_LIB OR NOT HDF4_INCLUDE_DIR)
++else(enable-hdf4)
++ set(MGL_HAVE_HDF4 0)
++endif(enable-hdf4)
++
++if(enable-hdf5)
++ set(MGL_HAVE_HDF5 1)
++ include(FindHDF5)
++ if(NOT HDF5_FOUND)
++ message(SEND_ERROR "Couldn't find HDF5 library.")
++ endif(NOT HDF5_FOUND)
++else(enable-hdf5)
++ set(MGL_HAVE_HDF5 0)
++endif(enable-hdf5)
++
++if(enable-jpeg)
++ set(MGL_HAVE_JPEG 1)
++ include(FindJPEG)
++ if(NOT JPEG_FOUND)
++ message(SEND_ERROR "Couldn't find JPEG library.")
++ endif(NOT JPEG_FOUND)
++else(enable-jpeg)
++ set(MGL_HAVE_JPEG 0)
++endif(enable-jpeg)
++
++
++if(enable-zlib)
++ set(MGL_HAVE_ZLIB 1)
++ include(FindZLIB)
++ if(NOT ZLIB_FOUND)
++ message(SEND_ERROR "Couldn't find ZLib library.")
++ endif(NOT ZLIB_FOUND)
++else(enable-zlib)
++ set(MGL_HAVE_ZLIB 0)
++endif(enable-zlib)
++
++if(enable-png)
++ set(MGL_HAVE_PNG 1)
++ if(NOT MGL_HAVE_ZLIB)
++ message(SEND_ERROR "You have to enable ZLib if you plan to use PNG export.")
++ endif(NOT MGL_HAVE_ZLIB)
++ include(FindPNG)
++ if(NOT PNG_FOUND)
++ message(SEND_ERROR "Couldn't find PNG library.")
++ endif(NOT PNG_FOUND)
++else(enable-png)
++ set(MGL_HAVE_PNG 0)
++endif(enable-png)
++
++
++if(enable-pdf)
++ set(MGL_HAVE_PDF 1)
++ if(NOT MGL_HAVE_PNG)
++ message(SEND_ERROR "You have to enable PNG if you plan to use PDF export.")
++ endif(NOT MGL_HAVE_PNG)
++ find_library(HPDF_LIB hpdf)
++ if(NOT HPDF_LIB)
++ message(SEND_ERROR "Couldn't find libHaru or libhpdf.")
++ endif(NOT HPDF_LIB)
++ find_path(HPDF_INCLUDE_DIR hpdf_u3d.h)
++ if(NOT HPDF_INCLUDE_DIR)
++ message(SEND_ERROR "Couldn't find headers of 3d-enabled version of libhpdf.")
++ endif(NOT HPDF_INCLUDE_DIR)
++else(enable-pdf)
++ set(MGL_HAVE_PDF 0)
++endif(enable-pdf)
++
++if(enable-gif)
++ set(MGL_HAVE_GIF 1)
++ include(FindGIF)
++ if(NOT GIF_FOUND)
++ message(SEND_ERROR "Couldn't find GIF library.")
++ endif(NOT GIF_FOUND)
++else(enable-gif)
++ set(MGL_HAVE_GIF 0)
++endif(enable-gif)
++
++if(enable-opengl)
++ set(MGL_HAVE_OPENGL 1)
++ include(FindOpenGL)
++ if(NOT OPENGL_FOUND)
++ message(SEND_ERROR "Couldn't find OpenGL libraries.")
++ endif(NOT OPENGL_FOUND)
++else(enable-opengl)
++ set(MGL_HAVE_OPENGL 0)
++endif(enable-opengl)
++
++if(enable-glut)
++ set(MGL_HAVE_GLUT 1)
++ if(NOT MGL_HAVE_OPENGL)
++ message(SEND_ERROR "You have to enable OpenGL if you plan to use GLUT.")
++ endif(NOT MGL_HAVE_OPENGL)
++ include(FindGLUT)
++ if(NOT GLUT_FOUND)
++ message(SEND_ERROR "Couldn't find GLUT library.")
++ endif(NOT GLUT_FOUND)
++else(enable-glut)
++ set(MGL_HAVE_GLUT 0)
++endif(enable-glut)
++
++if(enable-fltk)
++ set(MGL_HAVE_FLTK 1)
++ include(FindFLTK)
++ if(NOT FLTK_FOUND)
++ message(SEND_ERROR "Couldn't find FLTK library.")
++ endif(NOT FLTK_FOUND)
++else(enable-fltk)
++ set(MGL_HAVE_FLTK 0)
++endif(enable-fltk)
++
++if(enable-wx)
++ set(MGL_HAVE_WX 1)
++ FIND_PACKAGE(wxWidgets COMPONENTS base core gl)
++ if(NOT wxWidgets_FOUND)
++ message(SEND_ERROR "Couldn't find wxWidgets library.")
++ endif(NOT wxWidgets_FOUND)
++else(enable-wx)
++ set(MGL_HAVE_WX 0)
++endif(enable-wx)
++
++if(enable-qt)
++ set(MGL_HAVE_QT 1)
++ FIND_PACKAGE(Qt4 4.8 REQUIRED QtCore QtGui QtNetwork QtWebKit)
++ if(NOT QT4_FOUND)
++ message(SEND_ERROR "Couldn't find Qt4 library.")
++ endif(NOT QT4_FOUND)
++else(enable-qt)
++ set(MGL_HAVE_QT 0)
++endif(enable-qt)
++
++if(enable-python)
++ set(MGL_HAVE_SWIG 1)
++ set(MGL_HAVE_PYTHON 1)
++ FIND_PACKAGE(PythonInterp)
++ if(NOT PYTHONINTERP_FOUND)
++ message(SEND_ERROR "Couldn't find python interpreter.")
++ endif(NOT PYTHONINTERP_FOUND)
++ FIND_PACKAGE(PythonLibs)
++ if(NOT PYTHONLIBS_FOUND)
++ message(SEND_ERROR "Couldn't find python development libraries.")
++ endif(NOT PYTHONLIBS_FOUND)
++ execute_process(
++ COMMAND ${PYTHON_EXECUTABLE} -c "import numpy; print numpy.get_include()"
++ OUTPUT_VARIABLE NUMPY_INCLUDE_PATH
++ RESULT_VARIABLE NUMPY_ERR
++ OUTPUT_STRIP_TRAILING_WHITESPACE
++ )
++ if(NOT NUMPY_INCLUDE_PATH)
++ message(SEND_ERROR "Couldn't find numpy.")
++ endif(NOT NUMPY_INCLUDE_PATH)
++else(enable-python)
++ set(MGL_HAVE_PYTHON 0)
++endif(enable-python)
++
++if(enable-lua)
++ set(MGL_HAVE_SWIG 1)
++ set(MGL_HAVE_LUA 1)
++ INCLUDE(FindLua51)
++ if(NOT LUA51_FOUND)
++ message(SEND_ERROR "Couldn't find Lua 5.1 library.")
++ endif(NOT LUA51_FOUND)
++else(enable-lua)
++ set(MGL_HAVE_LUA 0)
++endif(enable-lua)
++
++
++if(enable-octave)
++ set(MGL_HAVE_SWIG 1)
++ set(MGL_HAVE_OCTAVE 1)
++ find_program(oct_prog octave-config)
++ if(NOT oct_prog)
++ message(SEND_ERROR "Couldn't find octave-config needed for octave interfaces compiling.")
++ endif(NOT oct_prog)
++ find_program(oct_exec octave)
++ if(NOT oct_exec)
++ message(SEND_ERROR "Couldn't find octave needed for octave interfaces compiling.")
++ endif(NOT oct_exec)
++ find_program(oct_mk mkoctfile)
++ if(NOT oct_mk)
++ message(SEND_ERROR "Couldn't find mkoctfile needed for octave interfaces compiling.")
++ endif(NOT oct_mk)
++ find_program(oct_tar tar)
++ if(NOT oct_tar)
++ message(SEND_ERROR "Couldn't find tar needed for octave interfaces creation.")
++ endif(NOT oct_tar)
++else(enable-octave)
++ set(MGL_HAVE_OCTAVE 0)
++endif(enable-octave)
++
++if(enable-doc-info)
++ set(MGL_HAVE_DOC_INFO 1)
++ find_program(findmi makeinfo)
++ if(NOT findmi)
++ message(SEND_ERROR "Couldn't find makeinfo needed for documentation building.")
++ endif(NOT findmi)
++else(enable-doc-info)
++ set(MGL_HAVE_DOC_INFO 0)
++endif(enable-doc-info)
++
++if(enable-doc-html OR enable-doc-site)
++ if(enable-texi2html)
++ find_program(findth texi2html)
++ if(NOT findth)
++ message(SEND_ERROR "Couldn't find texi2html needed for documentation building.")
++ endif(NOT findth)
++ else(enable-texi2html)
++ find_program(findth texi2any)
++ if(NOT findth)
++ message(SEND_ERROR "Couldn't find texi2any needed for documentation building.")
++ endif(NOT findth)
++ endif(enable-texi2html)
++endif(enable-doc-html OR enable-doc-site)
++
++if(enable-texi2html)
++ set(site_en ${CMAKE_BINARY_DIR}/texinfo/doc_en/doc_en.html)
++ set(site_ru ${CMAKE_BINARY_DIR}/texinfo/doc_ru/doc_ru.html)
++ set(th_opt )
++else(enable-texi2html)
++ set(th_opt --html)
++ set(site_en ${CMAKE_BINARY_DIR}/texinfo/doc_en/index.html)
++ set(site_ru ${CMAKE_BINARY_DIR}/texinfo/doc_ru/index.html)
++endif(enable-texi2html)
++
++if(enable-doc-html)
++ set(MGL_HAVE_DOC_HTML 1)
++else(enable-doc-html)
++ set(MGL_HAVE_DOC_HTML 0)
++endif(enable-doc-html)
++
++
++if(enable-doc-site)
++ set(MGL_HAVE_DOC_SITE 1)
++else(enable-doc-site)
++ set(MGL_HAVE_DOC_SITE 0)
++endif(enable-doc-site)
++
++if(enable-doc-pdf-ru)
++ set(MGL_HAVE_DOC_PDF_RU 1)
++ find_program(findtp texi2pdf)
++ if(NOT findtp)
++ message(SEND_ERROR "Couldn't find texi2pdf needed for documentation building.")
++ endif(NOT findtp)
++else(enable-doc-pdf-ru)
++ set(MGL_HAVE_DOC_PDF_RU 0)
++endif(enable-doc-pdf-ru)
++
++if(enable-doc-pdf-en)
++ set(MGL_HAVE_DOC_PDF_EN 1)
++ find_program(findtp texi2pdf)
++ if(NOT findtp)
++ message(SEND_ERROR "Couldn't find texi2pdf needed for documentation building.")
++ endif(NOT findtp)
++else(enable-doc-pdf-en)
++ set(MGL_HAVE_DOC_PDF_EN 0)
++endif(enable-doc-pdf-en)
++
++if(enable-doc-json)
++ set(MGL_HAVE_DOC_JSON 1)
++else(enable-doc-json)
++ set(MGL_HAVE_DOC_JSON 0)
++endif(enable-doc-json)
++
++if(enable-doc-prc)
++ set(MGL_HAVE_DOC_PRC 1)
++ if(NOT enable-pdf)
++ message(SEND_ERROR "You need to enable pdf support for MathGL.")
++ endif(NOT enable-pdf)
++else(enable-doc-prc)
++ set(MGL_HAVE_DOC_PRC 0)
++endif(enable-doc-prc)
++
++if(MGL_HAVE_SWIG)
++ FIND_PACKAGE(SWIG)
++ if(NOT SWIG_FOUND)
++ message(SEND_ERROR "Couldn't find swig needed for interfaces compiling.")
++ endif(NOT SWIG_FOUND)
++ INCLUDE(${SWIG_USE_FILE})
++endif(MGL_HAVE_SWIG)
++
++if(UNIX)
++ add_definitions(-DNO_COLOR_ARRAY)
++ add_definitions(-D_FORTIFY_SOURCE=2)
++endif(UNIX)
++
++if(WIN32)
++ add_definitions(-DWIN32)
++endif(WIN32)
++
++add_subdirectory( src )
++add_subdirectory( widgets )
++add_subdirectory( include )
++add_subdirectory( udav )
++add_subdirectory( json )
++#add_subdirectory( mgllab )
++add_subdirectory( lang )
++if(NOT MSVC AND NOT BORLAND)
++add_subdirectory( utils )
++add_subdirectory( examples )
++
++if(MGL_HAVE_DOC_HTML OR MGL_HAVE_DOC_SITE OR MGL_HAVE_DOC_INFO OR MGL_HAVE_DOC_PDF_RU OR MGL_HAVE_DOC_PDF_EN )
++add_subdirectory( texinfo )
++endif(MGL_HAVE_DOC_HTML OR MGL_HAVE_DOC_SITE OR MGL_HAVE_DOC_INFO OR MGL_HAVE_DOC_PDF_RU OR MGL_HAVE_DOC_PDF_EN )
++
++endif(NOT MSVC AND NOT BORLAND)
--- /dev/null
--- /dev/null
++/***************************************************************************
++ * fft.cpp is part of Math Graphic Library
++ * Copyright (C) 2007-2012 Alexey Balakin <mathgl.abalakin@gmail.ru> *
++ * *
++ * This program is free software; you can redistribute it and/or modify *
++ * it under the terms of the GNU Library General Public License as *
++ * published by the Free Software Foundation; either version 3 of the *
++ * License, or (at your option) any later version. *
++ * *
++ * This program is distributed in the hope that it will be useful, *
++ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
++ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
++ * GNU General Public License for more details. *
++ * *
++ * You should have received a copy of the GNU Library General Public *
++ * License along with this program; if not, write to the *
++ * Free Software Foundation, Inc., *
++ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
++ ***************************************************************************/
++#include "mgl2/datac.h"
++#include "mgl2/data.h"
++#include "mgl2/thread.h"
++#if MGL_HAVE_GSL
++#include <gsl/gsl_fft_complex.h>
++#include <gsl/gsl_dht.h>
++#include <gsl/gsl_sf.h>
++#endif
++//-----------------------------------------------------------------------------
++void MGL_EXPORT mglStartThreadT(void *(*func)(void *), long n, void *a, double *b, const void *v, void **w, const long *p, const void *re, const void *im)
++{
++ if(!func) return;
++#if MGL_HAVE_PTHREAD
++ if(mglNumThr<1) mgl_set_num_thr(0);
++ if(mglNumThr>1)
++ {
++ pthread_t *tmp=new pthread_t[mglNumThr];
++ mglThreadT *par=new mglThreadT[mglNumThr];
++ register long i;
++ for(i=0;i<mglNumThr;i++) // put parameters into the structure
++ { par[i].n=n; par[i].a=a; par[i].v=v; par[i].w=w; par[i].b=b;
++ par[i].p=p; par[i].re=re; par[i].im=im; par[i].id=i; }
++ for(i=0;i<mglNumThr;i++) pthread_create(tmp+i, 0, func, par+i);
++ for(i=0;i<mglNumThr;i++) pthread_join(tmp[i], 0);
++ delete []tmp; delete []par;
++ }
++ else
++#endif
++ {
++ mglNumThr = 1;
++ mglThreadT par;
++ par.n=n; par.a=a; par.b=b; par.v=v; par.w=w;
++ par.p=p; par.re=re; par.im=im; par.id=0;
++ func(&par);
++ }
++}
++//-----------------------------------------------------------------------------
++struct mglFFTdata
++{
++ long wnx,wny,wnz; // sizes for FFT
++ long hnx,hny,hnz; // sizes for Hankel
++ void *wtx,*wty,*wtz; // tables for FFT
++ void *htx,*hty,*htz; // tables for Hankel
++ mglFFTdata() { memset(this,0,sizeof(mglFFTdata)); }
++ ~mglFFTdata() { Clear(); }
++ void Clear()
++ {
++ if(wnx) { wnx=0; mgl_fft_free(wtx,0,0); }
++ if(wny) { wny=0; mgl_fft_free(wty,0,0); }
++ if(wnz) { wnz=0; mgl_fft_free(wtz,0,0); }
++#if MGL_HAVE_GSL
++ if(hnx) { hnx=0; gsl_dht_free((gsl_dht*)htx); }
++ if(hny) { hny=0; gsl_dht_free((gsl_dht*)hty); }
++ if(hnz) { hnz=0; gsl_dht_free((gsl_dht*)htz); }
++#endif
++ }
++} mgl_fft_data;
++void MGL_EXPORT mgl_clear_fft() { mgl_fft_data.Clear(); }
++//-----------------------------------------------------------------------------
++MGL_EXPORT void *mgl_fft_alloc_thr(long n)
++{
++#if MGL_HAVE_GSL
++ return gsl_fft_complex_workspace_alloc(n);
++#else
++ return new double[2*n];
++#endif
++}
++//-----------------------------------------------------------------------------
++MGL_EXPORT void *mgl_fft_alloc(long n, void **space, long nthr)
++{
++ if(space && nthr>0) for(long i=0;i<nthr;i++) space[i] = mgl_fft_alloc_thr(n);
++#if MGL_HAVE_GSL
++ return gsl_fft_complex_wavetable_alloc(n);
++#else
++ double *c = new double[2*n*n];
++#pragma omp parallel for collapse(2)
++ for(long i=0;i<n;i++) for(long j=0;j<n;j++)
++ { c[2*(i+n*j)]=cos(2*M_PI*i*j/n); c[2*(i+n*j)+1]=-sin(2*M_PI*i*j/n); }
++ return c;
++#endif
++}
++//-----------------------------------------------------------------------------
++void MGL_EXPORT mgl_fft_free_thr(void *ws)
++{
++#if MGL_HAVE_GSL
++ gsl_fft_complex_workspace_free((gsl_fft_complex_workspace*)ws);
++#else
++ delete []((double*)ws);
++#endif
++}
++//-----------------------------------------------------------------------------
++void MGL_EXPORT mgl_fft_free(void *wt, void **ws, long nthr)
++{
++ if(ws && nthr>0) for(long i=0;i<nthr;i++) mgl_fft_free_thr(ws[i]);
++#if MGL_HAVE_GSL
++ gsl_fft_complex_wavetable_free((gsl_fft_complex_wavetable*)wt);
++#else
++ delete []((double*)wt);
++#endif
++}
++//-----------------------------------------------------------------------------
++void MGL_EXPORT mgl_fft(double *x, long s, long n, const void *wt, void *ws, int inv)
++{
++#if MGL_HAVE_GSL
++<<<<<<< HEAD
++gsl_fft_complex_transform(x, s, n, (const gsl_fft_complex_wavetable*)wt, (gsl_fft_complex_workspace*)ws, inv?gsl_fft_backward:gsl_fft_forward);
++=======
++ if(inv) gsl_fft_complex_inverse(x, s, n, (const gsl_fft_complex_wavetable*)wt, (gsl_fft_complex_workspace*)ws);
++ else gsl_fft_complex_forward(x, s, n, (const gsl_fft_complex_wavetable*)wt, (gsl_fft_complex_workspace*)ws);
++>>>>>>> upstream/2.2
++#else // NOTE this is VERY slow!
++ const double *c = (const double *)wt;
++ double *d = (double *)ws, f = inv?1./n:1;
++ memset(d,0,2*n*sizeof(double));
++ if(inv)
++//#pragma omp parallel for // NOTE only 1st for can be used!
++ for(long i=0;i<n;i++) for(long j=0;j<n;j++)
++ {
++ register long ii = 2*(i+n*j), jj = 2*j*s;
++ d[2*i] += x[jj]*c[ii]+x[jj+1]*c[ii+1];
++ d[2*i+1]+= x[jj+1]*c[ii]-x[jj]*c[ii+1];
++ }
++ else
++//#pragma omp parallel for // NOTE only 1st for can be used!
++ for(long i=0;i<n;i++) for(long j=0;j<n;j++)
++ {
++ register long ii = 2*(i+n*j), jj = 2*j*s;
++ d[2*i] += x[jj]*c[ii]-x[jj+1]*c[ii+1];
++ d[2*i+1]+= x[jj+1]*c[ii]+x[jj]*c[ii+1];
++ }
++//#pragma omp parallel for
++ for(long j=0;j<n;j++)
++ { register long jj = 2*j*s; x[jj] = d[2*j]*f; x[jj+1] = d[2*j+1]*f; }
++#endif
++}
++//-----------------------------------------------------------------------------
++MGL_NO_EXPORT void* mgl_fftx(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long nx=t->p[0];
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ void *w = mgl_fft_alloc_thr(nx);
++#pragma omp for nowait
++ for(long i=t->id;i<t->n;i+=mglNumThr)
++ mgl_fft(t->b+2*nx*i, 1, nx, t->v, w, t->p[3]);
++ mgl_fft_free_thr(w);
++ }
++ return 0;
++}
++MGL_NO_EXPORT void* mgl_ffty(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ register long i,nx=t->p[0],ny=t->p[1];
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ void *w = mgl_fft_alloc_thr(nx);
++#pragma omp for nowait
++ for(i=t->id;i<t->n;i+=mglNumThr)
++ mgl_fft(t->b+2*(i%nx)+2*nx*ny*(i/nx), nx, ny, t->v, w, t->p[3]);
++ mgl_fft_free_thr(w);
++ }
++ return 0;
++}
++MGL_NO_EXPORT void* mgl_fftz(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ register long i,nx=t->p[0],ny=t->p[1],nz=t->p[2];
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ void *w = mgl_fft_alloc_thr(nx);
++#pragma omp for nowait
++ for(i=t->id;i<t->n;i+=mglNumThr)
++ mgl_fft(t->b+2*i, nx*ny, nz, t->v, w, t->p[3]);
++ mgl_fft_free_thr(w);
++ }
++ return 0;
++}
++void MGL_EXPORT mgl_datac_fft(HADT d, const char *dir)
++{
++ if(!dir || *dir==0) return;
++ long nx = d->nx, ny = d->ny, nz = d->nz;
++ void *wt=0;
++ bool clear=false;
++ long par[4]={nx,ny,nz,strchr(dir,'i')!=0};
++#if MGL_USE_DOUBLE
++ double *a = (double *)(d->a);
++#else
++ double *a = new double[2*nx*ny*nz]; // manually convert to double
++#pragma omp parallel for
++ for(long i=0;i<nx*ny*nz;i++)
++ { a[2*i] = real(d->a[i]); a[2*i+1] = imag(d->a[i]); }
++#endif
++ if(strchr(dir,'x') && nx>1)
++ {
++ if(mgl_fft_data.wnx==nx) wt = mgl_fft_data.wtx;
++ else { clear = true; wt = mgl_fft_alloc(nx,0,0); }
++ mglStartThreadT(mgl_fftx,ny*nz,0,a,wt,0,par);
++ if(mgl_fft_data.wnx==0)
++ { clear = false; mgl_fft_data.wtx = wt; mgl_fft_data.wnx=nx; }
++ }
++ if(strchr(dir,'y') && ny>1)
++ {
++ if(mgl_fft_data.wny==ny) wt = mgl_fft_data.wty;
++ else { clear = true; wt = mgl_fft_alloc(ny,0,0); }
++ mglStartThreadT(mgl_ffty,nx*nz,0,a,wt,0,par);
++ if(mgl_fft_data.wny==0)
++ { clear = false; mgl_fft_data.wty = wt; mgl_fft_data.wny=ny; }
++ }
++ if(strchr(dir,'z') && nz>1)
++ {
++ if(mgl_fft_data.wnz==nz) wt = mgl_fft_data.wtz;
++ else { clear = true; wt = mgl_fft_alloc(nz,0,0); }
++ mglStartThreadT(mgl_fftz,nx*ny,0,a,wt,0,par);
++ if(mgl_fft_data.wnz==0)
++ { clear = false; mgl_fft_data.wtz = wt; mgl_fft_data.wnz=nz; }
++ }
++ if(clear) mgl_fft_free(wt,0,0);
++#if !MGL_USE_DOUBLE
++#pragma omp parallel for
++ for(long i=0;i<nx*ny*nz;i++) d->a[i] = dual(a[2*i], a[2*i+1]);
++ delete []a;
++#endif
++}
++//-----------------------------------------------------------------------------
++void MGL_EXPORT mgl_data_fourier(HMDT re, HMDT im, const char *dir)
++{
++ if(!dir || *dir==0) return;
++ long nx = re->nx, ny = re->ny, nz = re->nz;
++ if(nx*ny*nz != im->nx*im->ny*im->nz || !dir || dir[0]==0) return;
++ bool clear=false;
++ void *wt=0;
++ long par[4]={nx,ny,nz,strchr(dir,'i')!=0};
++ double *a = new double[2*nx*ny*nz];
++#pragma omp parallel for
++ for(long i=0;i<nx*ny*nz;i++)
++ { a[2*i] = re->a[i]; a[2*i+1] = im->a[i]; }
++ if(strchr(dir,'x') && nx>1)
++ {
++ if(mgl_fft_data.wnx==nx) wt = mgl_fft_data.wtx;
++ else { clear = true; wt = mgl_fft_alloc(nx,0,0); }
++ mglStartThreadT(mgl_fftx,ny*nz,0,a,wt,0,par);
++ if(mgl_fft_data.wnx==0)
++ { mgl_fft_data.wtx = wt; clear = false; mgl_fft_data.wnx=nx; }
++ }
++ if(strchr(dir,'y') && ny>1)
++ {
++ if(mgl_fft_data.wny==ny) wt = mgl_fft_data.wty;
++ else { clear = true; wt = mgl_fft_alloc(ny,0,0); }
++ mglStartThreadT(mgl_ffty,nx*nz,0,a,wt,0,par);
++ if(mgl_fft_data.wny==0)
++ { mgl_fft_data.wty = wt; clear = false; mgl_fft_data.wny=ny; }
++ }
++ if(strchr(dir,'z') && nz>1)
++ {
++ if(mgl_fft_data.wnz==nz) wt = mgl_fft_data.wtz;
++ else { clear = true; wt = mgl_fft_alloc(nz,0,0); }
++ mglStartThreadT(mgl_fftz,nx*ny,0,a,wt,0,par);
++ if(mgl_fft_data.wnz==0)
++ { mgl_fft_data.wtz = wt; clear = false; mgl_fft_data.wnz=nz; }
++ }
++ if(clear) { mgl_fft_free(wt,0,0); }
++#pragma omp parallel for
++ for(long i=0;i<nx*ny*nz;i++)
++ { re->a[i] = a[2*i]; im->a[i] = a[2*i+1]; }
++ delete []a;
++}
++//-----------------------------------------------------------------------------
++MGL_NO_EXPORT void* mgl_envx(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long nx=t->p[0];
++ mreal *a = (mreal*)t->a;
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ double *b = new double[2*nx];
++ void *w = mgl_fft_alloc_thr(nx);
++#pragma omp for nowait
++ for(long i=t->id;i<t->n;i+=mglNumThr)
++ {
++ for(long j=0;j<nx;j++) { b[2*j] = a[j+i*nx]; b[2*j+1] = 0; }
++ mgl_fft(b, 1, nx, t->v, w, false);
++ for(long j=0;j<nx;j++) { b[j] *= 2.; b[j+nx] = 0; }
++ mgl_fft(b, 1, nx, t->v, w, true);
++ for(long j=0;j<nx;j++) a[j+i*nx] = hypot(b[2*j], b[2*j+1]);
++ }
++ mgl_fft_free_thr(w); delete []b;
++ }
++ return 0;
++}
++MGL_NO_EXPORT void* mgl_envy(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long nx=t->p[0],ny=t->p[1];
++ mreal *a = (mreal*)t->a;
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ double *b = new double[2*ny];
++ void *w = mgl_fft_alloc_thr(ny);
++#pragma omp for nowait
++ for(long i=t->id;i<t->n;i+=mglNumThr)
++ {
++ for(long j=0;j<ny;j++) { b[2*j] = a[(i%nx)+nx*(j+ny*(i/nx))]; b[2*j+1] = 0; }
++ mgl_fft(b, 1, ny, t->v, t->w[t->id], false);
++ for(long j=0;j<ny;j++) { b[j] *= 2.; b[j+ny] = 0; }
++ mgl_fft(b, 1, ny, t->v, t->w[t->id], true);
++ for(long j=0;j<ny;j++) a[(i%nx)+nx*(j+ny*(i/nx))] = hypot(b[2*j], b[2*j+1]);
++ }
++ mgl_fft_free_thr(w); delete []b;
++ }
++ return 0;
++}
++MGL_NO_EXPORT void* mgl_envz(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long nx=t->p[0],ny=t->p[1],nz=t->p[2],k=nx*ny;
++ mreal *a = (mreal*)t->a;
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ double *b = new double[2*nz];
++ void *w = mgl_fft_alloc_thr(nz);
++#pragma omp for nowait
++ for(long i=t->id;i<t->n;i+=mglNumThr)
++ {
++ for(long j=0;j<nz;j++) { b[2*j] = a[j*k+i]; b[2*j+1] = 0; }
++ mgl_fft(b, 1, nz, t->v, t->w[t->id], false);
++ for(long j=0;j<nz;j++) { b[j] *= 2.; b[j+nz] = 0; }
++ mgl_fft(b, 1, nz, t->v, t->w[t->id], true);
++ for(long j=0;j<nz;j++) a[j*k+i] = hypot(b[2*j], b[2*j+1]);
++ }
++ mgl_fft_free_thr(w); delete []b;
++ }
++ return 0;
++}
++void MGL_EXPORT mgl_data_envelop(HMDT d, char dir)
++{
++ long nx=d->nx,ny=d->ny,nz=d->nz,par[3]={nx,ny,nz};
++ bool clear=false;
++ void *wt=0;
++ if(dir=='x' && nx>1)
++ {
++ if(mgl_fft_data.wnx==nx) wt = mgl_fft_data.wtx;
++ else { clear = true; wt = mgl_fft_alloc(nx,0,0); }
++ mglStartThreadT(mgl_envx,ny*nz,d->a,0,wt,0,par);
++ if(mgl_fft_data.wnx==0)
++ { mgl_fft_data.wtx = wt; clear = false; mgl_fft_data.wnx=nx; }
++ }
++ if(dir=='y' && ny>1)
++ {
++ if( mgl_fft_data.wny==ny) wt = mgl_fft_data.wty;
++ else { clear = true; wt = mgl_fft_alloc(ny,0,0); }
++ mglStartThreadT(mgl_envy,nx*nz,d->a,0,wt,0,par);
++ if(mgl_fft_data.wny==0)
++ { mgl_fft_data.wty = wt; clear = false; mgl_fft_data.wny=ny; }
++ }
++ if(dir=='z' && nz>1)
++ {
++ if(mgl_fft_data.wnz==nz) wt = mgl_fft_data.wtz;
++ else { clear = true; wt = mgl_fft_alloc(nz,0,0); }
++ mglStartThreadT(mgl_envz,nx*ny,d->a,0,wt,0,par);
++ if(mgl_fft_data.wnz==0)
++ { mgl_fft_data.wtz = wt; clear = false; mgl_fft_data.wnz=nz; }
++ }
++ if(clear) mgl_fft_free(wt,0,0);
++}
++//-----------------------------------------------------------------------------
++MGL_NO_EXPORT void* mgl_stfa1(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long mx=t->p[0],mz=t->p[2],dn=t->p[3],dd=dn/2,ny=t->p[4];
++ mreal *d = (mreal*)t->a;
++ HCDT re = (HCDT)t->re, im = (HCDT)t->im;
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ double *a = new double[4*dn], ff;
++ void *w = mgl_fft_alloc_thr(2*dn);
++#pragma omp for nowait
++ for(long ii=t->id;ii<t->n;ii+=mglNumThr)
++ {
++ register long i = ii%mx, j = ii/mx, i0;
++ for(long k=0;k<2*dn;k++)
++ {
++ i0 = k-dd+j*dn; ff = 1;
++ if(i0<0) i0=0; else if(i0>=ny) i0=ny-1;
++ if(k<dd)
++ { ff = 0.5*(k-dd/2.)/dd; ff=0.5+ff*(3-ff*ff); }
++ else if(k>=dn+dd)
++ { ff = 0.5*(k-3.5*dd)/dd; ff=0.5-ff*(3-ff*ff); }
++ a[2*k] = re->v(i,i0)*ff; a[2*k+1] = im->v(i,i0)*ff;
++ }
++ mgl_fft(a, 1, 2*dn, t->v, w, false);
++ for(long k=0;k<dd;k++)
++ {
++ i0 = i+mx*(j+mz*k);
++ d[i0+mx*mz*dd] = hypot(a[4*k],a[4*k+1])/dn;
++ d[i0] = hypot(a[4*k+2*dn],a[4*k+2*dn+1])/dn;
++ }
++ }
++ mgl_fft_free_thr(w); delete []a;
++ }
++ return 0;
++}
++MGL_NO_EXPORT void* mgl_stfa2(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long mx=t->p[0],my=t->p[1],dn=t->p[3],dd=dn/2,nx=t->p[4];
++ mreal *d = (mreal*)t->a;
++ HCDT re = (HCDT)t->re, im = (HCDT)t->im;
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ double *a = new double[4*dn], ff;
++ void *w = mgl_fft_alloc_thr(2*dn);
++#pragma omp for nowait
++ for(long ii=t->id;ii<t->n;ii+=mglNumThr)
++ {
++ register long i = ii%my, j = ii/my, i0;
++ for(long k=0;k<2*dn;k++)
++ {
++ i0 = k-dd+i*dn; ff = 1;
++ if(i0<0) i0=0; else if(i0>=nx) i0=nx-1;
++ if(k<dd)
++ { ff = 0.5*(k-dd/2.)/dd; ff=0.5+ff*(3-ff*ff); }
++ else if(k>=3*dd)
++ { ff = 0.5*(k-3.5*dd)/dd; ff=0.5-ff*(3-ff*ff); }
++ a[2*k] = re->v(i0,j)*ff; a[2*k+1] = im->v(i0,j)*ff;
++ }
++ mgl_fft(a, 1, 2*dn, t->v, w, false);
++ for(long k=0;k<dd;k++)
++ {
++ i0 = i+my*(k+mx*j);
++ d[i0+dd*my] = hypot(a[4*k],a[4*k+1])/dn;
++ d[i0] = hypot(a[4*k+2*dn],a[4*k+2*dn+1])/dn;
++ }
++ }
++ mgl_fft_free_thr(w); delete []a;
++ }
++ return 0;
++}
++HMDT MGL_EXPORT mgl_data_stfa(HCDT re, HCDT im, long dn, char dir)
++{
++ if(dn<2) return 0;
++ dn = 2*(dn/2);
++ long nx = re->GetNx(), ny = re->GetNy();
++ if(nx*ny!=im->GetNx()*im->GetNy()) return 0;
++ void *wt = mgl_fft_alloc(2*dn,0,0);
++ long mx,my,mz;
++ mglData *d=new mglData;
++ if(dir=='y')
++ {
++ mx = nx; my = dn; mz = ny/dn;
++ mgl_data_create(d, mx, mz, my);
++ long par[5]={mx,my,mz,dn,ny};
++ mglStartThreadT(mgl_stfa1,mx*mz,d->a,0,wt,0,par,re,im);
++ }
++ else
++ {
++ mx = dn; my = nx/dn; mz = ny;
++ mgl_data_create(d, my, mx, mz);
++ long par[5]={mx,my,mz,dn,nx};
++ mglStartThreadT(mgl_stfa2,my*mz,d->a,0,wt,0,par,re,im);
++ }
++ mgl_fft_free(wt,0,0);
++ return d;
++}
++//-----------------------------------------------------------------------------
++MGL_NO_EXPORT void* mgl_sinx(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long nx=t->p[0];
++ mreal *a = (mreal*)t->a;
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ double *b = new double[2*nx], f=sqrt(2./nx);
++ void *w = mgl_fft_alloc_thr(nx);
++#pragma omp for nowait
++ for(long i=t->id;i<t->n;i+=mglNumThr)
++ {
++ register long k = i*nx; memset(b,0,2*nx*sizeof(double));
++ for(long j=1;j<nx;j++) b[2*j]=sin(M_PI*j/nx)*(a[j+k]+a[nx-j+k])+(a[j+k]-a[nx-j+k])*0.5;
++ mgl_fft(b,1,nx,t->v,w,false);
++ a[k]=0; a[k+1]=b[0]*f/2; // fill sinfft
++ for(long j=1;j<nx/2;j++)
++ {
++ a[k+2*j] = -b[2*j+1]*f;
++ a[k+2*j+1] = a[k+2*j-1]+b[2*j]*f;
++ }
++ if(nx%2) a[nx-1] = -b[nx]*f;
++ }
++ mgl_fft_free_thr(w); delete []b;
++ }
++ return 0;
++}
++MGL_NO_EXPORT void* mgl_siny(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long nx=t->p[0],ny=t->p[1];
++ mreal *a = (mreal*)t->a;
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ double *b = new double[2*ny], f=sqrt(2./ny);
++ void *w = mgl_fft_alloc_thr(ny);
++#pragma omp for nowait
++ for(long ii=t->id;ii<t->n;ii+=mglNumThr)
++ {
++ register long i = ii%nx, k = ii/nx; memset(b,0,2*ny*sizeof(double));
++ for(long j=1;j<ny;j++) b[2*j]=sin(M_PI*j/ny)*(a[i+nx*(ny*k+j)]+a[i+nx*(ny*k+ny-j)])+(a[i+nx*(ny*k+j)]-a[i+nx*(ny*k+ny-j)])*0.5;
++ mgl_fft(b,1,ny,t->v,w,false);
++ a[i+nx*ny*k]=0; a[i+nx*(ny*k+1)]=b[0]*f/2; // fill sinfft
++ for(long j=1;j<ny/2;j++)
++ {
++ a[i+nx*(ny*k+2*j)] = -b[2*j+1]*f;
++ a[i+nx*(ny*k+2*j+1)] = a[i+nx*(ny*k+2*j-1)]+b[2*j]*f;
++ }
++ if(ny%2) a[i+nx*(ny*k+ny-1)] = -b[ny]*f;
++ }
++ mgl_fft_free_thr(w); delete []b;
++ }
++ return 0;
++}
++MGL_NO_EXPORT void* mgl_sinz(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long nx=t->p[0],ny=t->p[1],nz=t->p[2],k=nx*ny;
++ mreal *a = (mreal*)t->a;
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ double *b = new double[2*nz], f=sqrt(2./nz);
++ void *w = mgl_fft_alloc_thr(nz);
++#pragma omp for nowait
++ for(long i=t->id;i<t->n;i+=mglNumThr)
++ {
++ memset(b,0,2*nz*sizeof(double));
++ for(long j=1;j<nz;j++) b[2*j]=sin(M_PI*j/nz)*(a[i+k*j]+a[i+k*(nz-j)])+(a[i+k*j]-a[i+k*(nz-j)])*0.5;
++ mgl_fft(b,1,nz,t->v,w,false);
++ a[i]=0; a[i+k]=b[0]*f/2; // fill sinfft
++ for(long j=1;j<nz/2;j++)
++ {
++ a[i+k*2*j] = -b[2*j+1]*f;
++ a[i+k*(2*j+1)] = a[i+k*(2*j-1)]+b[2*j]*f;
++ }
++ if(nz%2) a[i+k*nz-k] = -b[nz]*f;
++ }
++ mgl_fft_free_thr(w); delete []b;
++ }
++ return 0;
++}
++void MGL_EXPORT mgl_data_sinfft(HMDT d, const char *dir) // use DST-1
++{
++ if(!dir || *dir==0) return;
++ bool clear=false;
++ void *wt=0;
++ long nx=d->nx, ny=d->ny, nz=d->nz, par[3]={nx,ny,nz};
++ if(strchr(dir,'x') && nx>1)
++ {
++ if(mgl_fft_data.wnx==nx) wt = mgl_fft_data.wtx;
++ else { clear = true; wt = mgl_fft_alloc(nx,0,0); }
++ mglStartThreadT(mgl_sinx,ny*nz,d->a,0,wt,0,par);
++ if(mgl_fft_data.wnx==0)
++ { mgl_fft_data.wtx = wt; clear = false; mgl_fft_data.wnx=nx; }
++ }
++ if(strchr(dir,'y') && ny>1)
++ {
++ if(mgl_fft_data.wny==ny) wt = mgl_fft_data.wty;
++ else { clear = true; wt = mgl_fft_alloc(ny,0,0); }
++ mglStartThreadT(mgl_siny,nx*nz,d->a,0,wt,0,par);
++ if(mgl_fft_data.wny==0)
++ { mgl_fft_data.wty = wt; clear = false; mgl_fft_data.wny=ny; }
++ }
++ if(strchr(dir,'z') && nz>1)
++ {
++ if(mgl_fft_data.wnz==nz) wt = mgl_fft_data.wtz;
++ else { clear = true; wt = mgl_fft_alloc(nz,0,0); }
++ mglStartThreadT(mgl_sinz,nx*ny,d->a,0,wt,0,par);
++ if(mgl_fft_data.wnz==0)
++ { mgl_fft_data.wtz = wt; clear = false; mgl_fft_data.wnz=nz; }
++ }
++ if(clear) mgl_fft_free(wt,0,0);
++}
++//-----------------------------------------------------------------------------
++MGL_NO_EXPORT void* mgl_cosx(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long nx=t->p[0],nn=nx-1;
++ mreal *a = (mreal*)t->a;
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ double *b = new double[2*nx], f=sqrt(2./nn);
++ void *w = mgl_fft_alloc_thr(nn);
++#pragma omp for nowait
++ for(long i=t->id;i<t->n;i+=mglNumThr)
++ {
++ register long k = i*nx; memset(b,0,2*nx*sizeof(double));
++ for(long j=0;j<nn;j++) b[2*j]=(a[j+k]+a[nn-j+k])*0.5-sin(M_PI*j/nn)*(a[j+k]-a[nn-j+k]);
++ mgl_fft(b,1,nn,t->v,w,false);
++ double f1=0.5*(a[k]-a[nn+k]), s=-1;
++ a[nn+k]=0.5*(a[k]+a[nn+k]*(nn%2?-1:1));
++ for(long j=1;j<nn;j++)
++ {
++ f1 += a[j+k]*cos(M_PI*j/nn);
++ a[nn+k] += a[j+k]*s; s = -s;
++ }
++ a[k]=b[0]*f; a[1+k]=f1*f; a[nn+k]*=f; // fill cosfft
++ for(long j=1;j<nn/2;j++)
++ {
++ a[2*j+k] = b[2*j]*f;
++ a[2*j+1+k] = a[2*j-1+k]-b[2*j+1]*f;
++ }
++ if(nn%2) a[nn-1+k] = b[nn-1]*f;
++ }
++ mgl_fft_free_thr(w);
++ delete []b;
++ }
++ return 0;
++}
++MGL_NO_EXPORT void* mgl_cosy(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long nx=t->p[0],ny=t->p[1],nn=ny-1;
++ mreal *a = (mreal*)t->a;
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ double *b = new double[2*ny], f=sqrt(2./nn);
++ void *w = mgl_fft_alloc_thr(nn);
++#pragma omp for nowait
++ for(long ii=t->id;ii<t->n;ii+=mglNumThr)
++ {
++ register long i = ii%nx, k = ii/nx; memset(b,0,2*ny*sizeof(double));
++ for(long j=0;j<nn;j++) b[2*j]=(a[i+nx*(ny*k+j)]+a[i+nx*(ny*k+nn-j)])*0.5-sin(M_PI*j/nn)*(a[i+nx*(ny*k+j)]-a[i+nx*(ny*k+nn-j)]);
++ mgl_fft(b,1,nn,t->v,w,false);
++ double f1=0.5*(a[i+nx*ny*k]-a[i+nx*(ny*k+nn)]), s=-1;
++ a[i+nx*(ny*k+nn)]=0.5*(a[i+nx*ny*k]+a[i+nx*(ny*k+nn)]*(nn%2?-1:1));
++ for(long j=1;j<nn;j++)
++ {
++ f1 += a[i+nx*(ny*k+j)]*cos(M_PI*j/nn);
++ a[i+nx*(ny*k+nn)] += a[i+nx*(ny*k+j)]*s; s = -s;
++ }
++ a[i+nx*ny*k]=b[0]*f; a[i+nx*(ny*k+1)]=f1*f; a[i+nx*(ny*k+nn)]*=f; // fill cosfft
++ for(long j=1;j<nn/2;j++)
++ {
++ a[i+nx*(ny*k+2*j)] = b[2*j]*f;
++ a[i+nx*(ny*k+2*j+1)] = a[i+nx*(ny*k+2*j-1)]-b[2*j+1]*f;
++ }
++ if(nn%2) a[i+nx*(ny*k+nn-1)] = b[nn-1]*f;
++ }
++ mgl_fft_free_thr(w);
++ delete []b;
++ }
++ return 0;
++}
++MGL_NO_EXPORT void* mgl_cosz(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long nx=t->p[0],ny=t->p[1],nz=t->p[2],k=nx*ny,nn=nz-1;
++ mreal *a = (mreal*)t->a;
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ double *b = new double[2*nz], f=sqrt(2./nn);
++ void *w = mgl_fft_alloc_thr(nn);
++#pragma omp for nowait
++ for(long i=t->id;i<t->n;i+=mglNumThr)
++ {
++ memset(b,0,2*nz*sizeof(double));
++ for(long j=0;j<nn;j++) b[2*j]=(a[i+k*j]+a[i+k*(nn-j)])*0.5-sin(M_PI*j/nn)*(a[i+k*j]-a[i+k*(nn-j)]);
++ mgl_fft(b,1,nn,t->v,w,false);
++ double f1=0.5*(a[i]-a[i+k*nn]), s=-1;
++ a[i+k*nn]=0.5*(a[i]+a[i+k*nn]*(nn%2?-1:1));
++ for(long j=1;j<nn;j++)
++ {
++ f1 += a[i+k*j]*cos(M_PI*j/nn);
++ a[i+k*nn] += a[i+k*j]*s; s = -s;
++ }
++ a[i]=b[0]*f; a[i+k]=f1*f; a[i+k*nn]*=f; // fill cosfft
++ for(long j=1;j<nn/2;j++)
++ {
++ a[i+k*2*j] = b[2*j]*f;
++ a[i+k*2*j+k] = a[i+k*2*j-k]-b[2*j+1]*f;
++ }
++ if(nn%2) a[i+k*nn-k] = b[nn-1]*f;
++ }
++ mgl_fft_free_thr(w);
++ delete []b;
++ }
++ return 0;
++}
++void MGL_EXPORT mgl_data_cosfft(HMDT d, const char *dir)
++{
++ if(!dir || *dir==0) return;
++ bool clear=false;
++ void *wt=0;
++ long nx=d->nx, ny=d->ny, nz=d->nz, par[3]={nx,ny,nz};
++ if(strchr(dir,'x') && nx>1)
++ {
++ if(mgl_fft_data.wnx==nx-1) wt = mgl_fft_data.wtx;
++ else { clear = true; wt = mgl_fft_alloc(nx-1,0,0); }
++ mglStartThreadT(mgl_cosx,ny*nz,d->a,0,wt,0,par);
++ if(mgl_fft_data.wnx==0)
++ { mgl_fft_data.wtx = wt; clear = false; mgl_fft_data.wnx=nx-1; }
++ }
++ if(strchr(dir,'y') && ny>1)
++ {
++ if(mgl_fft_data.wny==ny-1) wt = mgl_fft_data.wty;
++ else { clear = true; wt = mgl_fft_alloc(ny-1,0,0); }
++ mglStartThreadT(mgl_cosy,nx*nz,d->a,0,wt,0,par);
++ if(mgl_fft_data.wny==0)
++ { mgl_fft_data.wty = wt; clear = false; mgl_fft_data.wny=ny-1; }
++ }
++ if(strchr(dir,'z') && nz>1)
++ {
++ if(mgl_fft_data.wnz==nz-1) wt = mgl_fft_data.wtz;
++ else { clear = true; wt = mgl_fft_alloc(nz-1,0,0); }
++ mglStartThreadT(mgl_cosz,nx*ny,d->a,0,wt,0,par);
++ if(mgl_fft_data.wnz==0)
++ { mgl_fft_data.wtz = wt; clear = false; mgl_fft_data.wnz=nz-1; }
++ }
++ if(clear) mgl_fft_free(wt,0,0);
++}
++//-----------------------------------------------------------------------------
++HMDT MGL_EXPORT mgl_transform_a(HCDT am, HCDT ph, const char *tr)
++{
++ long nx = am->GetNx(), ny = am->GetNy(), nz = am->GetNz();
++ if(nx*ny*nz != ph->GetNx()*ph->GetNy()*ph->GetNz() || !tr || tr[0]==0)
++ return 0;
++ mglData re(nx,ny,nz), im(nx,ny,nz);
++ const mglData *da=dynamic_cast<const mglData *>(am);
++ const mglData *dp=dynamic_cast<const mglData *>(ph);
++ if(da && dp)
++#pragma omp parallel for
++ for(long i=0;i<nx*ny*nz;i++)
++ { re.a[i] = da->a[i]*cos(dp->a[i]);
++ im.a[i] = da->a[i]*sin(dp->a[i]); }
++ else
++#pragma omp parallel for
++ for(long i=0;i<nx*ny*nz;i++)
++ { re.a[i] = am->vthr(i)*cos(ph->vthr(i));
++ im.a[i] = am->vthr(i)*sin(ph->vthr(i)); }
++ return mgl_transform(&re, &im, tr);
++}
++//-----------------------------------------------------------------------------
++HMDT MGL_EXPORT mgl_transform(HCDT re, HCDT im, const char *tr)
++{
++ if(!tr || *tr==0) return 0;
++ long nx = re->GetNx(), ny = re->GetNy(), nz = re->GetNz();
++ if(nx*ny*nz != im->GetNx()*im->GetNy()*im->GetNz() || !tr || tr[0]==0)
++ return 0;
++ mglData rr(re),ii(im);
++ if(strchr(tr,'i') && strchr(tr,'f')) // general case
++ {
++ if(tr[0]=='f') mgl_data_fourier(&rr,&ii,"x");
++ if(tr[0]=='i') mgl_data_fourier(&rr,&ii,"xi");
++ if(tr[1]=='f') mgl_data_fourier(&rr,&ii,"y");
++ if(tr[1]=='i') mgl_data_fourier(&rr,&ii,"yi");
++ if(tr[2]=='f') mgl_data_fourier(&rr,&ii,"z");
++ if(tr[2]=='i') mgl_data_fourier(&rr,&ii,"zi");
++ }
++ else if(strchr(tr,'f')) // do Fourier only once for speeding up
++ {
++ char str[4] = " ";
++ if(tr[0]=='f') str[0]='x';
++ if(tr[1]=='f') str[1]='y';
++ if(tr[2]=='f') str[2]='z';
++ mgl_data_fourier(&rr,&ii,str);
++ }
++ else if(strchr(tr,'i')) // do Fourier only once for speeding up
++ {
++ char str[5] = " i";
++ if(tr[0]=='i') str[0]='x';
++ if(tr[1]=='i') str[1]='y';
++ if(tr[2]=='i') str[2]='z';
++ mgl_data_fourier(&rr,&ii,str);
++ }
++ else if(strchr(tr,'s')) // do Fourier only once for speeding up
++ {
++ if(tr[0]=='s') { rr.SinFFT("x"); ii.SinFFT("x"); }
++ if(tr[1]=='s') { rr.SinFFT("y"); ii.SinFFT("y"); }
++ if(tr[2]=='s') { rr.SinFFT("z"); ii.SinFFT("z"); }
++ }
++ else if(strchr(tr,'c')) // do Fourier only once for speeding up
++ {
++ if(tr[0]=='c') { rr.CosFFT("x"); ii.CosFFT("x"); }
++ if(tr[1]=='c') { rr.CosFFT("y"); ii.CosFFT("y"); }
++ if(tr[2]=='c') { rr.CosFFT("z"); ii.CosFFT("z"); }
++ }
++ else if(strchr(tr,'h')) // do Fourier only once for speeding up
++ {
++ if(tr[0]=='h') { rr.Hankel("x"); ii.Hankel("x"); }
++ if(tr[1]=='h') { rr.Hankel("y"); ii.Hankel("y"); }
++ if(tr[2]=='h') { rr.Hankel("z"); ii.Hankel("z"); }
++ }
++ mglData *d = new mglData(nx, ny, nz);
++#pragma omp parallel for
++ for(long i=0;i<nx*ny*nz;i++) d->a[i] = hypot(rr.a[i],ii.a[i]);
++ return d;
++}
++//-----------------------------------------------------------------------------
++uintptr_t MGL_EXPORT mgl_transform_a_(uintptr_t *am, uintptr_t *ph, const char *tr, int l)
++{ char *s=new char[l+1]; memcpy(s,tr,l); s[l]=0;
++ uintptr_t res = uintptr_t(mgl_transform_a(_DA_(am),_DA_(ph),s));
++ delete []s; return res; }
++uintptr_t MGL_EXPORT mgl_transform_(uintptr_t *re, uintptr_t *im, const char *tr, int l)
++{ char *s=new char[l+1]; memcpy(s,tr,l); s[l]=0;
++ uintptr_t res = uintptr_t(mgl_transform(_DA_(re),_DA_(im),s));
++ delete []s; return res; }
++//-----------------------------------------------------------------------------
++void MGL_EXPORT mgl_data_envelop_(uintptr_t *d, const char *dir, int)
++{ mgl_data_envelop(_DT_,*dir); }
++//-----------------------------------------------------------------------------
++#if MGL_HAVE_GSL
++MGL_NO_EXPORT void* mgl_chnkx(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long nx=t->p[0];
++ dual *a = (dual*)t->a;
++ const gsl_dht *dht = (const gsl_dht*)t->v;
++ double mm = gsl_sf_bessel_zero_J0(nx+1);
++
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ double *b = new double[3*nx];
++#pragma omp for nowait
++ for(long i=t->id;i<t->n;i+=mglNumThr)
++ {
++ for(long j=0;j<nx;j++) b[j] = real(a[j+nx*i]);
++ gsl_dht_apply(dht,b,b+nx);
++ for(long j=0;j<nx;j++) b[j] = imag(a[j+nx*i]);
++ gsl_dht_apply(dht,b,b+2*nx);
++ for(long j=0;j<nx;j++) a[j+nx*i] = dual(b[j+nx]*mm,b[j+2*nx]*mm);
++ }
++ delete []b;
++ }
++ return 0;
++}
++MGL_NO_EXPORT void* mgl_chnky(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long nx=t->p[0],ny=t->p[1];
++ dual *a = (dual*)t->a;
++ const gsl_dht *dht = (const gsl_dht*)t->v;
++ double mm = gsl_sf_bessel_zero_J0(ny+1);
++
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ double *b = new double[3*ny];
++#pragma omp for nowait
++ for(long ii=t->id;ii<t->n;ii+=mglNumThr)
++ {
++ register long i = ii%nx, k = ii/nx;
++ for(long j=0;j<ny;j++) b[j] = real(a[i+nx*(j+ny*k)]);
++ gsl_dht_apply(dht,b,b+ny);
++ for(long j=0;j<ny;j++) b[j] = imag(a[i+nx*(j+ny*k)]);
++ gsl_dht_apply(dht,b,b+2*ny);
++ for(long j=0;j<ny;j++) a[i+nx*(j+ny*k)] = dual(b[j+ny]*mm,b[j+2*ny]*mm);
++ }
++ delete []b;
++ }
++ return 0;
++}
++MGL_NO_EXPORT void* mgl_chnkz(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long k=t->p[0]*t->p[1],nz=t->p[2];
++ dual *a = (dual*)t->a;
++ const gsl_dht *dht = (const gsl_dht*)t->v;
++ double mm = gsl_sf_bessel_zero_J0(nz+1);
++
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ double *b = new double[3*nz];
++#pragma omp for nowait
++ for(long i=t->id;i<t->n;i+=mglNumThr)
++ {
++ for(long j=0;j<nz;j++) b[j] = real(a[i+j*k]);
++ gsl_dht_apply(dht,b,b+nz);
++ for(long j=0;j<nz;j++) b[j] = imag(a[i+j*k]);
++ gsl_dht_apply(dht,b,b+2*nz);
++ for(long j=0;j<nz;j++) a[i+j*k] = dual(b[j+nz]*mm,b[j+2*nz]*mm);
++ }
++ delete []b;
++ }
++ return 0;
++}
++void MGL_EXPORT mgl_datac_hankel(HADT d, const char *dir)
++{
++ if(!dir || *dir==0) return;
++ gsl_dht *dht=0;
++ bool clear = false;
++ long nx=d->nx, ny=d->ny, nz=d->nz;
++ long par[3]={nx,ny,nz};
++ if(strchr(dir,'x') && nx>1)
++ {
++ if(mgl_fft_data.hnx==nx) dht = (gsl_dht *)mgl_fft_data.htx;
++ else { dht = gsl_dht_new(nx,0,1); clear = true; }
++ mglStartThreadT(mgl_chnkx,ny*nz,d->a,0,dht,0,par);
++ if(mgl_fft_data.hnx==0)
++ { mgl_fft_data.htx = dht; clear = false; mgl_fft_data.hnx=nx; }
++ }
++ if(strchr(dir,'y') && ny>1)
++ {
++ if(mgl_fft_data.hny==ny) dht = (gsl_dht *)mgl_fft_data.hty;
++ else { dht = gsl_dht_new(ny,0,1); clear = true; }
++ mglStartThreadT(mgl_chnky,nx*nz,d->a,0,dht,0,par);
++ if(mgl_fft_data.hny==0)
++ { mgl_fft_data.hty = dht; clear = false; mgl_fft_data.hny=ny; }
++ }
++ if(strchr(dir,'z') && nz>1)
++ {
++ if(mgl_fft_data.hnz==nz) dht = (gsl_dht *)mgl_fft_data.htz;
++ else { dht = gsl_dht_new(nz,0,1); clear = true; }
++ mglStartThreadT(mgl_chnkz,nx*ny,d->a,0,dht,0,par);
++ if(mgl_fft_data.hnz==0)
++ { mgl_fft_data.htz = dht; clear = false; mgl_fft_data.hnz=nz; }
++ }
++ if(clear) gsl_dht_free(dht);
++}
++#else
++void MGL_EXPORT mgl_datac_hankel(HADT , const char *){}
++#endif
++void MGL_EXPORT mgl_datac_hankel_(uintptr_t *d, const char *dir,int l)
++{ char *s=new char[l+1]; memcpy(s,dir,l); s[l]=0;
++ mgl_datac_hankel(_DC_,s); delete []s; }
++//-----------------------------------------------------------------------------
++#if MGL_HAVE_GSL
++MGL_NO_EXPORT void* mgl_hnkx(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long nx=t->p[0];
++ mreal *a = (mreal*)t->a;
++ const gsl_dht *dht = (const gsl_dht*)t->v;
++ double mm = gsl_sf_bessel_zero_J0(nx+1);
++
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ double *b = new double[2*nx];
++#pragma omp for nowait
++ for(long i=t->id;i<t->n;i+=mglNumThr)
++ {
++ for(long j=0;j<nx;j++) b[j] = a[j+nx*i];
++ gsl_dht_apply(dht,b,b+nx);
++ for(long j=0;j<nx;j++) a[j+nx*i] = b[j+nx]*mm;
++ }
++ delete []b;
++ }
++ return 0;
++}
++MGL_NO_EXPORT void* mgl_hnky(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long nx=t->p[0],ny=t->p[1];
++ mreal *a = (mreal*)t->a;
++ const gsl_dht *dht = (const gsl_dht*)t->v;
++ double mm = gsl_sf_bessel_zero_J0(ny+1);
++
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ double *b = new double[2*ny];
++#pragma omp for nowait
++ for(long ii=t->id;ii<t->n;ii+=mglNumThr)
++ {
++ register long i = ii%nx, k = ii/nx;
++ for(long j=0;j<ny;j++) b[j] = a[i+nx*(j+ny*k)];
++ gsl_dht_apply(dht,b,b+ny);
++ for(long j=0;j<ny;j++)a[i+nx*(j+ny*k)] = b[j+ny]*mm;
++ }
++ delete []b;
++ }
++ return 0;
++}
++MGL_NO_EXPORT void* mgl_hnkz(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long k=t->p[0]*t->p[1],nz=t->p[2];
++ mreal *a = (mreal*)t->a;
++ const gsl_dht *dht = (const gsl_dht*)t->v;
++ double mm = gsl_sf_bessel_zero_J0(nz+1);
++
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ double *b = new double[2*nz];
++#pragma omp for nowait
++ for(long i=t->id;i<t->n;i+=mglNumThr)
++ {
++ for(long j=0;j<nz;j++) b[j] = a[i+j*k];
++ gsl_dht_apply(dht,b,b+nz);
++ for(long j=0;j<nz;j++) a[i+j*k] = b[j+nz]*mm;
++ }
++ delete []b;
++ }
++ return 0;
++}
++void MGL_EXPORT mgl_data_hankel(HMDT d, const char *dir)
++{
++ if(!dir || *dir==0) return;
++ bool clear = false;
++ gsl_dht *dht=0;
++ long nx=d->nx, ny=d->ny, nz=d->nz;
++ long par[3]={nx,ny,nz};
++ if(strchr(dir,'x') && nx>1)
++ {
++ if(mgl_fft_data.hnx==nx) dht = (gsl_dht *)mgl_fft_data.htx;
++ else { dht = gsl_dht_new(nx,0,1); clear = true; }
++ mglStartThreadT(mgl_hnkx,ny*nz,d->a,0,dht,0,par);
++ if(mgl_fft_data.hnx==0)
++ { mgl_fft_data.htx = dht; clear = false; mgl_fft_data.hnx=nx; }
++ }
++ if(strchr(dir,'y') && ny>1)
++ {
++ if(mgl_fft_data.hny==ny) dht = (gsl_dht *)mgl_fft_data.hty;
++ else { dht = gsl_dht_new(ny,0,1); clear = true; }
++ mglStartThreadT(mgl_hnky,nx*nz,d->a,0,dht,0,par);
++ if(mgl_fft_data.hny==0)
++ { mgl_fft_data.hty = dht; clear = false; mgl_fft_data.hny=ny; }
++ }
++ if(strchr(dir,'z') && nz>1)
++ {
++ if(mgl_fft_data.hnz==nz) dht = (gsl_dht *)mgl_fft_data.htz;
++ else { dht = gsl_dht_new(nz,0,1); clear = true; }
++ mglStartThreadT(mgl_hnkz,nx*ny,d->a,0,dht,0,par);
++ if(mgl_fft_data.hnz==0)
++ { mgl_fft_data.htz = dht; clear = false; mgl_fft_data.hnz=nz; }
++ }
++ if(clear) gsl_dht_free(dht);
++}
++#else
++void MGL_EXPORT mgl_data_hankel(HMDT , const char *){}
++#endif
++void MGL_EXPORT mgl_data_hankel_(uintptr_t *d, const char *dir,int l)
++{ char *s=new char[l+1]; memcpy(s,dir,l); s[l]=0;
++ mgl_data_hankel(_DT_,s); delete []s; }
++//-----------------------------------------------------------------------------
++void MGL_EXPORT mgl_data_fill_sample(HMDT d, const char *how)
++{
++ if(!how || *how==0) return;
++ bool xx = strchr(how,'x');
++ long n=d->nx;
++ mreal *aa=d->a;
++ if(strchr(how,'h')) // Hankel
++ {
++#if MGL_HAVE_GSL
++ gsl_dht *dht = gsl_dht_new(n,0,1);
++#pragma omp parallel for
++ for(long i=0;i<n;i++)
++ aa[i] = xx ? gsl_dht_x_sample(dht, i) : gsl_dht_k_sample(dht, i);
++ gsl_dht_free(dht);
++#endif
++ }
++ else // Fourier
++ {
++ if(xx)
++#pragma omp parallel for
++ for(long i=0;i<n;i++) aa[i] = mreal(2*i-n)/n;
++ else
++#pragma omp parallel for
++ for(long i=0;i<n;i++) aa[i] = M_PI*(i<n/2 ? i:i-n);
++ }
++#pragma omp parallel for
++ for(long i=1;i<d->ny*d->nz;i++) memcpy(aa+i*n,aa,n*sizeof(mreal));
++}
++void MGL_EXPORT mgl_data_fill_sample_(uintptr_t *d, const char *how,int l)
++{ char *s=new char[l+1]; memcpy(s,how,l); s[l]=0;
++ mgl_data_fill_sample(_DT_,s); delete []s; }
++//-----------------------------------------------------------------------------
++void MGL_EXPORT mgl_datac_fft_(uintptr_t *d, const char *dir, int l)
++{ char *s=new char[l+1]; memcpy(s,dir,l); s[l]=0;
++ mgl_datac_fft(_DC_,s); delete []s; }
++void MGL_EXPORT mgl_data_fourier_(uintptr_t *re, uintptr_t *im, const char *dir, int l)
++{ char *s=new char[l+1]; memcpy(s,dir,l); s[l]=0;
++ mgl_data_fourier(_DM_(re),_DM_(im),s); delete []s; }
++uintptr_t MGL_EXPORT mgl_data_stfa_(uintptr_t *re, uintptr_t *im, int *dn, char *dir, int)
++{ return uintptr_t(mgl_data_stfa(_DA_(re),_DA_(im),*dn,*dir)); }
++void MGL_EXPORT mgl_data_cosfft_(uintptr_t *d, const char *dir,int l)
++{ char *s=new char[l+1]; memcpy(s,dir,l); s[l]=0;
++ mgl_data_cosfft(_DT_,s); delete []s; }
++void MGL_EXPORT mgl_data_sinfft_(uintptr_t *d, const char *dir,int l)
++{ char *s=new char[l+1]; memcpy(s,dir,l); s[l]=0;
++ mgl_data_sinfft(_DT_,s); delete []s; }
++//-----------------------------------------------------------------------------
++MGL_NO_EXPORT void* mgl_corx(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long nx=t->p[0];
++ double *a = (double *)t->a;
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ void *w = mgl_fft_alloc_thr(nx);
++#pragma omp for nowait
++ for(long i=t->id;i<t->n;i+=mglNumThr)
++ {
++ mgl_fft(t->b+2*nx*i, 1, nx, t->v, w, false);
++ mgl_fft(a+2*nx*i, 1, nx, t->v, w, false);
++ for(long j=0;j<nx;j++)
++ {
++ register long ii = 2*j+2*nx*i;
++ register double re = t->b[ii], im = t->b[ii+1];
++ t->b[ii] = re*a[ii] + im*a[ii+1];
++ t->b[ii+1] = im*a[ii] - re*a[ii+1];
++ }
++ mgl_fft(t->b+2*nx*i, 1, nx, t->v, w, true);
++ }
++ mgl_fft_free_thr(w);
++ }
++ return 0;
++}
++MGL_NO_EXPORT void* mgl_cory(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long nx=t->p[0],ny=t->p[1];
++ double *a = (double *)t->a;
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ void *w = mgl_fft_alloc_thr(ny);
++#pragma omp for nowait
++ for(long i=t->id;i<t->n;i+=mglNumThr)
++ {
++ register long k = 2*(i%nx)+2*nx*ny*(i/nx);
++ mgl_fft(t->b+k, nx, ny, t->v, w, false);
++ mgl_fft(a+k, nx, ny, t->v, w, false);
++ for(long j=0;j<ny;j++)
++ {
++ register long ii = 2*nx*j+k;
++ register double re = t->b[ii], im = t->b[ii+1];
++ t->b[ii] = re*a[ii] + im*a[ii+1];
++ t->b[ii+1] = im*a[ii] - re*a[ii+1];
++ }
++ mgl_fft(t->b+k, nx, ny, t->v, w, true);
++ }
++ mgl_fft_free_thr(w);
++ }
++ return 0;
++}
++MGL_NO_EXPORT void* mgl_corz(void *par)
++{
++ mglThreadT *t=(mglThreadT *)par;
++ long nx=t->p[0],ny=t->p[1],nz=t->p[2];
++ double *a = (double *)t->a;
++#if !MGL_HAVE_PTHREAD
++#pragma omp parallel
++#endif
++ {
++ void *w = mgl_fft_alloc_thr(nz);
++#pragma omp for nowait
++ for(long i=t->id;i<t->n;i+=mglNumThr)
++ {
++ mgl_fft(t->b+2*i, nx*ny, nz, t->v, w, false);
++ mgl_fft(a+2*i, nx*ny, nz, t->v, w, false);
++ for(long j=0;j<nz;j++)
++ {
++ register long ii = 2*nx*ny*j+2*i;
++ register double re = t->b[ii], im = t->b[ii+1];
++ t->b[ii] = re*a[ii] + im*a[ii+1];
++ t->b[ii+1] = im*a[ii] - re*a[ii+1];
++ }
++ mgl_fft(t->b+2*i, nx*ny, nz, t->v, w, true);
++ }
++ mgl_fft_free_thr(w);
++ }
++ return 0;
++}
++MGL_NO_EXPORT double *mgl_d_correl(HCDT d1, HCDT d2, const char *dir)
++{
++ if(!dir || *dir==0) return 0;
++ long nx = d1->GetNx(), ny = d1->GetNy(), nz = d1->GetNz(), nn=nx*ny*nz;
++ if(nx*ny*nz!=d2->GetNN()) return 0;
++ void *wt=0;
++ bool clear=false;
++ long par[3]={nx,ny,nz};
++
++ double *a = new double[2*nn]; memset(a,0,2*nn*sizeof(double));
++ double *b = new double[2*nn]; memset(b,0,2*nn*sizeof(double));
++ const mglDataC *dd1 = dynamic_cast<const mglDataC *>(d1);
++ const mglDataC *dd2 = dynamic_cast<const mglDataC *>(d2);
++ const mglData *rd1 = dynamic_cast<const mglData *>(d1);
++ const mglData *rd2 = dynamic_cast<const mglData *>(d2);
++ if(dd1)
++#pragma omp parallel for
++ for(long i=0;i<nn;i++)
++ { a[2*i] = real(dd1->a[i]); a[2*i+1] = imag(dd1->a[i]); }
++ else if(rd1)
++#pragma omp parallel for
++ for(long i=0;i<nn;i++) a[2*i] = rd1->a[i];
++ else
++#pragma omp parallel for
++ for(long i=0;i<nn;i++) a[2*i] = d1->vthr(i);
++ if(dd2)
++#pragma omp parallel for
++ for(long i=0;i<nn;i++)
++ { b[2*i] = real(dd2->a[i]); b[2*i+1] = imag(dd2->a[i]); }
++ else if(rd2)
++#pragma omp parallel for
++ for(long i=0;i<nn;i++) b[2*i] = rd2->a[i];
++ else
++#pragma omp parallel for
++ for(long i=0;i<nn;i++) b[2*i] = d2->vthr(i);
++
++ if(strchr(dir,'x') && nx>1)
++ {
++ if(mgl_fft_data.wnx==nx) wt = mgl_fft_data.wtx;
++ else { clear = true; wt = mgl_fft_alloc(nx,0,0); }
++ mglStartThreadT(mgl_corx,ny*nz,b,a,wt,0,par);
++ if(mgl_fft_data.wnx==0)
++ { mgl_fft_data.wtx = wt; clear = false; mgl_fft_data.wnx=nx; }
++ }
++ if(strchr(dir,'y') && ny>1)
++ {
++ if(mgl_fft_data.wny==ny) wt = mgl_fft_data.wty;
++ else { clear = true; wt = mgl_fft_alloc(ny,0,0); }
++ mglStartThreadT(mgl_cory,nx*nz,b,a,wt,0,par);
++ if(mgl_fft_data.wny==0)
++ { mgl_fft_data.wty = wt; clear = false; mgl_fft_data.wny=ny; }
++ }
++ if(strchr(dir,'z') && nz>1)
++ {
++ if(mgl_fft_data.wnz==nz) wt = mgl_fft_data.wtz;
++ else { clear = true; wt = mgl_fft_alloc(nz,0,0); }
++ mglStartThreadT(mgl_corz,nx*ny,b,a,wt,0,par);
++ if(mgl_fft_data.wnz==0)
++ { mgl_fft_data.wtz = wt; clear = false; mgl_fft_data.wnz=nz; }
++ }
++ if(clear) mgl_fft_free(wt,0,0);
++ return a;
++}
++//-----------------------------------------------------------------------------
++HADT MGL_EXPORT mgl_datac_correl(HCDT d1, HCDT d2, const char *dir)
++{
++ double *a = mgl_d_correl(d1,d2,dir);
++ if(!a) return 0;
++ const long nx = d1->GetNx(), ny = d1->GetNy(), nz = d1->GetNz();
++ mglDataC *res = new mglDataC(nx,ny,nz);
++#pragma omp parallel for
++ for(long i=0;i<nx*ny*nz;i++) res->a[i] = dual(a[2*i], a[2*i+1]);
++ delete []a; return res;
++}
++uintptr_t MGL_EXPORT mgl_datac_correl_(uintptr_t *d1, uintptr_t *d2, const char *dir,int l)
++{ char *s=new char[l+1]; memcpy(s,dir,l); s[l]=0;
++ uintptr_t res = uintptr_t(mgl_datac_correl(_DA_(d1),_DA_(d2),s));
++ delete []s; return res; }
++//-----------------------------------------------------------------------------
++HMDT MGL_EXPORT mgl_data_correl(HCDT d1, HCDT d2, const char *dir)
++{
++ double *a = mgl_d_correl(d1,d2,dir); // NOTE: this is not so effective but straightforward way
++ if(!a) return 0;
++ const long nx = d1->GetNx(), ny = d1->GetNy(), nz = d1->GetNz();
++ mglData *res = new mglData(nx,ny,nz);
++#pragma omp parallel for
++ for(long i=0;i<nx*ny*nz;i++) res->a[i] = a[2*i];
++ delete []a; return res;
++}
++uintptr_t MGL_EXPORT mgl_data_correl_(uintptr_t *d1, uintptr_t *d2, const char *dir,int l)
++{ char *s=new char[l+1]; memcpy(s,dir,l); s[l]=0;
++ uintptr_t res = uintptr_t(mgl_datac_correl(_DA_(d1),_DA_(d2),s));
++ delete []s; return res; }
++//-----------------------------------------------------------------------------
++//-----------------------------------------------------------------------------
++//-----------------------------------------------------------------------------
projects / mathgl.git / commitdiff