[mathgl] 02/04: Merge tag 'upstream/2.2'
Dimitrios Eftaxiopoulos
eftaxiop-guest at moszumanska.debian.org
Mon Dec 30 15:55:18 UTC 2013
This is an automated email from the git hooks/post-receive script.
eftaxiop-guest pushed a commit to branch master
in repository mathgl.
commit 05cd6af14b504b2d4bf1c0cec47b1ed6f8496ae8
Merge: 1dc510d 9beecfb
Author: Dimitrios Eftaxiopoulos <eftaxi12 at otenet.gr>
Date: Sat Dec 28 16:23:27 2013 +0200
Merge tag 'upstream/2.2'
Upstream version 2.2
Conflicts:
CMakeLists.txt
src/fft.cpp
CMakeLists.txt | 206 ++-
CMakeLists.txt => CMakeLists.txt.orig | 186 ++-
ChangeLog.txt | 32 +
brush.ods | Bin 10909 -> 19306 bytes
examples/full_test.cpp | 118 +-
examples/glut_example.cpp | 1 +
examples/mpi_test.cpp | 39 +
examples/qt_example.cpp | 10 +-
examples/samples.cpp | 497 +++++--
examples/wnd_samples.cpp | 50 +-
examples/wx_example.cpp | 105 +-
include/CMakeLists.txt | 1 +
include/config.h.in | 31 +-
include/mgl2/abstract.h | 105 ++
include/mgl2/addon.h | 4 +-
include/mgl2/base.h | 207 ++-
include/mgl2/base_cf.h | 26 +-
include/mgl2/canvas.h | 161 +-
include/mgl2/canvas_cf.h | 9 +-
include/mgl2/canvas_wnd.h | 2 +-
include/mgl2/cont.h | 6 +-
include/mgl2/data.h | 56 +-
include/mgl2/data_cf.h | 56 +-
include/mgl2/datac.h | 27 +-
include/mgl2/datac_cf.h | 18 +-
include/mgl2/define.h | 109 +-
include/mgl2/eval.h | 2 +-
include/mgl2/evalc.h | 2 +-
include/mgl2/fit.h | 2 +-
include/mgl2/fltk.h | 18 +-
include/mgl2/font.h | 31 +-
include/mgl2/glut.h | 2 +-
include/mgl2/mgl.h | 82 +-
include/mgl2/mpi.h | 3 +-
include/mgl2/opengl.h | 23 +-
include/mgl2/other.h | 7 +-
include/mgl2/parser.h | 28 +-
include/mgl2/pde.h | 56 +
include/mgl2/plot.h | 17 +-
include/mgl2/prim.h | 2 +-
include/mgl2/qmathgl.h | 19 +-
include/mgl2/qt.h | 4 +-
include/mgl2/surf.h | 2 +-
include/mgl2/thread.h | 86 ++
include/mgl2/type.h | 8 +-
include/mgl2/vect.h | 2 +-
include/mgl2/volume.h | 2 +-
include/mgl2/window.h | 7 +-
include/mgl2/wnd.h | 5 +-
include/mgl2/wx.h | 40 +-
json/main.js | 38 +-
json/mathgl.Graph.js | 859 ++++++-----
json/mathgl.WebkitBackend.js | 18 +-
lang/CMakeLists.txt | 13 +-
include/mgl2/data.h => lang/data.i | 976 ++++++------
lang/mathgl.i | 18 +-
include/mgl2/mgl.h => lang/mgl.i | 2638 ++++++++++++++++-----------------
lang/sed_rules | 13 +
include/mgl2/type.h => lang/type.i | 68 +-
mathgl-2x.cbp | 19 +-
mgllab/editor.cpp | 3 +-
src/CMakeLists.txt | 19 +-
src/addon.cpp | 67 +-
src/axis.cpp | 320 ++--
src/base.cpp | 443 +++---
src/base_cf.cpp | 48 +-
src/canvas.cpp | 547 +++----
src/canvas_cf.cpp | 4 +
src/complex.cpp | 611 +++++---
src/complex_io.cpp | 217 ++-
src/cont.cpp | 309 ++--
src/crust.cpp | 235 +--
src/data.cpp | 1411 ++++++++++--------
src/{data_new.cpp => data_ex.cpp} | 471 +++---
src/data_gr.cpp | 177 +++
src/data_io.cpp | 319 ++--
src/data_png.cpp | 99 +-
src/eval.cpp | 64 +-
src/evalc.cpp | 25 +-
src/evalp.cpp | 577 +++----
src/exec.cpp | 80 +-
src/export.cpp | 74 +-
src/export_2d.cpp | 129 +-
src/export_3d.cpp | 302 ++--
src/fft.cpp | 1141 +++++++++-----
src/fft.cpp.orig | 1315 ++++++++++++++++
src/fit.cpp | 280 ++--
src/font.cpp | 324 ++--
src/interp.hpp | 384 +++++
src/mpi.cpp | 12 +-
src/obj.cpp | 4 -
src/opengl.cpp | 170 ++-
src/other.cpp | 163 +-
src/parser.cpp | 482 +++---
src/pde.cpp | 380 ++---
src/pixel.cpp | 1341 +++++++++++++----
src/plot.cpp | 762 ++++++----
src/prc/writePRC.h | 7 -
src/prim.cpp | 218 +--
src/surf.cpp | 305 ++--
src/vect.cpp | 525 +++----
src/volume.cpp | 79 +-
src/window.cpp | 1 +
texinfo/CMakeLists.txt | 194 +--
texinfo/appendix_en.texi | 2 +
texinfo/appendix_ru.texi | 2 +
texinfo/concept_en.texi | 46 +-
texinfo/concept_ru.texi | 89 +-
texinfo/core_en.texi | 288 +++-
texinfo/core_ru.texi | 320 ++--
texinfo/data_en.texi | 130 +-
texinfo/data_ru.texi | 235 +--
texinfo/doc_en.texi | 126 +-
texinfo/doc_ru.texi | 124 +-
texinfo/ex_mgl_en.texi | 405 +++--
texinfo/ex_mgl_ru.texi | 419 ++++--
texinfo/example_en.texi | 481 ++++--
texinfo/example_ru.texi | 529 +++++--
texinfo/formats_en.texi | 7 +-
texinfo/formats_ru.texi | 7 +-
texinfo/index.html | 23 +-
texinfo/json.html | 2 +-
texinfo/mathgl.js | 74 +-
texinfo/mathgl_en.texi | 75 +-
texinfo/mathgl_ru.texi | 78 +-
texinfo/mgl_en.texi | 29 +-
texinfo/mgl_ru.texi | 31 +-
texinfo/other_en.texi | 6 +-
texinfo/other_ru.texi | 6 +-
texinfo/overview_en.texi | 10 +
texinfo/overview_ru.texi | 11 +
texinfo/parse_en.texi | 26 +-
texinfo/parse_ru.texi | 242 +--
texinfo/surf_cont_fog_g.png | Bin 30517 -> 0 bytes
texinfo/symbols_en.texi | 64 +-
texinfo/symbols_ru.texi | 68 +-
texinfo/time.texi | 112 ++
texinfo/time_big.texi | 112 ++
texinfo/title.html | 10 -
texinfo/toc_en.html | 44 -
texinfo/toc_fr.html | 45 -
texinfo/toc_ru.html | 54 -
texinfo/udav_en.texi | 6 +
texinfo/udav_ru.texi | 6 +
texinfo/version.texi | 2 +
texinfo/web_en.texi | 80 +-
texinfo/web_fr.texi | 4 +-
texinfo/web_ru.texi | 80 +-
texinfo/widget_en.texi | 178 ++-
texinfo/widget_ru.texi | 184 ++-
todo.txt | 55 +-
udav/dat_pnl.cpp | 15 +-
udav/info_dlg.cpp | 2 +-
udav/plot_pnl.cpp | 16 +-
udav/udav_wnd.cpp | 5 +-
widgets/CMakeLists.txt | 15 +-
widgets/fltk.cpp | 20 +-
widgets/glut.cpp | 42 +-
widgets/qt.cpp | 39 +-
widgets/wx.cpp | 303 +---
160 files changed, 16959 insertions(+), 10235 deletions(-)
diff --cc CMakeLists.txt.orig
index d91bf8c,0000000..1e69ca8
mode 100644,000000..100644
--- a/CMakeLists.txt.orig
+++ b/CMakeLists.txt.orig
@@@ -1,426 -1,0 +1,562 @@@
+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:
++ 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 1.3)
++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)
- # message(STATUS "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)
- # message(STATUS "Found libHaru library at: ${HPDF_LIB}")
- # message(STATUS "Found libHaru headers: ${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-all-widgets OR enable-wx) AND (NOT enable-lgpl) )
+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-all-widgets OR enable-wx) AND (NOT enable-lgpl) )
+else(enable-wx)
+ set(MGL_HAVE_WX 0)
- #endif((enable-all-widgets OR enable-wx) AND (NOT enable-lgpl) )
+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)
- set(MGL_HAVE_DOC 1)
++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)
- find_program(findth texi2html)
- if(NOT findth)
- message(SEND_ERROR "Couldn't find texi2html needed for documentation building.")
- endif(NOT findth)
++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)
- set(MGL_HAVE_DOC 0)
- endif(enable-doc)
++else(enable-doc-pdf-ru)
++ set(MGL_HAVE_DOC_PDF_RU 0)
++endif(enable-doc-pdf-ru)
+
- if(MGL_HAVE_PYTHON OR MGL_HAVE_OCTAVE)
++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_PYTHON OR MGL_HAVE_OCTAVE)
-
- #execute_process(
- # COMMAND ${oct_prog} -q --eval \'strcat\(octave_config_info\("canonical_host_type"\), "-", octave_config_info\("api_version"\)\)'|sed -e 's/ans = //'
- # OUTPUT_VARIABLE OCTAVE_ARCH
- # )
- #message(SEND_ERROR "${OCTAVE_ARCH}")
++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)
diff --cc src/fft.cpp.orig
index 0000000,0000000..57e9326
new file mode 100644
--- /dev/null
+++ b/src/fft.cpp.orig
@@@ -1,0 -1,0 +1,1315 @@@
++/***************************************************************************
++ * fft.cpp is part of Math Graphic Library
++ * Copyright (C) 2007-2012 Alexey Balakin <mathgl.abalakin at 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; }
++//-----------------------------------------------------------------------------
++//-----------------------------------------------------------------------------
++//-----------------------------------------------------------------------------
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/mathgl.git
More information about the debian-science-commits
mailing list