diff --git a/CMakeLists.txt b/CMakeLists.txt index f9e0b83de..633ff4eb1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -84,7 +84,9 @@ option(Tasmanian_ENABLE_MPI "Enable MPI support for Tasmanian (stable)" option(Tasmanian_ENABLE_SWIG "Regenerate Fortran 2003 bindings with SWIG (developers only)" OFF) option(Tasmanian_ENABLE_DOXYGEN "Enable Doxygen documentation for Tasmanian (complete)" OFF) -option(Tasmanian_ENABLE_FORTRAN90 "(deprecated, untested) Enable the old Fortran 90 interface of Tasmanian" OFF) +if (Tasmanian_ENABLE_FORTRAN90) + message(FATAL_ERROR "Tasmanian_ENABLE_FORTRAN90 and the Fortran 90 interface have been removed since v8.2") +endif() # treat those similar to options, but give values other than ON/OFF set(Tasmanian_MATLAB_WORK_FOLDER "" CACHE PATH "Enable MATLAB interface and use the path for the MATLAB work folder (stable)") @@ -147,9 +149,6 @@ add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/Addons" "${CMAKE_CURRENT_BINA add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/Tasgrid/" "${CMAKE_CURRENT_BINARY_DIR}/Tasgrid") if (Tasmanian_ENABLE_FORTRAN) - if (Tasmanian_ENABLE_FORTRAN90) - add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/InterfaceFortran" "${CMAKE_CURRENT_BINARY_DIR}/Fortran") - endif() add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/InterfaceSwig" "${CMAKE_CURRENT_BINARY_DIR}/Fortran03") endif() diff --git a/Doxygen/InterfaceFortran2003.md b/Doxygen/InterfaceFortran2003.md index 68fb6486a..eb11c8785 100644 --- a/Doxygen/InterfaceFortran2003.md +++ b/Doxygen/InterfaceFortran2003.md @@ -2,7 +2,7 @@ Tasmanian comes with Fortran 2003 interface generated by [Swig-Fortran](https://github.com/swig-fortran). The automatic generation makes the code easy to maintain and support, and the generation process is done off-line by the development team; hence, users are not required to install Swig and can use the interface with just a regular Fortran compiler. Furthermore, the 2003 standard support objects which allows for the OOP interface of Tasmanian C++ to be duplicated similar to Python and thus the on-line C++ Doxygen documentation is now relevant to Fortran as well. -The new module is called `Tasmanian` (without the sg suffix), and it is incompatible with the Fortran 90/95, i.e., use one or the other but not both. A `TasmanianSparseGrid` type is introduced as before, but the interface is fully object-oriented. Currently, the Tasmanian Sparse Grid module has been implemented together with the MPI Addons; the Dream module and the remainder of the Addon templates are not supported (possibly in the future). +The module is called `Tasmanian` and requires Fortran 2003 and a class `TasmanianSparseGrid` type is used. Currently, the Tasmanian Sparse Grid module has been implemented together with the MPI Addons; the Dream module and the remainder of the Addon templates are not supported (possibly in the future). Simple example: ``` diff --git a/InterfaceFortran/CMakeLists.txt b/InterfaceFortran/CMakeLists.txt deleted file mode 100644 index 1f54ed865..000000000 --- a/InterfaceFortran/CMakeLists.txt +++ /dev/null @@ -1,62 +0,0 @@ -######################################################################## -# Fortran librareis and command line tools -######################################################################## - -add_library(Tasmanian_libfortran90 TasmanianSG.f90 - tsgC2FortranBridge.f90 - tsgC2Fortran.cpp) - -target_include_directories(Tasmanian_libfortran90 PUBLIC $) -target_link_libraries(Tasmanian_libfortran90 Tasmanian_addons) - -set_target_properties(Tasmanian_libfortran90 PROPERTIES OUTPUT_NAME "tasmanianfortran90" - SOVERSION ${Tasmanian_VERSION_MAJOR} - VERSION ${PROJECT_VERSION}) -Tasmanian_rpath_target(TARGET Tasmanian_libfortran90 COMPONENTS SparseGrids DREAM) - -install(TARGETS Tasmanian_libfortran90 - EXPORT "${Tasmanian_export_name}" - RUNTIME DESTINATION "bin" - LIBRARY DESTINATION "lib" - ARCHIVE DESTINATION "lib") - - -######################################################################## -# add the fortran tester and examples executables -######################################################################## -add_executable(Tasmanian_fortester fortester.f90) - -set_target_properties(Tasmanian_fortester PROPERTIES OUTPUT_NAME "fortester" LINKER_LANGUAGE Fortran) -Tasmanian_rpath_target(TARGET Tasmanian_fortester USE_CURRENT COMPONENTS SparseGrids DREAM) - -target_link_libraries(Tasmanian_fortester Tasmanian_libfortran90) - - -######################################################################## -# handle the MPI dependence and MPI tests -######################################################################## -if (Tasmanian_ENABLE_MPI) - target_link_libraries(Tasmanian_libfortran90 MPI::MPI_Fortran) - - add_executable(Tasmanian_mpif90 mpitester.f90) - set_target_properties(Tasmanian_mpif90 PROPERTIES OUTPUT_NAME "mpitester" LINKER_LANGUAGE Fortran) - target_link_libraries(Tasmanian_mpif90 Tasmanian_libfortran90) - add_test(MPIFortranGridIO ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 3 ${MPIEXEC_PREFLAGS} ${CMAKE_CURRENT_BINARY_DIR}/mpitester ${MPIEXEC_POSTFLAGS}) - set_tests_properties(MPIFortranGridIO PROPERTIES RUN_SERIAL ON) - Tasmanian_set_test_properties(TESTS MPIFortranGridIO) -endif() - - -######################################################################## -# Testing -######################################################################## -add_test(Fortran90 fortester) -Tasmanian_set_test_properties(TESTS Fortran90) - - -######################################################################## -# Installation -######################################################################## -install(FILES "${CMAKE_CURRENT_BINARY_DIR}/tasmaniansg.mod" - DESTINATION include - PERMISSIONS OWNER_EXECUTE OWNER_WRITE OWNER_READ GROUP_EXECUTE GROUP_READ WORLD_EXECUTE WORLD_READ) diff --git a/InterfaceFortran/Examples/example_sparse_grids.f90 b/InterfaceFortran/Examples/example_sparse_grids.f90 deleted file mode 100644 index 549f892b5..000000000 --- a/InterfaceFortran/Examples/example_sparse_grids.f90 +++ /dev/null @@ -1,1038 +0,0 @@ -!================================================================================================================================================================================== -! Copyright (c) 2017, Miroslav Stoyanov -! -! This file is part of -! Toolkit for Adaptive Stochastic Modeling And Non-Intrusive ApproximatioN: TASMANIAN -! -! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: -! -! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. -! -! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions -! and the following disclaimer in the documentation and/or other materials provided with the distribution. -! -! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse -! or promote products derived from this software without specific prior written permission. -! -! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, -! INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. -! IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, -! OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, -! OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -! -! UT-BATTELLE, LLC AND THE UNITED STATES GOVERNMENT MAKE NO REPRESENTATIONS AND DISCLAIM ALL WARRANTIES, BOTH EXPRESSED AND IMPLIED. -! THERE ARE NO EXPRESS OR IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, OR THAT THE USE OF THE SOFTWARE WILL NOT INFRINGE ANY PATENT, -! COPYRIGHT, TRADEMARK, OR OTHER PROPRIETARY RIGHTS, OR THAT THE SOFTWARE WILL ACCOMPLISH THE INTENDED RESULTS OR THAT THE SOFTWARE OR ITS USE WILL NOT RESULT IN INJURY OR DAMAGE. -! THE USER ASSUMES RESPONSIBILITY FOR ALL LIABILITIES, PENALTIES, FINES, CLAIMS, CAUSES OF ACTION, AND COSTS AND EXPENSES, CAUSED BY, RESULTING FROM OR ARISING OUT OF, -! IN WHOLE OR IN PART THE USE, STORAGE OR DISPOSAL OF THE SOFTWARE. -!================================================================================================================================================================================== -PROGRAM TasmanianSGExample - USE TasmanianSG -IMPLICIT NONE - type(TasmanianSparseGrid) :: grid, grid1, grid2, grid3 - INTEGER :: dims, outs, level - INTEGER :: N1, N2, N3 - INTEGER :: N, i, j, verm, vern - DOUBLE PRECISION :: err1, err2, err3, err4, exact - REAL :: cpuStart, cpuEnd, stages(2,3) - INTEGER :: conformal(3) - CHARACTER, pointer :: string(:) - DOUBLE PRECISION, pointer :: points(:,:), weights(:) - DOUBLE PRECISION :: x, y, integ, E - DOUBLE PRECISION, allocatable :: transformA(:), transformB(:), values(:,:), tvalues(:,:) - DOUBLE PRECISION, allocatable :: res(:), res2(:,:) - DOUBLE PRECISION :: desired_x(2) - DOUBLE PRECISION :: randPoints(4,1000), randPoints2(2,1000), randPoints3(3,1000) - DOUBLE PRECISION :: PI = 4.D0*DATAN(1.D0) - - -! This is the sound "Glaucodon Ballaratensis" makes :) -! WRITE(*,*) "Ghurrrrrphurrr" - -! ============ reference table of rules ============ ! -! 1: clenshaw-curtis 2: clenshaw-curtis-zero -! 3: chebyshev 4: chebyshev-odd -! 5: gauss-legendre 6: gauss-legendreodd -! 7: gauss-patterson 8: leja -! 9: lejaodd 10: rleja -! 11: rleja-odd 12: rleja-double-2 -! 13: rleja-double-4 14: rleja-shifted -! 15: rleja-shifted-even 16: rleja-shifted-double -! 17: max-lebesgue 18: max-lebesgue-odd -! 19: min-lebesgue 20: min-lebesgue-odd -! 21: min-delta 22: min-delta-odd -! 23: gauss-chebyshev-1 24: gauss-chebyshev-1-odd -! 25: gauss-chebyshev-2 26: gauss-chebyshev-2-odd -! 27: fejer2 28: gauss-gegenbauer -! 29: gauss-gegenbauer-odd 30: gauss-jacobi -! 31: gauss-jacobi-odd 32: gauss-laguerre -! 33: gauss-laguerre-odd 34: gauss-hermite -! 35: gauss-hermite-odd 36: custom-tabulated -! 37: localp 38: localp-zero -! 39: semi-localp 40: wavelet - -! ============ reference table of grid types ============ ! -! 1: level 2: curved -! 3: iptotal 4: ipcurved -! 5: qptotal 6: qpcurved -! 7: hyperbolic 8: iphyperbolic -! 9: qphyperbolic 10: tensor -! 11: iptensor 12: qptensor - -! ============ reference table of refinement types ============ ! -! 1: classic 2: parents first -! 3: directional 4: FDS (both parents and directions) - -! ============ reference table of acceleration types ============ ! -! 0: none 1: CPU BLAS 2: GPU cuBLAS -! 3: GPU CUDA 4: GPU MAGMA - -! ==================================================================== ! -! EXAMPLE 1: integrate: f(x,y) = exp(-x^2) * cos(y) over [-1,1] x [-1,1] -! using classical Smolyak grid with Clenshaw-Curtis points and weights - - - verm = tsgGetVersionMajor() - vern = tsgGetVersionMinor() - - ! WARNING: do not DEALLOCATE the string pointer, it is const char* - string => tsgGetLicense() - WRITE(*,*) "-------------------------------------------------------------------------------------------------" - WRITE(*,*) "Tasmanian Sparse Grids Fortran Module (TasmanianSG)" - WRITE(*,*) "-------------------------------------------------------------------------------------------------" - WRITE(*,"(A,I2,A,I1)") " Tasmanian Sparse Grid module, version: ", verm, ".", vern - WRITE(*,"(A,40A)") " license: ", string - WRITE(*,*) - - WRITE(*,*) "-------------------------------------------------------------------------------------------------" - WRITE(*,*) "Example 1: integrate f(x,y) = exp(-x^2) * cos(y), using clenshaw-curtis level nodes" - - dims = 2 - level = 6 - -! before you use a grid, you must allocate a new grid - call tsgAllocateGrid(grid) - - CALL tsgMakeGlobalGrid(grid, dims, 0, level, tsg_level, tsg_clenshaw_curtis) - points => tsgGetPoints(grid) - weights => tsgGetQuadratureWeights(grid) - - N = tsgGetNumPoints(grid) - integ = 0.0 - - DO i = 1, N - x = points(1, i) - y = points(2, i) - integ = integ + weights(i) * exp(-x*x) * cos(y) - END DO - - E = abs(integ - 2.513723354063905D+00) - - WRITE(*,"(A,I4)") " at level: ", level - WRITE(*,"(A,I4,A)") " the grid has: ", N, " points" - WRITE(*,"(A,E25.16)") " integral: ", integ - WRITE(*,"(A,E25.16)") " error: ", E - WRITE(*,*) - - level = 7 -! no need to ask for a new ID when remaking an existing grid - CALL tsgMakeGlobalGrid(grid, dims, 0, level, tsg_level, tsg_clenshaw_curtis) - -! do not forget to release the memory associated with points and weights - DEALLOCATE(points) - DEALLOCATE(weights) - points => tsgGetPoints(grid) - weights => tsgGetQuadratureWeights(grid) - - N = tsgGetNumPoints(grid) - integ = 0.0 - - DO i = 1, N - x = points(1, i) - y = points(2, i) - integ = integ + weights(i) * exp(-x*x) * cos(y) - END DO - - E = abs(integ - 2.513723354063905D+00) - WRITE(*,"(A,I4)") " at level: ", level - WRITE(*,"(A,I4,A)") " the grid has: ", N, " points" - WRITE(*,"(A,E25.16)") " integral: ", integ - WRITE(*,"(A,E25.16)") " error: ", E - WRITE(*,*) - - DEALLOCATE(points) - DEALLOCATE(weights) - -! after calling tsgDeallocateGrid(), we can no longer use this grid -! until we call tsgAllocateGrid() again - CALL tsgDeallocateGrid(grid) - -! ==================================================================== ! -! EXAMPLE 2: integrate: f(x,y) = exp(-x^2) * cos(y) -! over (x,y) in [-5,5] x [-2,3] -! using Gauss-Patterson rules chosen to integrate exactly polynomials of -! total degree up to degree specified by prec - - WRITE(*,*) "-------------------------------------------------------------------------------------------------" - WRITE(*,*) "Example 2: integrate f(x,y) = exp(-x^2) * cos(y) over [-5,5] x [-2,3] using Gauss-Patterson nodes" - - dims = 2 - level = 20 - - ALLOCATE(transformA(dims)) - ALLOCATE(transformB(dims)) - transformA(1) = -5.0 - transformA(2) = -2.0 - transformB(1) = 5.0 - transformB(2) = 3.0 - -! need new grid, since we freed this earlier - call tsgAllocateGrid(grid) - -! gauss-patterson = 7, type_qptotal = 5 - CALL tsgMakeGlobalGrid(grid, dims, 0, level, tsg_qptotal, tsg_gauss_patterson) - CALL tsgSetDomainTransform(grid, transformA, transformB) - - points => tsgGetPoints(grid) - weights => tsgGetQuadratureWeights(grid) - - N = tsgGetNumPoints(grid) - integ = 0.0 - - DO i = 1, N - x = points(1, i) - y = points(2, i) - integ = integ + weights(i) * exp(-x*x) * cos(y) - END DO - - E = abs(integ - 1.861816427518323D+00) - WRITE(*,"(A,I4)") " at precision: ", level - WRITE(*,"(A,I4,A)") " the grid has: ", N, " points" - WRITE(*,"(A,E25.16)") " integral: ", integ - WRITE(*,"(A,E25.16)") " error: ", E - WRITE(*,*) - - level = 40 -! no need to ask for a new ID when remaking an existing grid - CALL tsgMakeGlobalGrid(grid, dims, 0, level, tsg_qptotal, tsg_gauss_patterson) - CALL tsgSetDomainTransform(grid, transformA, transformB) - -! do not forget to release the memory associated with points and weights - DEALLOCATE(points) - DEALLOCATE(weights) - points => tsgGetPoints(grid) - weights => tsgGetQuadratureWeights(grid) - - N = tsgGetNumPoints(grid) - integ = 0.0 - - DO i = 1, N - x = points(1, i) - y = points(2, i) - integ = integ + weights(i) * exp(-x*x) * cos(y) - END DO - - E = abs(integ - 1.861816427518323D+00) - WRITE(*,"(A,I4)") " at precision: ", level - WRITE(*,"(A,I4,A)") " the grid has: ", N, " points" - WRITE(*,"(A,E25.16)") " integral: ", integ - WRITE(*,"(A,E25.16)") " error: ", E - WRITE(*,*) - - DEALLOCATE(points) - DEALLOCATE(weights) - ! keep transformA and transformB for the next example - CALL tsgDeallocateGrid(grid) - -! ==================================================================== ! -! EXAMPLE 3: integrate: f(x,y) = exp(-x^2) * cos(y) -! over (x,y) in [-5,5] x [-2,3] -! using different rules - - call tsgAllocateGrid(grid1) - call tsgAllocateGrid(grid2) - call tsgAllocateGrid(grid3) - - WRITE(*,*) "-------------------------------------------------------------------------------------------------" - WRITE(*,*) "Example 3: integrate f(x,y) = exp(-x^2) * cos(y) over [-5,5] x [-2,3] using different rules" - WRITE(*,*) - - WRITE(*,*) " Clenshaw-Curtis Gauss-Legendre Gauss-Patterson" - WRITE(*,*) " precision points error points error points error" - - DO level = 9, 30, 4 - CALL tsgMakeGlobalGrid(grid1, dims, 0, level, tsg_qptotal, tsg_clenshaw_curtis) - CALL tsgSetDomainTransform(grid1, transformA, transformB) - CALL tsgMakeGlobalGrid(grid2, dims, 0, level, tsg_qptotal, tsg_gauss_legendre) - CALL tsgSetDomainTransform(grid2, transformA, transformB) - CALL tsgMakeGlobalGrid(grid3, dims, 0, level, tsg_qptotal, tsg_gauss_patterson) - CALL tsgSetDomainTransform(grid3, transformA, transformB) - - points => tsgGetPoints(grid1) - weights => tsgGetQuadratureWeights(grid1) - N1 = tsgGetNumPoints(grid1) - integ = 0.0 - DO i = 1, N1 - x = points(1, i) - y = points(2, i) - integ = integ + weights(i) * exp(-x*x) * cos(y) - END DO - err1 = abs(integ - 1.861816427518323D+00) - DEALLOCATE(points) - DEALLOCATE(weights) - - points => tsgGetPoints(grid2) - weights => tsgGetQuadratureWeights(grid2) - N2 = tsgGetNumPoints(grid2) - integ = 0.0 - DO i = 1, N2 - x = points(1, i) - y = points(2, i) - integ = integ + weights(i) * exp(-x*x) * cos(y) - END DO - err2 = abs(integ - 1.861816427518323D+00) - DEALLOCATE(points) - DEALLOCATE(weights) - - points => tsgGetPoints(grid3) - weights => tsgGetQuadratureWeights(grid3) - N3 = tsgGetNumPoints(grid3) - integ = 0.0 - DO i = 1, N3 - x = points(1, i) - y = points(2, i) - integ = integ + weights(i) * exp(-x*x) * cos(y) - END DO - err3 = abs(integ - 1.861816427518323D+00) - DEALLOCATE(points) - DEALLOCATE(weights) - - WRITE(*,"(I10,I10,E11.3,I9,E11.3,I9,E11.3)") level, N1, err1, N2, err2, N3, err3 - - END DO - WRITE(*,*) - - CALL tsgDeallocateGrid(grid1) - CALL tsgDeallocateGrid(grid2) - CALL tsgDeallocateGrid(grid3) - - DEALLOCATE(transformA) - DEALLOCATE(transformB) - - -! ==================================================================== ! -! EXAMPLE 4: interpolate: f(x,y) = exp(-x^2) * cos(y) -! with a rule that exactly interpolates polynomials of total degree - - call tsgAllocateGrid(grid) - - WRITE(*,*) "-------------------------------------------------------------------------------------------------" - WRITE(*,*) "Example 4: interpolate f(x,y) = exp(-x^2) * cos(y), using clenshaw-curtis iptotal rule" - WRITE(*,*) - - dims = 2 - outs = 1 - level = 10 - - desired_x(1) = 0.3 - desired_x(2) = 0.7 - - ! desired value - exact = exp(-desired_x(1)**2) * cos(desired_x(2)) - - CALL tsgMakeGlobalGrid(grid, dims, outs, level, tsg_iptotal, tsg_clenshaw_curtis) - - N = tsgGetNumNeeded(grid) - points => tsgGetNeededPoints(grid) - ALLOCATE(values(outs,N)) - - DO i = 1, N - x = points(1, i) - y = points(2, i) - values(1,i) = exp(-x**2) * cos(y) - END DO - - CALL tsgLoadNeededPoints(grid, values) - DEALLOCATE(values) - DEALLOCATE(points) - - ALLOCATE(res(outs)) ! will DEALLOCATE later - CALL tsgEvaluate(grid, desired_x, res) - E = abs(res(1) - exact) - - WRITE(*,"(A,I4)") " using polynomials of total degree: ", level - WRITE(*,"(A,I4,A)") " the grid has: ", N, " points" - WRITE(*,"(A,E25.16)") " interpolant at (0.3,0.7): ", res(1) - WRITE(*,"(A,E25.16)") " error: ", E - WRITE(*,*) - - ! do the same with level = 12 - level = 12 - - CALL tsgMakeGlobalGrid(grid, dims, outs, level, tsg_iptotal, tsg_clenshaw_curtis) - - N = tsgGetNumNeeded(grid) - points => tsgGetNeededPoints(grid) - ALLOCATE(values(outs,N)) - - DO i = 1, N - x = points(1, i) - y = points(2, i) - values(1,i) = exp(-x**2) * cos(y) - END DO - - CALL tsgLoadNeededPoints(grid, values) - DEALLOCATE(values) - DEALLOCATE(points) - - CALL tsgEvaluate(grid, desired_x, res) - E = abs(res(1) - exact) - - WRITE(*,"(A,I4)") " using polynomials of total degree: ", level - WRITE(*,"(A,I4,A)") " the grid has: ", N, " points" - WRITE(*,"(A,E25.16)") " interpolant at (0.3,0.7): ", res(1) - WRITE(*,"(A,E25.16)") " error: ", E - WRITE(*,*) - DEALLOCATE(res) - - CALL tsgDeallocateGrid(grid) - -! prepare random smaples for future tests - call random_seed() - - call random_number(randPoints) - -! ==================================================================== ! -! EXAMPLE 5: -! interpolate: f(x1,x2,x3,x4) = exp(-x1^2) * cos(x2) * exp(-x3^2) * cos(x4) -! with Global and Sequence Leja rules -! Skipping example 5: it takes time and slows down testing, but it dosn't show/test any API not covered in other cases - -! dims = 4 -! outs = 1 -! level = 15 -! -! CALL cpu_time(cpuStart) -! CALL tsgMakeGlobalGrid(grid, dims, outs, level, tsg_level, tsg_leja) -! CALL cpu_time(cpuEnd) -! stages(1,1) = cpuEnd - cpuStart -! -! N = tsgGetNumPoints(grid) -! -! WRITE(*,*) "-------------------------------------------------------------------------------------------------" -! WRITE(*,*) "Example 5: interpolate f(x1,x2,x3,x4) = exp(-x1^2) * cos(x2) * exp(-x3^2) * cos(x4)" -! WRITE(*,*) " comparign the performance of Global and Sequence grids with leja nodes" -! WRITE(*,"(A,I4)") " using polynomials of total degree up to: ", level -! WRITE(*,"(A,I4,A)") " the grids have: ", N, " points" -! WRITE(*,*) " both grids are evaluated at 1000 random points " -! WRITE(*,*) -! -! points => tsgGetNeededPoints(grid) -! ALLOCATE(values(outs,N)) -! DO i = 1, N -! values(1,i) = exp(-points(1,i)**2) * cos(points(2,i)) * exp(-points(3,i)**2) * cos(points(4,i)) -! END DO -! DEALLOCATE(points) -! -! CALL cpu_time(cpuStart) -! CALL tsgLoadNeededPoints(grid, values) -! CALL cpu_time(cpuEnd) -! stages(1,2) = cpuEnd - cpuStart -! -! ALLOCATE(res2(outs,1000)) ! 2-D result -! -! CALL cpu_time(cpuStart) -! CALL tsgEvaluateBatch(grid, randPoints, 1000, res2) -! CALL cpu_time(cpuEnd) -! stages(1,3) = cpuEnd - cpuStart -! -! CALL cpu_time(cpuStart) -! CALL tsgMakeSequenceGrid(grid, dims, outs, level, tsg_level, tsg_leja) -! CALL cpu_time(cpuEnd) -! stages(2,1) = cpuEnd - cpuStart -! -! ! points are the same, no need to recompue values -! CALL cpu_time(cpuStart) -! CALL tsgLoadNeededPoints(grid, values) -! CALL cpu_time(cpuEnd) -! stages(2,2) = cpuEnd - cpuStart -! -! DEALLOCATE(values) -! -! CALL cpu_time(cpuStart) -! CALL tsgEvaluateBatch(grid, randPoints, 1000, res2) -! CALL cpu_time(cpuEnd) -! stages(2,3) = cpuEnd - cpuStart -! -! WRITE(*,*) "Stage Global Grid Sequence Grid" -! WRITE(*,"(A,E20.8,E20.8)") " make grid ", stages(1,1), stages(2,1) -! WRITE(*,"(A,E20.8,E20.8)") " load needed", stages(1,2), stages(2,2) -! WRITE(*,"(A,E20.8,E20.8)") " evaluate ", stages(1,3), stages(2,3) -! -! CALL tsgDeallocateGrid(grid) -! DEALLOCATE(res2) - -! ==================================================================== ! -! EXAMPLE 6: -! interpolate: f(x,y) = exp(-x^2) * cos(y) -! using different refinement schemes - - ALLOCATE(tvalues(1,1000)) ! true values of f(x,y) - DO i = 1, 1000 - randPoints2(1,i) = randPoints(1,i) - randPoints2(2,i) = randPoints(2,i) - tvalues(1,i) = exp(-randPoints(1,i)**2) * cos(randPoints(2,i)) - END DO - - call tsgAllocateGrid(grid1) - call tsgAllocateGrid(grid2) - call tsgAllocateGrid(grid3) - - dims = 2 - outs = 1 - - CALL tsgMakeGlobalGrid(grid1, dims, outs, 3, tsg_iptotal, tsg_leja) - CALL tsgMakeGlobalGrid(grid2, dims, outs, 3, tsg_iptotal, tsg_leja) - CALL tsgMakeGlobalGrid(grid3, dims, outs, 3, tsg_iptotal, tsg_leja) - - N = tsgGetNumNeeded(grid1) - points => tsgGetNeededPoints(grid1) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = exp(-points(1,i)**2) * cos(points(2,i)) - END DO - CALL tsgLoadNeededPoints(grid1, values) - CALL tsgLoadNeededPoints(grid2, values) - CALL tsgLoadNeededPoints(grid3, values) - DEALLOCATE(values) - DEALLOCATE(points) - - WRITE(*,*) "-------------------------------------------------------------------------------------------------" - WRITE(*,*) "Example 6: interpolate: f(x,y) = exp(-x^2) * cos(y)" - WRITE(*,*) " using leja nodes and different refinement schemes " - WRITE(*,*) " the error is estimated as the maximum from 1000 random points" - - WRITE(*,*) " Total Degree Curved Surplus" - WRITE(*,*) "iteration points error points error points error" - - ALLOCATE(res2(outs,1000)) ! 2-D result - - ! iptotal: 3, ipcurved: 4 - DO j = 1, 10 - CALL tsgSetAnisotropicRefinement(grid1, tsg_iptotal, 10, 1) - - N = tsgGetNumNeeded(grid1) - points => tsgGetNeededPoints(grid1) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = exp(-points(1,i)**2) * cos(points(2,i)) - END DO - CALL tsgLoadNeededPoints(grid1, values) - DEALLOCATE(values) - DEALLOCATE(points) - - CALL tsgEvaluateBatch(grid1, randPoints2, 1000, res2) - err1 = 0.0 - DO i = 1, 1000 - IF(abs(res2(1,i) - tvalues(1,i)) .GT. err1)then - err1 = abs(res2(1,i) - tvalues(1,i)) - ENDIF - END DO - - CALL tsgSetAnisotropicRefinement(grid2, tsg_ipcurved, 10, 1) - - N = tsgGetNumNeeded(grid2) - points => tsgGetNeededPoints(grid2) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = exp(-points(1,i)**2) * cos(points(2,i)) - END DO - CALL tsgLoadNeededPoints(grid2, values) - DEALLOCATE(values) - DEALLOCATE(points) - - CALL tsgEvaluateBatch(grid2, randPoints2, 1000, res2) - err2 = 0.0 - DO i = 1, 1000 - IF(abs(res2(1,i) - tvalues(1,i)) .GT. err2)then - err2 = abs(res2(1,i) - tvalues(1,i)) - ENDIF - END DO - - CALL tsgSetGlobalSurplusRefinement(grid3, 1.D-10, 1) - - N = tsgGetNumNeeded(grid3) - points => tsgGetNeededPoints(grid3) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = exp(-points(1,i)**2) * cos(points(2,i)) - END DO - CALL tsgLoadNeededPoints(grid3, values) - DEALLOCATE(values) - DEALLOCATE(points) - - CALL tsgEvaluateBatch(grid3, randPoints2, 1000, res2) - err3 = 0.0 - DO i = 1, 1000 - IF(abs(res2(1,i) - tvalues(1,i)) .GT. err3)then - err3 = abs(res2(1,i) - tvalues(1,i)) - ENDIF - END DO - - N1 = tsgGetNumPoints(grid1) - N2 = tsgGetNumPoints(grid2) - N3 = tsgGetNumPoints(grid3) - - WRITE(*,"(I9,I9,E12.4,I9,E12.4,I9,E12.4)") j, N1, err1, N2, err2, N3, err3 - - END DO - WRITE(*,*) - CALL tsgDeallocateGrid(grid1) - CALL tsgDeallocateGrid(grid2) - CALL tsgDeallocateGrid(grid3) - -! ==================================================================== ! -! EXAMPLE 7: -! interpolate: f(x,y) = exp(-x^2) * cos(y) -! using localp and semilocalp grids - - WRITE(*,*) "-------------------------------------------------------------------------------------------------" - WRITE(*,*) "Example 7: interpolate: f(x,y) = exp(-x^2) * cos(y)" - WRITE(*,*) " using localp and semi-localp rules with depth 7" - WRITE(*,*) " the error is estimated as the maximum from 1000 random points" - WRITE(*,*) - - call tsgAllocateGrid(grid1) - call tsgAllocateGrid(grid2) - - dims = 2 - outs = 1 - - CALL tsgMakeLocalPolynomialGrid(grid1, dims, outs, 7, 2, tsg_localp) - CALL tsgMakeLocalPolynomialGrid(grid2, dims, outs, 7, 2, tsg_semi_localp) - - N = tsgGetNumNeeded(grid1) - points => tsgGetNeededPoints(grid1) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = exp(-points(1,i)**2) * cos(points(2,i)) - END DO - CALL tsgLoadNeededPoints(grid1, values) - DEALLOCATE(values) - DEALLOCATE(points) - - CALL tsgEvaluateBatch(grid1, randPoints2, 1000, res2) - err1 = 0.0 - DO i = 1, 1000 - IF(abs(res2(1,i) - tvalues(1,i)) .GT. err1)then - err1 = abs(res2(1,i) - tvalues(1,i)) - ENDIF - END DO - - points => tsgGetNeededPoints(grid2) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = exp(-points(1,i)**2) * cos(points(2,i)) - END DO - CALL tsgLoadNeededPoints(grid2, values) - DEALLOCATE(values) - DEALLOCATE(points) - - CALL tsgEvaluateBatch(grid2, randPoints2, 1000, res2) - err2 = 0.0 - DO i = 1, 1000 - IF(abs(res2(1,i) - tvalues(1,i)) .GT. err2)then - err2 = abs(res2(1,i) - tvalues(1,i)) - ENDIF - END DO - - WRITE(*,"(A,I5)") " Number of points: ", N - WRITE(*,"(A,E12.4)") " Error for rule_localp: ", err1 - WRITE(*,"(A,E12.4)") " Error for rule_semilocalp: ", err2 - WRITE(*,*) " Note: semi-localp wins this competition because the function is very smooth" - WRITE(*,*) - -! ==================================================================== ! -! EXAMPLE 8: -! interpolate: f(x,y) = cos(0.5 * pi * x) * cos(0.5 * pi * y) -! using localp and semilocalp grids - -! remake the true values for the next example - DO i = 1, 1000 - tvalues(1,i) = cos(0.5 * PI * randPoints2(1,i)) & - * cos(0.5 * PI * randPoints2(2,i)) - END DO - - WRITE(*,*) "-------------------------------------------------------------------------------------------------" - WRITE(*,*) "Example 8: interpolate f(x,y) = cos(0.5 * pi * x) * cos(0.5 * pi * y)" - WRITE(*,*) " using localp and localp-zero rules with depths 7 and 6" - WRITE(*,*) " the error is estimated as the maximum from 1000 random points" - WRITE(*,*) - - dims = 2 - outs = 1 - - CALL tsgMakeLocalPolynomialGrid(grid1, dims, outs, 7, 2, tsg_localp) - CALL tsgMakeLocalPolynomialGrid(grid2, dims, outs, 6, 2, tsg_localp_zero) - - N = tsgGetNumNeeded(grid1) - points => tsgGetNeededPoints(grid1) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = cos(0.5 * PI * points(1,i)) & - * cos(0.5 * PI * points(2,i)) - END DO - CALL tsgLoadNeededPoints(grid1, values) - DEALLOCATE(values) - DEALLOCATE(points) - - CALL tsgEvaluateBatch(grid1, randPoints2, 1000, res2) - err1 = 0.0 - DO i = 1, 1000 - IF(abs(res2(1,i) - tvalues(1,i)) .GT. err1)then - err1 = abs(res2(1,i) - tvalues(1,i)) - ENDIF - END DO - - N = tsgGetNumNeeded(grid2) - points => tsgGetNeededPoints(grid2) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = cos(0.5 * PI * points(1,i)) * cos(0.5 * PI * points(2,i)) - END DO - CALL tsgLoadNeededPoints(grid2, values) - DEALLOCATE(values) - DEALLOCATE(points) - - CALL tsgEvaluateBatch(grid2, randPoints2, 1000, res2) - err2 = 0.0 - DO i = 1, 1000 - IF(abs(res2(1,i) - tvalues(1,i)) .GT. err2)then - err2 = abs(res2(1,i) - tvalues(1,i)) - ENDIF - END DO - - N1 = tsgGetNumPoints(grid1) - N2 = tsgGetNumPoints(grid2) - - WRITE(*,"(A,I5,E12.4)") " For rule_localp Number of points: ", N1, err1 - WRITE(*,"(A,I5,E12.4)") " For rule_localp0 Number of points: ", N2, err2 - WRITE(*,*) " Note: localp-zero wins this competition because the function is zero at the boundary" - WRITE(*,*) - -! ==================================================================== ! -! EXAMPLE 9: -! interpolate: f(x,y) = exp(-x) / (1 + 100 * exp(-10 * y)) -! using different refinement schemes - -! remake the true values for the next example - DO i = 1, 1000 - tvalues(1,i) = exp(-randPoints2(1,i)) / (1.0 + 100.0 * exp(-10.0 * randPoints2(2,i))) - END DO - - dims = 2 - outs = 1 - - CALL tsgMakeLocalPolynomialGrid(grid1, dims, outs, 2, -1, tsg_localp) - CALL tsgMakeLocalPolynomialGrid(grid2, dims, outs, 2, -1, tsg_localp) - - N = tsgGetNumNeeded(grid1) - points => tsgGetNeededPoints(grid1) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = exp(-points(1,i)) / (1.0 + 100.0 * exp(-10.0 * points(2,i))) - END DO - CALL tsgLoadNeededPoints(grid1, values) - CALL tsgLoadNeededPoints(grid2, values) - DEALLOCATE(values) - DEALLOCATE(points) - - WRITE(*,*) "-------------------------------------------------------------------------------------------------" - WRITE(*,*) "Example 9: interpolate f(x,y) = exp(-x) / (1 + 100 * exp(-10 * y))" - WRITE(*,*) " the error is estimated as the maximum from 1000 random points" - WRITE(*,*) " tolerance is set at 1.E-5 and maximal order polynomials are used" - WRITE(*,*) - - WRITE(*,*) " Classic FDS" - WRITE(*,*) "iteration points error points error" - - DO j = 1, 7 - ! 1 below corresponds to classic refinement - CALL tsgSetLocalSurplusRefinement(grid1, 1.D-5, tsg_classic) - - N = tsgGetNumNeeded(grid1) - points => tsgGetNeededPoints(grid1) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = exp(-points(1,i)) / (1.0 + 100.0 * exp(-10.0 * points(2,i))) - END DO - CALL tsgLoadNeededPoints(grid1, values) - DEALLOCATE(values) - DEALLOCATE(points) - - CALL tsgEvaluateBatch(grid1, randPoints2, 1000, res2) - err1 = 0.0 - DO i = 1, 1000 - IF(abs(res2(1,i) - tvalues(1,i)) .GT. err1)then - err1 = abs(res2(1,i) - tvalues(1,i)) - ENDIF - END DO - - CALL tsgSetLocalSurplusRefinement(grid2, 1.D-5, tsg_fds) - - N = tsgGetNumNeeded(grid2) - points => tsgGetNeededPoints(grid2) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = exp(-points(1,i)) / (1.0 + 100.0 * exp(-10.0 * points(2,i))) - END DO - CALL tsgLoadNeededPoints(grid2, values) - DEALLOCATE(values) - DEALLOCATE(points) - - CALL tsgEvaluateBatch(grid2, randPoints2, 1000, res2) - err2 = 0.0 - DO i = 1, 1000 - IF(abs(res2(1,i) - tvalues(1,i)) .GT. err2)then - err2 = abs(res2(1,i) - tvalues(1,i)) - ENDIF - END DO - - N1 = tsgGetNumPoints(grid1) - N2 = tsgGetNumPoints(grid2) - - WRITE(*,"(I9,I9,E12.4,I9,E12.4)") j, N1, err1, N2, err2 - END DO - WRITE(*,*) - -! ==================================================================== ! -! EXAMPLE 10: -! interpolate: f(x,y) = exp(-x) / (1 + 100 * exp(-10 * y)) -! using local polynomails and wavelets - - dims = 2 - outs = 1 - - CALL tsgMakeLocalPolynomialGrid(grid1, dims, outs, 3, 1, tsg_localp) - CALL tsgMakeWaveletGrid(grid2, dims, outs, 1, 1) - - N = tsgGetNumNeeded(grid1) - points => tsgGetNeededPoints(grid1) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = exp(-points(1,i)) / (1.0 + 100.0 * exp(-10.0 * points(2,i))) - END DO - CALL tsgLoadNeededPoints(grid1, values) - DEALLOCATE(values) - DEALLOCATE(points) - - N = tsgGetNumNeeded(grid2) - points => tsgGetNeededPoints(grid2) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = exp(-points(1,i)) / (1.0 + 100.0 * exp(-10.0 * points(2,i))) - END DO - CALL tsgLoadNeededPoints(grid2, values) - DEALLOCATE(values) - DEALLOCATE(points) - - WRITE(*,*) "-------------------------------------------------------------------------------------------------" - WRITE(*,*) "Example 10: interpolate f(x,y) = exp(-x) / (1 + 100 * exp(-10 * y))" - WRITE(*,*) " the error is estimated as the maximum from 1000 random points" - WRITE(*,*) " using local polynomials and wavelets" - WRITE(*,*) - - WRITE(*,*) " Polynomials Wavelets" - WRITE(*,*) "iteration points error points error" - - DO j = 1, 8 - CALL tsgSetLocalSurplusRefinement(grid1, 1.D-5, tsg_fds) - - N = tsgGetNumNeeded(grid1) - points => tsgGetNeededPoints(grid1) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = exp(-points(1,i)) / (1.0 + 100.0 * exp(-10.0 * points(2,i))) - END DO - CALL tsgLoadNeededPoints(grid1, values) - DEALLOCATE(values) - DEALLOCATE(points) - - CALL tsgEvaluateBatch(grid1, randPoints2, 1000, res2) - err1 = 0.0 - DO i = 1, 1000 - IF(abs(res2(1,i) - tvalues(1,i)) .GT. err1)then - err1 = abs(res2(1,i) - tvalues(1,i)) - ENDIF - END DO - - CALL tsgSetLocalSurplusRefinement(grid2, 1.D-5, tsg_fds) - - N = tsgGetNumNeeded(grid2) - points => tsgGetNeededPoints(grid2) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = exp(-points(1,i)) / (1.0 + 100.0 * exp(-10.0 * points(2,i))) - END DO - CALL tsgLoadNeededPoints(grid2, values) - DEALLOCATE(values) - DEALLOCATE(points) - - CALL tsgEvaluateBatch(grid2, randPoints2, 1000, res2) - err2 = 0.0 - DO i = 1, 1000 - IF(abs(res2(1,i) - tvalues(1,i)) .GT. err2)then - err2 = abs(res2(1,i) - tvalues(1,i)) - ENDIF - END DO - - N1 = tsgGetNumPoints(grid1) - N2 = tsgGetNumPoints(grid2) - - WRITE(*,"(I9,I9,E12.4,I9,E12.4)") j, N1, err1, N2, err2 - END DO - WRITE(*,*) - - CALL tsgDeallocateGrid(grid1) - CALL tsgDeallocateGrid(grid2) - -! ==================================================================== ! -! EXAMPLE 11: interpolate: f(x,y,z) = 1/((1+4x^2)*(1+5y^2)*(1+6z^2)) -! using classical and conformal transformation - - WRITE(*,*) "-------------------------------------------------------------------------------------------------" - WRITE(*,*) "Example 11: interpolate f(x,y,z) = 1/((1+4x^2)*(1+5y^2)*(1+6z^2))" - WRITE(*,*) " using conformal transformation" - WRITE(*,*) " the error is estimated as the maximum from 1000 random points" - - dims = 3 - outs = 1 - level = 12 - conformal(1) = 4 - conformal(2) = 4 - conformal(3) = 4 - - DO i = 1, 1000 - randPoints3(1,i) = randPoints(1,i) - randPoints3(2,i) = randPoints(2,i) - randPoints3(3,i) = randPoints(3,i) - tvalues(1,i) = 1.0 / ((1.0 + 4.0 * randPoints3(1,i)**2) * & - (1.0 + 5.0 * randPoints3(2,i)**2) * & - (1.0 + 6.0 * randPoints3(3,i)**2)) - END DO - - call tsgAllocateGrid(grid1) - - CALL tsgMakeGlobalGrid(grid1, dims, outs, level, tsg_iptotal, tsg_clenshaw_curtis) - - N = tsgGetNumNeeded(grid1) - points => tsgGetNeededPoints(grid1) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = 1.0 / ((1.0 + 4.0 * points(1,i)**2) * & - (1.0 + 5.0 * points(2,i)**2) * & - (1.0 + 6.0 * points(3,i)**2)) - END DO - CALL tsgLoadNeededPoints(grid1, values) - DEALLOCATE(values) - DEALLOCATE(points) - - CALL tsgEvaluateBatch(grid1, randPoints3, 1000, res2) - err1 = 0.0 - DO i = 1, 1000 - IF(abs(res2(1,i) - tvalues(1,i)) .GT. err1)then - err1 = abs(res2(1,i) - tvalues(1,i)) - ENDIF - END DO - N1 = tsgGetNumPoints(grid1) - - CALL tsgMakeGlobalGrid(grid1, dims, outs, level, tsg_iptotal, tsg_clenshaw_curtis) - CALL tsgSetConformalTransformASIN(grid1, conformal) - - N = tsgGetNumNeeded(grid1) - points => tsgGetNeededPoints(grid1) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = 1.0 / ((1.0 + 4.0 * points(1,i)**2) * & - (1.0 + 5.0 * points(2,i)**2) * & - (1.0 + 6.0 * points(3,i)**2)) - END DO - CALL tsgLoadNeededPoints(grid1, values) - DEALLOCATE(values) - DEALLOCATE(points) - - CALL tsgEvaluateBatch(grid1, randPoints3, 1000, res2) - err2 = 0.0 - DO i = 1, 1000 - IF(abs(res2(1,i) - tvalues(1,i)) .GT. err2)then - err2 = abs(res2(1,i) - tvalues(1,i)) - ENDIF - END DO - - CALL tsgMakeLocalPolynomialGrid(grid1, dims, outs, level-4, 2, tsg_localp) - - N = tsgGetNumNeeded(grid1) - points => tsgGetNeededPoints(grid1) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = 1.0 / ((1.0 + 4.0 * points(1,i)**2) * & - (1.0 + 5.0 * points(2,i)**2) * & - (1.0 + 6.0 * points(3,i)**2)) - END DO - CALL tsgLoadNeededPoints(grid1, values) - DEALLOCATE(values) - DEALLOCATE(points) - - CALL tsgEvaluateBatch(grid1, randPoints3, 1000, res2) - err3 = 0.0 - DO i = 1, 1000 - IF(abs(res2(1,i) - tvalues(1,i)) .GT. err3)then - err3 = abs(res2(1,i) - tvalues(1,i)) - ENDIF - END DO - N2 = tsgGetNumPoints(grid1) - - CALL tsgMakeLocalPolynomialGrid(grid1, dims, outs, level-4, 2, tsg_localp) - CALL tsgSetConformalTransformASIN(grid1, conformal) - - N = tsgGetNumNeeded(grid1) - points => tsgGetNeededPoints(grid1) - ALLOCATE(values(outs,N)) - DO i = 1, N - values(1,i) = 1.0 / ((1.0 + 4.0 * points(1,i)**2) * & - (1.0 + 5.0 * points(2,i)**2) * & - (1.0 + 6.0 * points(3,i)**2)) - END DO - CALL tsgLoadNeededPoints(grid1, values) - DEALLOCATE(values) - DEALLOCATE(points) - - CALL tsgEvaluateBatch(grid1, randPoints3, 1000, res2) - err4 = 0.0 - DO i = 1, 1000 - IF(abs(res2(1,i) - tvalues(1,i)) .GT. err4)then - err4 = abs(res2(1,i) - tvalues(1,i)) - ENDIF - END DO - - WRITE(*,*) "Grid Type nodes error regular error conformal" - WRITE(*,"(A,I8,E18.4,E18.4)") " Global ", N1, err1, err2 - WRITE(*,"(A,I8,E18.4,E18.4)") " Localp ", N2, err3, err4 - WRITE(*,*) - WRITE(*,*) "Note: conformal maps address specific problems with the region of analyticity of a function" - WRITE(*,*) " the map can accelerate or slow down convergence depending on the problem" - WRITE(*,*) - - CALL tsgDeallocateGrid(grid1) -! ==================================================================== ! - -! cleanup - DEALLOCATE(res2) - DEALLOCATE(tvalues) - - WRITE(*,*) "-------------------------------------------------------------------------------------------------" - WRITE(*,*) - -END PROGRAM TasmanianSGExample - diff --git a/InterfaceFortran/TasmanianSG.f90 b/InterfaceFortran/TasmanianSG.f90 deleted file mode 100644 index cd93d624f..000000000 --- a/InterfaceFortran/TasmanianSG.f90 +++ /dev/null @@ -1,827 +0,0 @@ -!================================================================================================================================================================================== -! Copyright (c) 2017, Miroslav Stoyanov -! -! This file is part of -! Toolkit for Adaptive Stochastic Modeling And Non-Intrusive ApproximatioN: TASMANIAN -! -! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: -! -! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. -! -! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions -! and the following disclaimer in the documentation and/or other materials provided with the distribution. -! -! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse -! or promote products derived from this software without specific prior written permission. -! -! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, -! INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. -! IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, -! OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, -! OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -! -! UT-BATTELLE, LLC AND THE UNITED STATES GOVERNMENT MAKE NO REPRESENTATIONS AND DISCLAIM ALL WARRANTIES, BOTH EXPRESSED AND IMPLIED. -! THERE ARE NO EXPRESS OR IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, OR THAT THE USE OF THE SOFTWARE WILL NOT INFRINGE ANY PATENT, -! COPYRIGHT, TRADEMARK, OR OTHER PROPRIETARY RIGHTS, OR THAT THE SOFTWARE WILL ACCOMPLISH THE INTENDED RESULTS OR THAT THE SOFTWARE OR ITS USE WILL NOT RESULT IN INJURY OR DAMAGE. -! THE USER ASSUMES RESPONSIBILITY FOR ALL LIABILITIES, PENALTIES, FINES, CLAIMS, CAUSES OF ACTION, AND COSTS AND EXPENSES, CAUSED BY, RESULTING FROM OR ARISING OUT OF, -! IN WHOLE OR IN PART THE USE, STORAGE OR DISPOSAL OF THE SOFTWARE. -!================================================================================================================================================================================== -Module TasmanianSG -implicit none -PUBLIC ! default include all symbols, list only exceptions as PRIVATE :: - - type TasmanianSparseGrid - double complex :: pntr - end type TasmanianSparseGrid - - integer, parameter :: tsg_clenshaw_curtis = 1, tsg_clenshaw_curtis_zero = 2, & - tsg_chebyshev = 3, tsg_chebyshev_odd = 4, & - tsg_gauss_legendre = 5, tsg_gauss_legendreodd = 6, & - tsg_gauss_patterson = 7, tsg_leja = 8, & - tsg_lejaodd = 9, tsg_rleja = 10, & - tsg_rleja_odd = 11, tsg_rleja_double_2 = 12, & - tsg_rleja_double_4 = 13, tsg_rleja_shifted = 14, & - tsg_rleja_shifted_even = 15, tsg_rleja_shifted_double = 16, & - tsg_max_lebesgue = 17, tsg_max_lebesgue_odd = 18, & - tsg_min_lebesgue = 19, tsg_min_lebesgue_odd = 20, & - tsg_min_delta = 21, tsg_min_delta_odd = 22, & - tsg_gauss_chebyshev_1 = 23, tsg_gauss_chebyshev_1_odd = 24, & - tsg_gauss_chebyshev_2 = 25, tsg_gauss_chebyshev_2_odd = 26, & - tsg_fejer2 = 27, tsg_gauss_gegenbauer = 28, & - tsg_gauss_gegenbauer_odd = 29, tsg_gauss_jacobi = 30, & - tsg_gauss_jacobi_odd = 31, tsg_gauss_laguerre = 32, & - tsg_gauss_laguerre_odd = 33, tsg_gauss_hermite = 34, & - tsg_gauss_hermite_odd = 35, tsg_custom_tabulated = 36, & - tsg_localp = 37, tsg_localp_zero = 38, & - tsg_semi_localp = 39, tsg_wavelet = 40, & - tsg_fourier = 41, tsg_localpb = 42, & - tsg_level = 1, tsg_curved = 2, & - tsg_iptotal = 3, tsg_ipcurved = 4, & - tsg_qptotal = 5, tsg_qpcurved = 6, & - tsg_hyperbolic = 7, tsg_iphyperbolic = 8, & - tsg_qphyperbolic = 9, tsg_tensor = 10, & - tsg_iptensor = 11, tsg_qptensor = 12, & - tsg_classic = 1, tsg_parents_first = 2, & - tsg_directional = 3, tsg_fds = 4, & - tsg_stable = 5, & - tsg_acc_none = 0, tsg_acc_cpu_blas = 1, tsg_acc_gpu_cublas = 2, & - tsg_acc_gpu_rocblas = 2, tsg_acc_gpu_hip = 3, & - tsg_acc_gpu_cuda = 3, tsg_acc_gpu_magma = 4 - - integer :: rows, cols, length - double precision, pointer :: matrix(:,:), vector(:) - character, pointer :: string(:) - -PRIVATE :: rows, cols, length, matrix, vector, string - -contains -!======================================================================= -subroutine tsgAllocateGrid(grid) - type(TasmanianSparseGrid) :: grid - call tsgalloc(grid%pntr); -end subroutine tsgAllocateGrid -subroutine tsgDeallocateGrid(grid) - type(TasmanianSparseGrid) :: grid - call tsgfree(grid%pntr); -end subroutine tsgDeallocateGrid -!======================================================================= -function tsgIsGlobal(grid) result(res) - type(TasmanianSparseGrid) :: grid - integer :: flag - logical :: res - call tsgisg(grid%pntr, flag) - if (flag .ne. 0) then - res = .true. - else - res = .false. - endif -end function -!======================================================================= -function tsgIsSequence(grid) result(res) - type(TasmanianSparseGrid) :: grid - integer :: flag - logical :: res - call tsgiss(grid%pntr, flag) - if (flag .ne. 0) then - res = .true. - else - res = .false. - endif -end function -!======================================================================= -function tsgIsLocalPolynomial(grid) result(res) - type(TasmanianSparseGrid) :: grid - integer :: flag - logical :: res - call tsgisl(grid%pntr, flag) - if (flag .ne. 0) then - res = .true. - else - res = .false. - endif -end function -!======================================================================= -function tsgIsWavelet(grid) result(res) - type(TasmanianSparseGrid) :: grid - integer :: flag - logical :: res - call tsgisw(grid%pntr, flag) - if (flag .ne. 0) then - res = .true. - else - res = .false. - endif -end function -!======================================================================= -function tsgIsFourier(grid) result(res) - type(TasmanianSparseGrid) :: grid - integer :: flag - logical :: res - call tsgisf(grid%pntr, flag) - if (flag .ne. 0) then - res = .true. - else - res = .false. - endif -end function -!======================================================================= -! class TasmanianSparseGrid wrapper functions -!======================================================================= -function tsgGetVersionMajor() result(ver) - integer :: ver - call tsggvm(ver) -end function tsgGetVersionMajor -!======================================================================= -function tsgGetVersionMinor() result(ver) - integer :: ver - call tsggvn(ver) -end function tsgGetVersionMinor -!======================================================================= -function tsgGetLicense() result(lic) - character, pointer :: lic(:) - call tsggli() - lic => string -end function tsgGetLicense -!======================================================================= -subroutine tsgMakeGlobalGrid(grid, dims, outs, depth, gtype, rule, & - aweights, alpha, beta, customRuleFilename, levelLimits) - type(TasmanianSparseGrid) :: grid - integer, intent(in) :: dims, outs, depth, gtype, rule - integer, optional, target :: aweights(:), levelLimits(dims) - double precision, optional :: alpha, beta - character(len=*), optional :: customRuleFilename - integer :: opt_flags(5) - character(len=80) :: cfn = char(0) - double precision :: al, be - integer, pointer :: aw(:) => null() - integer, pointer :: ll(:) => null() - - opt_flags = 0 - if ( present(aweights) ) then - opt_flags(1) = 1 - aw => aweights - endif - if ( present(alpha) ) then - opt_flags(2) = 1 - al = alpha - endif - if ( present(beta) ) then - opt_flags(3) = 1 - be = beta - endif - if ( present(customRuleFilename) ) then - opt_flags(4) = 1 - cfn = customRuleFilename//char(0) - endif - if ( present(levelLimits) ) then - opt_flags(5) = 1 - ll => levelLimits - endif - - call tsgmg(grid%pntr, dims, outs, depth, gtype, rule, opt_flags, aw, al, be, cfn, ll) -end subroutine tsgMakeGlobalGrid -!======================================================================= -subroutine tsgMakeSequenceGrid(grid, dims, outs, depth, gtype, rule, aweights, levelLimits) - type(TasmanianSparseGrid) :: grid - integer :: dims, outs, depth, gtype, rule - integer, optional, target :: aweights(:), levelLimits(dims) - integer :: opt_flags(2) - integer, pointer :: aw(:) => null() - integer, pointer :: ll(:) => null() - - opt_flags = 0 - if ( present(aweights) ) then - opt_flags(1) = 1 - aw => aweights - endif - if ( present(levelLimits) ) then - opt_flags(2) = 1 - ll => levelLimits - endif - - call tsgms(grid%pntr, dims, outs, depth, gtype, rule, opt_flags, aw, ll) -end subroutine tsgMakeSequenceGrid -!======================================================================= -subroutine tsgMakeLocalPolynomialGrid(grid, dims, outs, depth, order, rule, levelLimits) - type(TasmanianSparseGrid) :: grid - integer :: dims, outs, depth - integer, optional :: order, rule - integer, optional, target :: levelLimits(dims) - integer :: opt_flags(3), or, ru - integer, pointer :: ll(:) => null() - - opt_flags = 0 - if ( present(order) ) then - opt_flags(1) = 1 - or = order - endif - if ( present(rule) ) then - opt_flags(2) = 1 - ru = rule - endif - if ( present(levelLimits) ) then - opt_flags(3) = 1 - ll => levelLimits - endif - - call tsgml(grid%pntr, dims, outs, depth, opt_flags, or, ru, ll) -end subroutine tsgMakeLocalPolynomialGrid -!======================================================================= -subroutine tsgMakeWaveletGrid(grid, dims, outs, depth, order, levelLimits) - type(TasmanianSparseGrid) :: grid - integer :: dims, outs, depth - integer, optional, target :: levelLimits(dims) - integer, optional :: order - integer :: opt_flags(2), or - integer, pointer :: ll(:) => null() - - opt_flags = 0 - if ( present(order) ) then - opt_flags(1) = 1 - or = order - endif - if ( present(levelLimits) ) then - opt_flags(2) = 1 - ll => levelLimits - endif - - call tsgmw(grid%pntr, dims, outs, depth, opt_flags, or, ll) -end subroutine tsgMakeWaveletGrid -!======================================================================= -subroutine tsgMakeFourierGrid(grid, dims, outs, depth, gtype, aweights, levelLimits) - type(TasmanianSparseGrid) :: grid - integer :: dims, outs, depth, gtype - integer, optional, target :: aweights(:), levelLimits(dims) - integer :: opt_flags(2) - integer, pointer :: aw(:) => null() - integer, pointer :: ll(:) => null() - - opt_flags = 0 - if ( present(aweights) ) then - opt_flags(1) = 1 - aw => aweights - endif - if ( present(levelLimits) ) then - opt_flags(2) = 1 - ll => levelLimits - endif - - call tsgmf(grid%pntr, dims, outs, depth, gtype, opt_flags, aw, ll) -end subroutine tsgMakeFourierGrid -!======================================================================= -subroutine tsgCopyGrid(grid, source) - type(TasmanianSparseGrid) :: grid, source - call tsgcp(grid%pntr, source) -end subroutine tsgCopyGrid -!======================================================================= -subroutine tsgUpdateGlobalGrid(grid, depth, gtype, aweights) - type(TasmanianSparseGrid) :: grid - integer, intent(in) :: depth, gtype - integer, optional, target :: aweights(:) - integer :: opt_flags = 0 - integer, pointer :: aw(:) => null() - if ( present(aweights) ) then - opt_flags = 1 - aw => aweights - endif - call tsgug(grid%pntr, depth, gtype, opt_flags, aw) -end subroutine tsgUpdateGlobalGrid -!======================================================================= -subroutine tsgUpdateSequenceGrid(grid, depth, gtype, aweights) - type(TasmanianSparseGrid) :: grid - integer, intent(in) :: depth, gtype - integer, optional, target :: aweights(:) - integer :: opt_flags = 0 - integer, pointer :: aw(:) => null() - if ( present(aweights) ) then - opt_flags = 1 - aw => aweights - endif - call tsgus(grid%pntr, depth, gtype, opt_flags, aw) -end subroutine tsgUpdateSequenceGrid -!======================================================================= -function tsgRead(grid, filename) result(res) - type(TasmanianSparseGrid) :: grid - character(len=*), intent(in) :: filename - logical :: res - integer :: N, i, flag - character, allocatable :: cfilename(:) - N = len(trim(filename)) - allocate(cfilename(N+1)) - do i = 1, N - cfilename(i) = filename(i:i) - end do - cfilename(N+1) = CHAR(0) - call tsgrea(grid%pntr, cfilename, flag) - res = (flag .ne. 0) - deallocate(cfilename) -end function tsgRead -!======================================================================= -subroutine tsgWrite(grid, filename, useBinary) - type(TasmanianSparseGrid) :: grid - logical, intent(in), optional :: useBinary - integer :: ubin - character(len=*), intent(in) :: filename - integer :: N, i - character, allocatable :: cfilename(:) - N = len(trim(filename)) - allocate(cfilename(N+1)) - do i = 1, N - cfilename(i) = filename(i:i) - end do - cfilename(N+1) = CHAR(0) - if(present(useBinary))then - if(useBinary)then - ubin = 1 - else - ubin = 0 - endif - else - ubin = 1 - endif - call tsgwri(grid%pntr, cfilename, ubin) - deallocate(cfilename) -end subroutine tsgWrite -!======================================================================= -function tsgGetAlpha(grid) result(alpha) - type(TasmanianSparseGrid) :: grid - double precision :: alpha - call tsggal(grid%pntr, alpha) -end function tsgGetAlpha -!======================================================================= -function tsgGetBeta(grid) result(beta) - type(TasmanianSparseGrid) :: grid - double precision :: beta - call tsggbe(grid%pntr, beta) -end function tsgGetBeta -!======================================================================= -function tsgGetOrder(grid) result(order) - type(TasmanianSparseGrid) :: grid - integer :: order - call tsggor(grid%pntr, order) -end function tsgGetOrder -!======================================================================= -function tsgGetNumDimensions(grid) result(dims) - type(TasmanianSparseGrid) :: grid - integer :: dims - call tsggnd(grid%pntr, dims) -end function tsgGetNumDimensions -!======================================================================= -function tsgGetNumOutputs(grid) result(outs) - type(TasmanianSparseGrid) :: grid - integer :: outs - call tsggno(grid%pntr, outs) -end function tsgGetNumOutputs -!======================================================================= -function tsgGetRule(grid) result(rule) - type(TasmanianSparseGrid) :: grid - integer :: rule - call tsggru(grid%pntr, rule) -end function tsgGetRule -!======================================================================= -function tsgGetNumLoaded(grid) result(num) - type(TasmanianSparseGrid) :: grid - integer :: num - call tsggnl(grid%pntr, num) -end function tsgGetNumLoaded -!======================================================================= -function tsgGetNumNeeded(grid) result(num) - type(TasmanianSparseGrid) :: grid - integer :: num - call tsggnn(grid%pntr, num) -end function tsgGetNumNeeded -!======================================================================= -function tsgGetNumPoints(grid) result(num) - type(TasmanianSparseGrid) :: grid - integer :: num - call tsggnp(grid%pntr, num) -end function tsgGetNumPoints -!======================================================================= -function tsgGetLoadedPoints(grid) result(p) - type(TasmanianSparseGrid) :: grid - double precision, pointer :: p(:,:) - integer :: rows, cols - rows = tsgGetNumDimensions(grid) - cols = tsgGetNumLoaded(grid) - allocate(p(rows,cols)) - call tsgglp(grid%pntr, p) -end function tsgGetLoadedPoints -!======================================================================= -function tsgGetNeededPoints(grid) result(p) - type(TasmanianSparseGrid) :: grid - double precision, pointer :: p(:,:) - integer :: rows, cols - rows = tsgGetNumDimensions(grid) - cols = tsgGetNumNeeded(grid) - allocate(p(rows,cols)) - call tsggdp(grid%pntr, p) -end function tsgGetNeededPoints -!======================================================================= -function tsgGetPoints(grid) result(p) - type(TasmanianSparseGrid) :: grid - double precision, pointer :: p(:,:) - integer :: rows, cols - rows = tsgGetNumDimensions(grid) - cols = tsgGetNumPoints(grid) - allocate(p(rows,cols)) - call tsggpp(grid%pntr, p) -end function tsgGetPoints -!======================================================================= -subroutine tsgGetLoadedPointsStatic(grid, points) - type(TasmanianSparseGrid) :: grid - double precision :: points(:,:) - call tsgglp(grid%pntr, points) -end subroutine tsgGetLoadedPointsStatic -!======================================================================= -subroutine tsgGetNeededPointsStatic(grid, points) - type(TasmanianSparseGrid) :: grid - double precision :: points(:,:) - call tsggdp(grid%pntr, points) -end subroutine tsgGetNeededPointsStatic -!======================================================================= -subroutine tsgGetPointsStatic(grid, points) - type(TasmanianSparseGrid) :: grid - double precision :: points(:,:) - call tsggpp(grid%pntr, points) -end subroutine tsgGetPointsStatic -!======================================================================= -function tsgGetQuadratureWeights(grid) result(w) - type(TasmanianSparseGrid) :: grid - integer :: length - double precision, pointer :: w(:) - length = tsgGetNumPoints(grid) - allocate(w(length)) - call tsggqw(grid%pntr, w) -end function tsgGetQuadratureWeights -!======================================================================= -subroutine tsgGetQuadratureWeightsStatic(grid, weights) - type(TasmanianSparseGrid) :: grid - double precision :: weights(*) - call tsggqw(grid%pntr, weights) -end subroutine tsgGetQuadratureWeightsStatic -!======================================================================= -function tsgGetInterpolationWeights(grid, x) result(w) - type(TasmanianSparseGrid) :: grid - integer :: length - double precision :: x(*) - double precision, pointer :: w(:) - length = tsgGetNumPoints(grid) - allocate(w(length)) - call tsggiw(grid%pntr, x, w) -end function tsgGetInterpolationWeights -!======================================================================= -subroutine tsgGetInterpolationWeightsStatic(grid, x, weights) - type(TasmanianSparseGrid) :: grid - double precision :: x(:) - double precision :: weights(:) - call tsggiw(grid%pntr, x, weights) -end subroutine tsgGetInterpolationWeightsStatic -!======================================================================= -subroutine tsgLoadNeededPoints(grid, values) - type(TasmanianSparseGrid) :: grid - double precision :: values(:,:) - call tsglnp(grid%pntr, values) -end subroutine tsgLoadNeededPoints -!======================================================================= -subroutine tsgEvaluate(grid, x, y) - type(TasmanianSparseGrid) :: grid - double precision, intent(in) :: x(:) - double precision :: y(:) - call tsgeva(grid%pntr, x, y) -end subroutine tsgEvaluate -!======================================================================= -subroutine tsgEvaluateFast(grid, x, y) - type(TasmanianSparseGrid) :: grid - double precision, intent(in) :: x(:) - double precision :: y(:) - call tsgevf(grid%pntr, x, y) -end subroutine tsgEvaluateFast -!======================================================================= -subroutine tsgEvaluateBatch(grid, x, numX, y) - type(TasmanianSparseGrid) :: grid - integer :: numX - double precision :: x(:,:), y(:,:) - call tsgevb(grid%pntr, x, numX, y) -end subroutine tsgEvaluateBatch -!======================================================================= -subroutine tsgEvaluateHierarchicalFunctions(grid, x, numX, y) - type(TasmanianSparseGrid) :: grid - integer :: numX - double precision :: x(:,:), y(:,:) - if ( .not. tsgIsFourier(grid) ) then - call tsgehf(grid%pntr, x, numX, y) - else - write(*,*) "ERROR: called tsgEvaluateHierarchicalFunctions() on a Fourier grid, " - write(*,*) " use tsgEvaluateComplexHierarchicalFunctions() instead" - endif -end subroutine tsgEvaluateHierarchicalFunctions -!======================================================================= -subroutine tsgEvaluateComplexHierarchicalFunctions(grid, x, numX, y) - type(TasmanianSparseGrid) :: grid - integer :: numX - double precision :: x(:,:) - double complex :: y(:,:) - double precision, allocatable :: y_c_style(:) - integer :: i, j - if ( tsgIsFourier(grid) ) then - allocate(y_c_style(2*size(y,1)*size(y,2))) - call tsgehf(grid%pntr, x, numX, y_c_style) - do j = 1,size(y,2) - do i = 1,size(y,1) - y(i,j) = cmplx( y_c_style( 2*size(y,1)*(j-1) + 2*i-1 ), y_c_style( 2*size(y,1)*(j-1) + 2*i), kind(y) ) - enddo - enddo - else - write(*,*) "ERROR: called tsgEvaluateComplexHierarchicalFunctions() on a non-Fourier grid, " - write(*,*) " use tsgEvaluateHierarchicalFunctions() instead" - endif -end subroutine tsgEvaluateComplexHierarchicalFunctions -!======================================================================= -subroutine tsgEvaluateSparseHierarchicalFunctions(grid, x, numX, pntr, indx, y) - type(TasmanianSparseGrid) :: grid - integer :: numX - double precision :: x(:,:) - integer, pointer :: pntr(:), indx(:) - double precision, pointer :: y(:) - integer :: numNZ - if ( .not. tsgIsFourier(grid) ) then - call tsgehz(grid%pntr, x, numX, numNZ) - allocate( pntr(numX+1), indx(numNZ), y(numNZ) ) - call tsgehs(grid%pntr, x, numX, pntr, indx, y) - else - write(*,*) "ERROR: called tsgEvaluateSparseHierarchicalFunctions() on a Fourier grid" - endif -end subroutine tsgEvaluateSparseHierarchicalFunctions -!======================================================================= -function tsgGetHierarchicalCoefficients(grid) result(c) - type(TasmanianSparseGrid) :: grid - double precision, pointer :: c(:) - if ( .not. tsgIsFourier(grid) ) then - allocate(c(tsgGetNumOutputs(grid)*tsgGetNumPoints(grid))) - call tsgghc(grid%pntr, c) - else - write(*,*) "ERROR: called tsgGetHierarchicalCoefficients() on a Fourier grid, " - write(*,*) " use tsgGetComplexHierarchicalCoefficients() instead" - endif -end function tsgGetHierarchicalCoefficients -!======================================================================= -subroutine tsgGetHierarchicalCoefficientsStatic(grid, c) - type(TasmanianSparseGrid) :: grid - double precision :: c(:) - if ( .not. tsgIsFourier(grid) ) then - call tsgghc(grid%pntr, c) - else - write(*,*) "ERROR: called tsgGetHierarchicalCoefficientsStatic() on a Fourier grid, " - write(*,*) " use tsgGetComplexHierarchicalCoefficientsStatic() instead" - endif -end subroutine tsgGetHierarchicalCoefficientsStatic -!======================================================================= -function tsgGetComplexHierarchicalCoefficients(grid) result(c) - type(TasmanianSparseGrid) :: grid - double complex, pointer :: c(:) - double precision, allocatable :: c_real(:) - integer :: i - if ( tsgIsFourier(grid) ) then - allocate(c(tsgGetNumOutputs(grid)*tsgGetNumPoints(grid))) - allocate(c_real(2*tsgGetNumOutputs(grid)*tsgGetNumPoints(grid))) - call tsgghc(grid%pntr, c_real) - do i = 1,size(c) - c(i) = cmplx( c_real(i), c_real(i + size(c)), kind(c) ) - enddo - deallocate(c_real) - else - write(*,*) "ERROR: called tsgGetComplexHierarchicalCoefficients() on a non-Fourier grid, " - write(*,*) " use tsgGetHierarchicalCoefficients() instead" - endif -end function tsgGetComplexHierarchicalCoefficients -!======================================================================= -subroutine tsgGetComplexHierarchicalCoefficientsStatic(grid, c) - type(TasmanianSparseGrid) :: grid - double complex :: c(:) - double precision :: c_real(2*size(c)) - integer :: i - if ( tsgIsFourier(grid) ) then - call tsgghc(grid%pntr, c_real) - do i = 1,size(c) - c(i) = cmplx( c_real(i), c_real(i + size(c)), kind(c) ) - enddo - else - write(*,*) "ERROR: called tsgGetComplexHierarchicalCoefficientsStatic() on a non-Fourier grid, " - write(*,*) " use tsgGetHierarchicalCoefficientsStatic() instead" - endif -end subroutine tsgGetComplexHierarchicalCoefficientsStatic -!======================================================================= -subroutine tsgSetHierarchicalCoefficients(grid,c) - type(TasmanianSparseGrid) :: grid - double precision :: c(:) - call tsgshc(grid%pntr, c) -end subroutine tsgSetHierarchicalCoefficients -!======================================================================= -function tsgGetHierarchicalSupport(grid) result(c) - type(TasmanianSparseGrid) :: grid - double precision, pointer :: c(:,:) - allocate(c(tsgGetNumDimensions(grid), tsgGetNumPoints(grid))) - call tsghsu(grid%pntr, c) -end function tsgGetHierarchicalSupport -!======================================================================= -subroutine tsgIntegrate(grid, q) - type(TasmanianSparseGrid) :: grid - double precision :: q(:) - call tsgint(grid%pntr, q) -end subroutine tsgIntegrate -!======================================================================= -subroutine tsgSetDomainTransform(grid, transformA, transformB) - type(TasmanianSparseGrid) :: grid - double precision :: transformA(*), transformB(*) - call tsgsdt(grid%pntr, transformA, transformB) -end subroutine tsgSetDomainTransform -!======================================================================= -function tsgIsSetDomainTransform(grid) result(res) - type(TasmanianSparseGrid) :: grid - logical :: res - call tsgidt(grid%pntr, res) -end function tsgIsSetDomainTransform -!======================================================================= -subroutine tsgClearDomainTransform(grid) - type(TasmanianSparseGrid) :: grid - call tsgcdt(grid%pntr) -end subroutine tsgClearDomainTransform -!======================================================================= -subroutine tsgGetDomainTransform(grid, transformA, transformB) - type(TasmanianSparseGrid) :: grid - double precision :: transformA(*), transformB(*) - call tsggdt(grid%pntr, transformA, transformB) -end subroutine tsgGetDomainTransform -!======================================================================= -subroutine tsgSetAnisotropicRefinement(grid, gtype, minGrowth, output, levelLimits) - type(TasmanianSparseGrid) :: grid - integer :: gtype, minGrowth, output - integer, optional, target :: levelLimits(:) - integer :: opt_flags = 0 - integer, pointer :: ll(:) => null() - if (present(levelLimits)) then - opt_flags = 1 - ll => levelLimits - endif - call tsgsar(grid%pntr, gtype, minGrowth, output-1, opt_flags, ll) -end subroutine tsgSetAnisotropicRefinement -!======================================================================= -function tsgEstimateAnisotropicCoefficients(grid, gtype, output) result(coeff) - type(TasmanianSparseGrid) :: grid - integer :: gtype, output, N - integer, pointer :: coeff(:) - N = tsgGetNumDimensions(grid) - if ((gtype .EQ. tsg_curved) .OR. (gtype .EQ. tsg_ipcurved) .OR. (gtype .EQ. tsg_qpcurved))then - N = N * 2 - endif - allocate(coeff(N)) - call tsgeac(grid%pntr, gtype, output-1, coeff) -end function tsgEstimateAnisotropicCoefficients -!======================================================================= -subroutine tsgSetGlobalSurplusRefinement(grid, tolerance, output, levelLimits) - type(TasmanianSparseGrid) :: grid - integer :: output - integer, optional, target :: levelLimits(:) - double precision :: tolerance - integer :: opt_flags = 0 - integer, pointer :: ll(:) => null() - if (present(levelLimits)) then - opt_flags = 1 - ll => levelLimits - endif - call tsgssr(grid%pntr, tolerance, output-1, opt_flags, ll) -end subroutine tsgSetGlobalSurplusRefinement -!======================================================================= -subroutine tsgSetLocalSurplusRefinement(grid, tolerance, rtype, output, levelLimits) - type(TasmanianSparseGrid) :: grid - integer :: rtype - integer, optional :: output - integer, optional, target :: levelLimits(:) - double precision :: tolerance - integer :: opt_flags(2), theout - integer, pointer :: ll(:) => null() - - opt_flags = 0 - if (present(output)) then - opt_flags(1) = 1 - theout = output-1 - endif - if (present(levelLimits)) then - opt_flags(2) = 1 - ll => levelLimits - endif - - call tsgshr(grid%pntr, tolerance, rtype, opt_flags, theout, ll) -end subroutine tsgSetLocalSurplusRefinement -!======================================================================= -subroutine tsgClearRefinement(grid) - type(TasmanianSparseGrid) :: grid - call tsgcre(grid%pntr) -end subroutine tsgClearRefinement -!======================================================================= -subroutine tsgMergeRefinement(grid) - type(TasmanianSparseGrid) :: grid - call tsgmre(grid%pntr) -end subroutine tsgMergeRefinement -!======================================================================= -subroutine tsgSetConformalTransformASIN(grid, truncate) - type(TasmanianSparseGrid) :: grid - integer :: truncate(:) - call tsgsca(grid%pntr, truncate) -end subroutine tsgSetConformalTransformASIN -!======================================================================= -function tsgIsSetConformalTransformASIN(grid) result(isset) - type(TasmanianSparseGrid) :: grid - integer :: res - logical :: isset - call tsgica(grid%pntr, res) - if (res .EQ. 0)then - isset = .FALSE. - else - isset = .TRUE. - endif -end function tsgIsSetConformalTransformASIN -!======================================================================= -subroutine tsgClearConformalTransform(grid) - type(TasmanianSparseGrid) :: grid - call tsgcct(grid%pntr) -end subroutine tsgClearConformalTransform -!======================================================================= -function tsgGetConformalTransformASIN(grid) result(truncate) - type(TasmanianSparseGrid) :: grid - integer :: N - integer, allocatable :: truncate(:) - N = tsgGetNumDimensions(grid) - allocate(truncate(N)) - call tsggca(grid%pntr, truncate) -end function tsgGetConformalTransformASIN -!======================================================================= -subroutine tsgPrintStats(grid) - type(TasmanianSparseGrid) :: grid - call tsgpri(grid%pntr) -end subroutine tsgPrintStats -!======================================================================= -! Addon/MPI methods ! -!======================================================================= -subroutine tsgMPIGridSend(grid, destination, tag, comm, ierr) - type(TasmanianSparseGrid) :: grid - integer :: destination, tag, comm, ierr - call tsgmpigsend(grid, destination, tag, comm, ierr) -end subroutine tsgMPIGridSend -subroutine tsgMPIGridRecv(grid, source, tag, comm, stats, ierr) - type(TasmanianSparseGrid) :: grid - integer :: source, tag, comm, stats(:), ierr - call tsgmpigrecv(grid, source, tag, comm, ierr) -end subroutine tsgMPIGridRecv -subroutine tsgMPIGridBcast(grid, root, comm, ierr) - type(TasmanianSparseGrid) :: grid - integer :: root, comm, ierr - call tsgmpigcast(grid, root, comm, ierr) -end subroutine tsgMPIGridBcast -!======================================================================= -! do NOT call THOSE FUNCTIONS DIRECTLY ! -!======================================================================= -subroutine tsgReceiveString(l, S) - integer :: l - character, target :: S(l) - length = l - string => S(1:length) -end subroutine tsgReceiveString -!======================================================================= -subroutine tsgReceiveVector(s, V) - integer :: s - double precision, target :: V(s) - length = s - vector => V(1:length) -end subroutine tsgReceiveVector -!======================================================================= -subroutine tsgReceiveMatrix(r, c, M) - integer :: r, c - double precision, target :: M(r,c) - rows = r - cols = c - matrix => M(1:rows,1:cols) -end subroutine tsgReceiveMatrix -!======================================================================= -end MODULE diff --git a/InterfaceFortran/fortester.f90 b/InterfaceFortran/fortester.f90 deleted file mode 100644 index f71b2ffab..000000000 --- a/InterfaceFortran/fortester.f90 +++ /dev/null @@ -1,1594 +0,0 @@ -!================================================================================================================================================================================== -! Copyright (c) 2018, Miroslav Stoyanov -! -! This file is part of -! Toolkit for Adaptive Stochastic Modeling And Non-Intrusive ApproximatioN: TASMANIAN -! -! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: -! -! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. -! -! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions -! and the following disclaimer in the documentation and/or other materials provided with the distribution. -! -! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse -! or promote products derived from this software without specific prior written permission. -! -! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, -! INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. -! IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, -! OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, -! OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -! -! UT-BATTELLE, LLC AND THE UNITED STATES GOVERNMENT MAKE NO REPRESENTATIONS AND DISCLAIM ALL WARRANTIES, BOTH EXPRESSED AND IMPLIED. -! THERE ARE NO EXPRESS OR IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, OR THAT THE USE OF THE SOFTWARE WILL NOT INFRINGE ANY PATENT, -! COPYRIGHT, TRADEMARK, OR OTHER PROPRIETARY RIGHTS, OR THAT THE SOFTWARE WILL ACCOMPLISH THE INTENDED RESULTS OR THAT THE SOFTWARE OR ITS USE WILL NOT RESULT IN INJURY OR DAMAGE. -! THE USER ASSUMES RESPONSIBILITY FOR ALL LIABILITIES, PENALTIES, FINES, CLAIMS, CAUSES OF ACTION, AND COSTS AND EXPENSES, CAUSED BY, RESULTING FROM OR ARISING OUT OF, -! IN WHOLE OR IN PART THE USE, STORAGE OR DISPOSAL OF THE SOFTWARE. -!================================================================================================================================================================================== -program FORTESTER - USE TasmanianSG - implicit none - -integer, parameter :: seed = 86456 -double precision, parameter :: pi = 4.0D+0 * atan(1.0D+0) - -character, pointer :: licence(:) -character(14), parameter :: filename = 'fortranio.test' - -type(TasmanianSparseGrid) :: grid, grid_II - -double precision, pointer :: rnd(:,:) -double precision, pointer :: points(:,:), weights(:), pointsb(:,:), weightsb(:) - - -integer :: i_a, i_b, i_c, i_d, i, j -double precision :: d_a, d_b, d_c, d_d, int, int2 -logical :: bool - -integer, allocatable :: int_1d_a(:), int_1d_b(:), int_1d_c(:), int_1d_d(:) -double precision, allocatable :: double_1d_a(:), double_1d_b(:), double_1d_c(:), double_1d_d(:) -double complex, allocatable :: dcmplx_1d_a(:), dcmplx_1d_b(:), dcmplx_1d_c(:), dcmplx_1d_d(:) - -integer, allocatable :: int_2d_a(:,:), int_2d_b(:,:), int_2d_c(:,:), int_2d_d(:,:) -double precision, allocatable :: double_2d_a(:,:), double_2d_b(:,:), double_2d_c(:,:), double_2d_d(:,:) -double complex, allocatable :: dcmplx_2d_a(:,:), dcmplx_2d_b(:,:), dcmplx_2d_c(:,:), dcmplx_2d_d(:,:) - -integer, pointer :: int_pnt_1d_a(:), int_pnt_1d_b(:), int_pnt_1d_c(:) -double precision, pointer :: double_pnt_1d_a(:), double_pnt_1d_b(:), double_pnt_1d_c(:) -double complex, pointer :: dcmplx_pnt_1d_a(:), dcmplx_pnt_1d_b(:), dcmplx_pnt_1d_c(:) - - - -write(*,'(a)') 'Testing TASMANIAN FORTRAN interface' -write(*,*) - -! read library meta-data -licence => tsgGetLicense() -write(*,"(A,I2,A,I1)") "Tasmanian Sparse Grid module: ", tsgGetVersionMajor(), ".", tsgGetVersionMinor() -write(*,"(A,40A)") "Licence: ", licence - - -! initialize random number generator -call srand(seed) - - -call tsgAllocateGrid(grid) -call tsgAllocateGrid(grid_II) - - -!======================================================================= -! tsgIs***() -!======================================================================= -! test type of grid -call tsgMakeGlobalGrid(grid, 2, 0, 1, tsg_level, tsg_clenshaw_curtis) -if ( (.not. tsgIsGlobal(grid)) .or. tsgIsSequence(grid) .or. tsgIsLocalPolynomial(grid) & - .or. tsgIsFourier(grid) .or. tsgIsWavelet(grid) ) then - write(*,*) "Mismatch in tsgIsGlobal" - stop 1 -endif -call tsgMakeSequenceGrid(grid, 2, 1, 3, tsg_level, tsg_min_lebesgue) -if ( .not. tsgIsSequence(grid) ) then - write(*,*) "Mismatch in tsgIsSequence" - stop 1 -endif -call tsgMakeLocalPolynomialGrid(grid, 3, 1, 2, 2, tsg_localp) -if ( (.not. tsgIsLocalPolynomial(grid)) .or. tsgIsSequence(grid) .or. tsgIsGlobal(grid) & - .or. tsgIsFourier(grid) .or. tsgIsWavelet(grid) ) then - write(*,*) "Mismatch in tsgIsLocalPolynomial" - stop 1 -endif -call tsgMakeWaveletGrid(grid, 3, 1, 2, 1) -if ( .not. tsgIsWavelet(grid) ) then - write(*,*) "Mismatch in tsgIsWavelet" - stop 1 -endif -call tsgMakeFourierGrid(grid, 3, 1, 1, tsg_level) -if ( .not. tsgIsFourier(grid) ) then - write(*,*) "Mismatch in tsgIsFourier" - stop 1 -endif -! test transform -if ( tsgIsSetDomainTransform(grid) ) then - write(*,*) "Mismatch in tsgIsSetDomainTransform" - stop 1 -endif -call tsgMakeGlobalGrid(grid, 3, 0, 2, tsg_level, tsg_clenshaw_curtis) -call tsgSetDomainTransform(grid, transformA=(/ 3.0d0, -7.0d0, -12.0d0 /), transformB=(/ 5.0d0, -6.0d0, 17.0d0 /)) -points => tsgGetPoints(grid) -weights => tsgGetQuadratureWeights(grid) -if ( .not. tsgIsSetDomainTransform(grid) ) then - write(*,*) "Mismatch in tsgIsSetDomainTransform" - stop 1 -endif -deallocate(points, weights) - - - - - - - -!======================================================================= -! tsgMakeGlobalGrid() -!======================================================================= -! Check for correct points and weights -allocate( pointsb(2,5), weightsb(5) ) -pointsb = reshape((/0.0D+0, 0.0D+0, 0.0D+0, -1.0D+0, 0.0D+0, 1.0D+0, -1.0D+0, 0.0D+0, 1.0D+0, 0.0D+0/), shape(pointsb)) -weightsb = reshape((/4.0D+0/3.0D+0, 2.0D+0/3.0D+0, 2.0D+0/3.0D+0, 2.0D+0/3.0D+0, 2.0D+0/3.0D+0/), shape(weightsb)) -call tsgMakeGlobalGrid(grid, 2, 0, 1, tsg_level, tsg_clenshaw_curtis) -points => tsgGetPoints(grid) -weights => tsgGetQuadratureWeights(grid) -if ( norm(pointsb-points,size(points)) > 1.d-11 .OR. norm1d(weightsb-weights) > 1.d-11 ) then - write(*,*) "Mismatch in tsgMakeGlobal: core case 1", sum(abs(pointsb-points)), sum(abs(weightsb-weights)) - stop 1 -end if -deallocate( pointsb, weightsb, points, weights ) - - -call tsgMakeGlobalGrid(grid, 2, 0, 2, tsg_level, tsg_clenshaw_curtis) -points => tsgGetPoints(grid) -weights => tsgGetQuadratureWeights(grid) -if ((abs(points(2, 4) + sqrt(0.5D+0)) > 1.0D-11) .OR. & - (abs(points(2, 5) - sqrt(0.5D+0)) > 1.0D-11) .OR. & - (abs(weights(7) - 1.0D+0/9.0D+0) > 1.0D-11))then - write(*,*) "Mismatch in tsgMakeGlobal: core case 2", & - abs(points(2, 4) - sqrt(0.5D+0)), abs(points(2, 5) - sqrt(0.5D+0)), abs(weights(7) - 1.0D+0/9.0D+0) - stop 1 -end if -deallocate(points, weights) - - -call tsgMakeGlobalGrid(grid, 3, 0, 4, tsg_level, tsg_fejer2) -points => tsgGetPoints(grid) -weights => tsgGetQuadratureWeights(grid) -if ( (abs(sum(points)) > 1.0D-11) .OR. (abs(sum(weights)-8.0D+0) > 1.0D-11) )then - write(*,*) "Mismatch in tsgMakeGlobal: core case 3" - write(*,*) abs(sum(points)), abs(sum(weights) - 8.0D+0) - stop 1 -end if -deallocate(points, weights) - - -allocate( pointsb(1,4), weightsb(4) ) -pointsb(1,:) = (/0.0D+0, 1.0D+0, -1.0D+0, sqrt(1.0D+0/3.0D+0)/) -weightsb = (/4.0D+0/3.0D+0, 1.0D+0/3.0D+0, 1.0D+0/3.0D+0, 0.0D+0/) -call tsgMakeGlobalGrid(grid, 1, 0, 3, tsg_level, tsg_leja) -points => tsgGetPoints(grid) -weights => tsgGetQuadratureWeights(grid) -if ( (norm2d(pointsb-points) > 1.0D-11) .OR. (norm1d(weightsb-weights) > 1.0D-11) )then - write(*,*) "Mismatch in tsgMakeGlobal: core case 4", norm2d(pointsb-points), norm1d(weightsb-weights) - stop 1 -end if -deallocate(pointsb, weightsb, points, weights) - - -! test transform -call tsgMakeGlobalGrid(grid, 3, 0, 2, tsg_level, tsg_clenshaw_curtis) -call tsgSetDomainTransform(grid, transformA=(/ 3.0d0, -7.0d0, -12.0d0 /), transformB=(/ 5.0d0, -6.0d0, 17.0d0 /)) -points => tsgGetPoints(grid) -weights => tsgGetQuadratureWeights(grid) -if ( (abs(sum(weights)-58.d0) > 1.d-11) .or. & - (abs(maxval(points(1,:))-5.d0) > 1.d-11) .or. (abs(minval(points(1,:))-3.d0) > 1.d-11) .or. & - (abs(maxval(points(2,:))+6.d0) > 1.d-11) .or. (abs(minval(points(2,:))+7.d0) > 1.d-11) .or. & - (abs(maxval(points(3,:))-17.d0) > 1.d-11) .or. (abs(minval(points(3,:))+12.d0) > 1.d-11) ) then - write(*,*) "Mismatch in tsgMakeGlobal: transform " - stop 1 -endif -deallocate(points, weights) - - -! test alpha/beta -call tsgMakeGlobalGrid(grid, 1, 0, 4, tsg_level, tsg_gauss_hermite, alpha=2.0D+0) -weights => tsgGetQuadratureWeights(grid) -if ( (abs(sum(weights)-0.5d0*sqrt(pi)) > 1.D-11) .OR. (abs(tsgGetBeta(grid)) > 1.D-11) & - .OR. (abs(tsgGetAlpha(grid) - 2.0D+0) > 1.D-11) .OR. (tsgGetOrder(grid) /= -1)) then - write(*,*) "Mismatch in tsgMakeGlobal: alpha/beta", sum(weights), tsgGetBeta(grid) - stop 1 -endif -deallocate(weights) - - -! test anisotropy -allocate( pointsb(2,4) ) -pointsb = reshape((/0.0D+0, 0.0D+0, 0.0D+0, 1.0D+0, 0.0D+0, -1.0D+0, 1.0D+0, 0.0D+0/), shape(pointsb)) -call tsgMakeGlobalGrid(grid, 2, 0, 2, tsg_level, tsg_leja, aweights=(/2,1/)) -points => tsgGetPoints(grid) -weights => tsgGetQuadratureWeights(grid) -if ( (norm2d(pointsb-points) > 1.D-11) .OR. (abs(sum(weights)-4.0D+0) > 1.D-11))then - write(*,*) "Mismatch in tsgMakeGlobal: anisotropy", norm2d(pointsb-points), sum(weights) - stop 1 -endif -deallocate(pointsb, points, weights) - - -! make custom rule file -open(unit=17,iostat=i,file='tasmanianFortranTestCustomRule.table',status='replace') -if (i .ne. 0) then - print*, "Cannot open file with custom rule" - stop 1 -endif -write(17,*) "description: test Gauss-Legendre" -write(17,*) "levels: 5" -write(17,*) "1 1 2 3 3 5 4 7 5 9" -do i = 0,4 - call tsgMakeGlobalGrid(grid, 1, 0, i, tsg_level, tsg_gauss_legendre) - points => tsgGetPoints(grid) - weights => tsgGetQuadratureWeights(grid) - do j = 1,tsgGetNumPoints(grid) - write(17,*) weights(j), points(1,j) - enddo - deallocate(points,weights) -enddo -close(unit=17) -! test custom rule -call tsgMakeGlobalGrid(grid, 2, 0, 4, tsg_level, tsg_custom_tabulated, & - customRuleFilename="tasmanianFortranTestCustomRule.table") -call tsgMakeGlobalGrid(grid_II, 2, 0, 4, tsg_level, tsg_gauss_legendre) -points => tsgGetPoints(grid) -pointsb => tsgGetPoints(grid_II) -weights => tsgGetQuadratureWeights(grid) -weightsb => tsgGetQuadratureWeights(grid_II) -if ( (norm1d(weights-weightsb) > 1.d-11) .or. (norm2d(points-pointsb) > 1.d-11) ) then - write(*,*) "Mismatch in tsgMakeGlobal: custom rule" -endif -deallocate(points,pointsb,weights,weightsb) -! delete custom rule file -open(unit=17,file='tasmanianFortranTestCustomRule.table',status='replace') -close(unit=17,status='delete') - - -! test conformal -do i = 7,8 - call tsgMakeGlobalGrid(grid, 2, 0, i, tsg_qptotal, tsg_gauss_patterson) - call tsgMakeGlobalGrid(grid_II, 2, 0, i, tsg_qptotal, tsg_gauss_patterson) - call tsgSetConformalTransformASIN(grid_II, (/4,4/)) - points => tsgGetPoints(grid) - pointsb => tsgGetPoints(grid_II) - weights => tsgGetQuadratureWeights(grid) - weightsb => tsgGetQuadratureWeights(grid_II) - - int = sum( weights * ( 1.d0/((1.d0+5.d0*points(1,:)**2)*(1.d0+5.d0*points(2,:)**2 )) ) ) - int2 = sum( weightsb * ( 1.d0/((1.d0+5.d0*pointsb(1,:)**2)*(1.d0+5.d0*pointsb(2,:)**2)) ) ) - - deallocate(points, pointsb, weights, weightsb) - if ( abs(int-1.028825601981092d0**2) < abs(int2 - 1.028825601981092d0**2) ) then - write(*,*) "Mismatch in tsgMakeGlobal: conformal map" - stop 1 - endif -enddo - - -! test level limits -call tsgMakeGlobalGrid(grid, 2, 0, 20, tsg_qptotal, tsg_clenshaw_curtis, (/1,1/), 0.0D+0, 0.0D+0, levelLimits=(/1,3/)) -call tsgMakeGlobalGrid(grid_II, 2, 0, 1, tsg_tensor, tsg_clenshaw_curtis, (/1,3/)) -points => tsgGetPoints(grid) -pointsb => tsgGetPoints(grid_II) -weights => tsgGetQuadratureWeights(grid) -weightsb => tsgGetQuadratureWeights(grid_II) -if ( (tsgGetNumPoints(grid_II) .ne. tsgGetNumPoints(grid)) .or. & - (norm2d(points-pointsb) > 1.0D-11) .or. (norm1d(weights-weightsb) > 1.0D-11) ) then - write(*,*) "Mismatch in tsgMakeGlobal: level limits", tsgGetNumPoints(grid_II), norm2d(points-pointsb) - stop 1 -end if -deallocate(points, pointsb, weights, weightsb) - - - - - - - - -!======================================================================= -! tsgCopyGrid() -!======================================================================= -call tsgMakeGlobalGrid(grid, 2, 1, 3, tsg_iptotal, tsg_clenshaw_curtis) -call tsgSetConformalTransformASIN(grid, (/4,4/)) -points => tsgGetPoints(grid) -call tsgCopyGrid(grid_II,grid) -call tsgDeallocateGrid(grid) -pointsb => tsgGetPoints(grid_II) -if ( norm2d(points-pointsb) > 1.0D-11 ) then - write(*,*) "Mismatch in tsgMakeGlobal: copy grid" - stop 1 -end if -call tsgAllocateGrid(grid) -deallocate(points, pointsb) - - - - -!======================================================================= -! tsgRead/Write() -!======================================================================= -call tsgMakeFourierGrid(grid, 3, 1, 3, tsg_level, (/1, 2, 3/)) -points => tsgGetPoints(grid) -call tsgWrite(grid, filename) -call tsgDeallocateGrid(grid) -call tsgAllocateGrid(grid) ! reset the grid -if (.not. tsgRead(grid, filename)) then - write(*,*) "File read failed: 1" - stop 1 -end if -pointsb => tsgGetPoints(grid) -if ( norm2d(points-pointsb) > 1.0D-11 ) then - write(*,*) "Mismatch in tsgMakeGlobal: read/write 1" - stop 1 -end if -deallocate(pointsb) - -call tsgWrite(grid, filename, .false.) -call tsgDeallocateGrid(grid) -call tsgAllocateGrid(grid) ! reset the grid -if (.not. tsgRead(grid, filename)) then - write(*,*) "File read failed: 2" - stop 1 -end if -pointsb => tsgGetPoints(grid) -if ( norm2d(points-pointsb) > 1.0D-11 ) then - write(*,*) "Mismatch in tsgMakeGlobal: read/write 2" - stop 1 -end if -deallocate(points, pointsb) - -call tsgWrite(grid, filename, .true.) -call tsgDeallocateGrid(grid) -call tsgAllocateGrid(grid) ! reset the grid -if (.not. tsgRead(grid, filename)) then - write(*,*) "File read failed: 3" - stop 1 -end if -if ( tsgGetRule(grid) /= tsg_fourier ) then - write(*,*) "Mismatch in tsgMakeGlobal: read/write 3" - stop 1 -end if - -if (tsgRead(grid, 'nonafile')) then - write(*,*) "File read failed: 4" - stop 1 -end if - - - - -!======================================================================= -! tsgMakeSequenceGrid() -!======================================================================= -call tsgMakeSequenceGrid(grid, 2, 1, 3, tsg_level, tsg_min_lebesgue) -points => tsgGetPoints(grid) -allocate(pointsb(2,10)) -pointsb = reshape((/ 0.D0, 0.D0, 0.D0, 1.D0, 0.D0, -1.D0, 0.D0, sqrt(1.D0/3.D0), 1.D0, 0.D0, & - 1.D0, 1.D0, 1.D0, -1.D0, -1.D0, 0.D0, -1.D0, 1.D0, sqrt(1.D0/3.D0), 0.D0 /), shape(pointsb)) -if ( norm2d(points-pointsb) > 1.0D-11 ) then - write(*,*) "Mismatch in tsgMakeSequenceGrid: core case 1" - stop 1 -endif -deallocate(points, pointsb) - - -! test transform -call tsgMakeSequenceGrid(grid, 3, 1, 2, tsg_level, tsg_rleja) -call tsgSetDomainTransform(grid, transformA=(/3.0d0,-7.0d0,-12.0d0/), transformB=(/5.0d0,-6.0d0,17.0d0/)) -points => tsgGetPoints(grid) -weights => tsgGetQuadratureWeights(grid) -if ( (abs(sum(weights)-58.d0) > 1.d-11) .or. & - (abs(maxval(points(1,:))-5.d0) > 1.d-11) .or. (abs(minval(points(1,:))-3.d0) > 1.d-11) .or. & - (abs(maxval(points(2,:))+6.d0) > 1.d-11) .or. (abs(minval(points(2,:))+7.d0) > 1.d-11) .or. & - (abs(maxval(points(3,:))-17.d0) > 1.d-11) .or. (abs(minval(points(3,:))+12.d0) > 1.d-11) ) then - write(*,*) "Mismatch in tsgMakeSequenceGrid: transform " - stop 1 -endif - -allocate( double_1d_a(tsgGetNumPoints(grid)) ) -call tsgGetQuadratureWeightsStatic(grid, double_1d_a) -if ( (abs(sum(weights)-58.d0) > 1.d-11) ) then - write(*,*) "Mismatch in tsgGetQuadratureWeightsStatic" - stop 1 -endif -deallocate( double_1d_a ) - -allocate(double_1d_a(3), double_1d_b(3)) -call tsgGetDomainTransform(grid, double_1d_a, double_1d_b) -if ( (norm1d(double_1d_a - (/ 3.0d0, -7.0d0, -12.0d0 /)) > 1.d-11) .or. & - (norm1d(double_1d_b - (/ 5.0d0, -6.0d0, 17.0d0 /)) > 1.d-11) ) then - write(*,*) "Mismatch in tsgGetDomainTransform " - stop 1 -endif -deallocate(points, weights, double_1d_a, double_1d_b) - -call tsgClearDomainTransform(grid) -points => tsgGetPoints(grid) -if ( tsgIsSetDomainTransform(grid) .or. & - (abs(maxval(points(1,:))-1.d0) > 1.d-11) .or. (abs(minval(points(1,:))+1.d0) > 1.d-11) .or. & - (abs(maxval(points(2,:))-1.d0) > 1.d-11) .or. (abs(minval(points(2,:))+1.d0) > 1.d-11) .or. & - (abs(maxval(points(3,:))-1.d0) > 1.d-11) .or. (abs(minval(points(3,:))+1.d0) > 1.d-11) ) then - write(*,*) "Mismatch in tsgMakeSequenceGrid: cleared transform " - stop 1 -endif -deallocate(points) - - -! test anisotropy -allocate(pointsb(2,4)) -pointsb = reshape((/0.0D+0, 0.0D+0, 0.0D+0, 1.0D+0, 0.0D+0, -1.0D+0, 1.0D+0, 0.0D+0/), shape(pointsb)) -call tsgMakeSequenceGrid(grid, 2, 1, 2, tsg_level, tsg_leja, aweights=(/2,1/)) -points => tsgGetPoints(grid) -weights => tsgGetQuadratureWeights(grid) -if ( (norm2d(pointsb-points) > 1.D-11) .or. (abs(sum(weights)-4.0D+0) > 1.D-11) ) then - write(*,*) "Mismatch in tsgMakeSequenceGrid: anisotropy", norm2d(pointsb-points), sum(weights) - stop 1 -endif -deallocate(pointsb, points, weights) - - -! test conformal -allocate( pointsb(2,100), double_2d_a(1,100), double_2d_b(1,100), double_2d_c(1,100) ) -rnd => random(2,100) -pointsb = -1.d0 + 2.d0 * rnd -double_2d_a(1,:) = 1.d0/((1.d0+5.d0*pointsb(1,:)**2)*(1.d0+5.d0*pointsb(2,:)**2)) -do i = 7,8 - call tsgMakeSequenceGrid(grid, 2, 1, i, tsg_iptotal, tsg_rleja) - points => tsgGetPoints(grid) - allocate(double_2d_d(1, tsgGetNumPoints(grid))) - double_2d_d(1,:) = 1.d0/((1.d0+5.d0*points(1,:)**2)*(1.d0+5.d0*points(2,:)**2)) - call tsgLoadNeededPoints(grid, double_2d_d) - call tsgEvaluateBatch(grid,pointsb,100,double_2d_b) - deallocate(points, double_2d_d) - - call tsgMakeSequenceGrid(grid_II, 2, 1, i, tsg_iptotal, tsg_rleja) - call tsgSetConformalTransformASIN(grid_II, (/4,4/)) - points => tsgGetPoints(grid_II) - allocate(double_2d_d(1, tsgGetNumPoints(grid_II))) - double_2d_d(1,:) = 1.d0/((1.d0+5.d0*points(1,:)**2)*(1.d0+5.d0*points(2,:)**2)) - call tsgLoadNeededPoints(grid_II, double_2d_d) - call tsgEvaluateBatch(grid_II,pointsb,100,double_2d_c) - deallocate(points, double_2d_d) - - if ( norm2d(double_2d_a-double_2d_b) < norm2d(double_2d_a-double_2d_c) ) then - write(*,*) "Mismatch in tsgMakeSequenceGrid: conformal map" - stop 1 - endif -enddo -deallocate(rnd, pointsb, double_2d_a, double_2d_b, double_2d_c) - - -! test level limits -call tsgMakeSequenceGrid(grid, 2, 1, 20, tsg_iptotal, tsg_min_delta, levelLimits=(/1,3/)) -call tsgMakeSequenceGrid(grid_II, 2, 1, 1, tsg_tensor, tsg_min_delta, aweights=(/1,3/)) -points => tsgGetPoints(grid) -pointsb => tsgGetPoints(grid_II) -weights => tsgGetQuadratureWeights(grid) -weightsb => tsgGetQuadratureWeights(grid_II) -if ( (tsgGetNumPoints(grid_II) .ne. tsgGetNumPoints(grid)) .or. & - (norm2d(points-pointsb) > 1.0D-11) .or. (norm1d(weights-weightsb) > 1.0D-11) ) then - write(*,*) "Mismatch in tsgMakeSequenceGrid: level limits", tsgGetNumPoints(grid_II), norm2d(points-pointsb) - stop 1 -end if -deallocate(points, pointsb, weights, weightsb) - - - - - - - - -!======================================================================= -! tsgMakeLocalPolynomialGrid() -!======================================================================= -! test transform -call tsgMakeLocalPolynomialGrid(grid, 3, 1, 2, 2, tsg_localp) -call tsgSetDomainTransform(grid, transformA=(/3.0d0,-7.0d0,-12.0d0/), transformB=(/5.0d0,-6.0d0,17.0d0/)) -points => tsgGetPoints(grid) -weights => tsgGetQuadratureWeights(grid) -if ( (abs(sum(weights)-58.d0) > 1.d-11) .or. & - (abs(maxval(points(1,:))-5.d0) > 1.d-11) .or. (abs(minval(points(1,:))-3.d0) > 1.d-11) .or. & - (abs(maxval(points(2,:))+6.d0) > 1.d-11) .or. (abs(minval(points(2,:))+7.d0) > 1.d-11) .or. & - (abs(maxval(points(3,:))-17.d0) > 1.d-11) .or. (abs(minval(points(3,:))+12.d0) > 1.d-11) ) then - write(*,*) "Mismatch in tsgMakeLocalPolynomialGrid: transform " - stop 1 -endif -deallocate(points, weights) - - -! test conformal -allocate( pointsb(2,100), double_2d_a(1,100), double_2d_b(1,100), double_2d_c(1,100) ) -rnd => random(2,100) -pointsb = -1.d0 + 2.d0 * rnd -double_2d_a(1,:) = 1.d0/((1.d0+5.d0*pointsb(1,:)**2)*(1.d0+5.d0*pointsb(2,:)**2)) -do i = 3,4 - call tsgMakeLocalPolynomialGrid(grid, 2, 1, i, 2, tsg_semi_localp) - points => tsgGetPoints(grid) - allocate(double_2d_d(1, tsgGetNumPoints(grid))) - double_2d_d(1,:) = 1.d0/((1.d0+5.d0*points(1,:)**2)*(1.d0+5.d0*points(2,:)**2)) - call tsgLoadNeededPoints(grid, double_2d_d) - call tsgEvaluateBatch(grid,pointsb,100,double_2d_b) - deallocate(points, double_2d_d) - - call tsgMakeLocalPolynomialGrid(grid_II, 2, 1, i, 2, tsg_semi_localp) - call tsgSetConformalTransformASIN(grid_II, (/4,4/)) - points => tsgGetPoints(grid_II) - allocate(double_2d_d(1, tsgGetNumPoints(grid_II))) - double_2d_d(1,:) = 1.d0/((1.d0+5.d0*points(1,:)**2)*(1.d0+5.d0*points(2,:)**2)) - call tsgLoadNeededPoints(grid_II, double_2d_d) - call tsgEvaluateBatch(grid_II,pointsb,100,double_2d_c) - deallocate(points, double_2d_d) - - if ( norm2d(double_2d_a-double_2d_b) < norm2d(double_2d_a-double_2d_c) ) then - write(*,*) "Mismatch in tsgMakeSequenceGrid: conformal map" - stop 1 - endif -enddo -deallocate(rnd, pointsb, double_2d_a, double_2d_b, double_2d_c) - - -! test level limits -call tsgMakeLocalPolynomialGrid(grid, 3, 1, 3, 2, tsg_semi_localp, levelLimits=(/1,2,3/)) -points => tsgGetPoints(grid) -if ( minval(abs(points(1,:)-0.5)) < 1.e-8 ) then - write(*,*) "Mismatch in tsgMakeLocalPolynomialGrid: level limits, dim 1" - stop 1 -endif -if ( minval(abs(points(2,:)-0.75)) < 1.e-8 ) then - write(*,*) "Mismatch in tsgMakeLocalPolynomialGrid: level limits, dim 2" - stop 1 -endif -if ( minval(abs(points(3,:)-0.125)) < 1.e-8 ) then - write(*,*) "Mismatch in tsgMakeLocalPolynomialGrid: level limits, dim 3" - stop 1 -endif -deallocate(points) - - - - - - - - -!======================================================================= -! tsgMakeWaveletGrid() -!======================================================================= -! test transform -call tsgMakeWaveletGrid(grid, 3, 1, 2, 1) -call tsgSetDomainTransform(grid, transformA=(/3.0d0,-7.0d0,-12.0d0/), transformB=(/5.0d0,-6.0d0,17.0d0/)) -points => tsgGetPoints(grid) -weights => tsgGetQuadratureWeights(grid) -if ( (abs(sum(weights)-58.d0) > 1.d-11) .or. & - (abs(maxval(points(1,:))-5.d0) > 1.d-11) .or. (abs(minval(points(1,:))-3.d0) > 1.d-11) .or. & - (abs(maxval(points(2,:))+6.d0) > 1.d-11) .or. (abs(minval(points(2,:))+7.d0) > 1.d-11) .or. & - (abs(maxval(points(3,:))-17.d0) > 1.d-11) .or. (abs(minval(points(3,:))+12.d0) > 1.d-11) ) then - write(*,*) "Mismatch in tsgMakeWaveletGrid: transform " - stop 1 -endif -deallocate(points, weights) - - -! correctness of 1-D -call tsgMakeWaveletGrid(grid, 1, 1, 2, 1) -call tsgMakeLocalPolynomialGrid(grid_II, 1, 1, 3, 1, tsg_localp) -points => tsgGetPoints(grid) -pointsb => tsgGetPoints(grid_II) -if ( norm2d(points-pointsb) > 1.d-11 ) then - write(*,*) "Mismatch in tsgMakeWaveletGrid: points " - stop 1 -endif -deallocate(points, pointsb) - - -! test conformal -do i = 3,4 - call tsgMakeWaveletGrid(grid, 2, 1, i, 1) - call tsgMakeWaveletGrid(grid_II, 2, 1, i, 1) - call tsgSetConformalTransformASIN(grid_II, (/4,4/)) - points => tsgGetPoints(grid) - pointsb => tsgGetPoints(grid_II) - weights => tsgGetQuadratureWeights(grid) - weightsb => tsgGetQuadratureWeights(grid_II) - - int = sum( weights * ( 1.d0/((1.d0+5.d0*points(1,:)**2)*(1.d0+5.d0*points(2,:)**2)) ) ) - int2 = sum( weightsb * ( 1.d0/((1.d0+5.d0*pointsb(1,:)**2)*(1.d0+5.d0*pointsb(2,:)**2)) ) ) - - if ( abs(int-1.028825601981092d0**2) < abs(int2 - 1.028825601981092d0**2) ) then - write(*,*) "Mismatch in tsgMakeWaveletGrid: conformal map" - stop 1 - endif - deallocate(points, pointsb, weights, weightsb) -enddo - - -! test level limits -call tsgMakeWaveletGrid(grid, 3, 1, 2, 1, levelLimits=(/0,1,2/)) -points => tsgGetPoints(grid) -if ( minval(abs(points(1,:)-0.5)) < 1.e-8 ) then - write(*,*) "Mismatch in tsgMakeWaveletGrid: level limits, dim 1" - stop 1 -endif -if ( minval(abs(points(2,:)-0.75)) < 1.e-8 ) then - write(*,*) "Mismatch in tsgMakeWaveletGrid: level limits, dim 2" - stop 1 -endif -if ( minval(abs(points(3,:)-0.125)) < 1.e-8 ) then - write(*,*) "Mismatch in tsgMakeWaveletGrid: level limits, dim 3" - stop 1 -endif -deallocate(points) - - - - - - - - -!======================================================================= -! tsgMakeFourierGrid() -!======================================================================= -! test transform -call tsgMakeFourierGrid(grid, 3, 1, 1, tsg_level) -call tsgSetDomainTransform(grid, transformA=(/3.0d0,-7.0d0,-12.0d0/), transformB=(/5.0d0,-6.0d0,17.0d0/)) -points => tsgGetPoints(grid) -weights => tsgGetQuadratureWeights(grid) -if ( (abs(sum(weights)-58.d0) > 1.d-11) .or. & - (abs(maxval(points(1,:))-13./3.d0) > 1.d-11) .or. (abs(minval(points(1,:))-3.d0) > 1.d-11) .or. & - (abs(maxval(points(2,:))+19./3.d0) > 1.d-11) .or. (abs(minval(points(2,:))+7.d0) > 1.d-11) .or. & - (abs(maxval(points(3,:))-22./3.d0) > 1.d-11) .or. (abs(minval(points(3,:))+12.d0) > 1.d-11) ) then - write(*,*) "Mismatch in tsgMakeFourierGrid: transform " - stop 1 -endif -deallocate(points, weights) - - -! correctness of 1-D -call tsgMakeFourierGrid(grid, 1, 1, 2, tsg_level) -allocate(pointsb(1,9)) -pointsb = reshape( (/ 0.d0, 1.0/3.d0, 2.0/3.d0, 1.0/9.d0, 2.0/9.d0, & - 4.0/9.d0, 5.0/9.d0, 7.0/9.d0, 8.0/9.d0 /), (/1, 9/) ) -points => tsgGetPoints(grid) -if ( norm2d(points-pointsb) > 1.d-11 ) then - write(*,*) "Mismatch in tsgMakeFourierGrid: points " - stop 1 -endif -deallocate(points, pointsb) - - -! test level limits -call tsgMakeFourierGrid(grid, 3, 1, 3, tsg_level, levelLimits=(/0,1,2/)) -points => tsgGetPoints(grid) -if ( maxval(abs(points(1,:))) > 1.e-8 ) then - write(*,*) "Mismatch in tsgMakeFourierGrid: level limits, dim 1" - stop 1 -endif -if ( minval(abs(points(2,:)-1./9.d0)) < 1.e-8 ) then - write(*,*) "Mismatch in tsgMakeFourierGrid: level limits, dim 2" - stop 1 -endif -if ( minval(abs(points(3,:)-1./27.d0)) < 1.e-8 ) then - write(*,*) "Mismatch in tsgMakeFourierGrid: level limits, dim 3" - stop 1 -endif -deallocate(points) - - - -write(*,*) "tsgMake* functions: PASS" - - - - -!======================================================================= -! tsgGet***Points() -!======================================================================= -allocate(pointsb(2,5)) -pointsb = reshape( (/ 0.d0, 0.d0, 0.d0, -1.d0, 0.d0, & - 1.d0, -1.d0, 0.d0, 1.d0, 0.d0 /), (/ 2, 5 /) ) -call tsgMakeGlobalGrid(grid, 2, 1, 1, tsg_level, tsg_clenshaw_curtis) - -points => tsgGetPoints(grid) -if ( norm2d(points-pointsb) > 1.d-11 ) then - write(*,*) "Mismatch in tsgGetPoints: core case 1" - stop 1 -endif -deallocate(points) - -points => tsgGetNeededPoints(grid) -if ( norm2d(points-pointsb) > 1.d-11 ) then - write(*,*) "Mismatch in tsgGetNeededPoints: core case 1" - stop 1 -endif -deallocate(points) - -allocate(points(tsgGetNumDimensions(grid),tsgGetNumPoints(grid))) -call tsgGetPointsStatic(grid, points ) -if ( norm2d(points-pointsb) > 1.d-11 ) then - write(*,*) "Mismatch in tsgGetPointsStatic: core case 1" - stop 1 -endif -deallocate(points) - -allocate(points(tsgGetNumDimensions(grid),tsgGetNumPoints(grid))) -call tsgGetNeededPointsStatic(grid, points ) -if ( norm2d(points-pointsb) > 1.d-11 ) then - write(*,*) "Mismatch in tsgGetNeededPointsStatic: core case 1" - stop 1 -endif -deallocate(points, pointsb) - - - - - - - - -!======================================================================= -! tsgLoadNeededPoints(), tsgGetLoadedPoints() -!======================================================================= -call tsgMakeGlobalGrid(grid, 2, 2, 4, tsg_level, tsg_min_delta) -points => tsgGetPoints(grid) -pointsb => tsgGetNeededPoints(grid) -if ( norm2d(points-pointsb) > 1.d-11 ) then - write(*,*) "Mismatch in tsgLoadNeededPoints: tsgGetNeededPoints" - stop 1 -endif -deallocate(pointsb) -allocate(double_2d_a(2,tsgGetNumPoints(grid))) -double_2d_a(1,:) = exp( -points(1,:)**2 - points(2,:)**2 ) -double_2d_a(2,:) = cos( -points(1,:) - 2.d0 * points(2,:)**2 ) -call tsgLoadNeededPoints(grid, double_2d_a) -pointsb => tsgGetLoadedPoints(grid) -if ( norm2d(points-pointsb) > 1.d-11 ) then - write(*,*) "Mismatch in tsgLoadNeededPoints: tsgGetLoadedPoints" - stop 1 -endif -deallocate(pointsb) -allocate(pointsb(tsgGetNumDimensions(grid),tsgGetNumPoints(grid))) -call tsgGetLoadedPointsStatic(grid, pointsb) -if ( norm2d(points-pointsb) > 1.d-11 ) then - write(*,*) "Mismatch in tsgLoadNeededPoints: tsgGetLoadedPointsStatic" - stop 1 -endif -if ( tsgGetNumNeeded(grid) .ne. 0 ) then - write(*,*) "Mismatch in tsgLoadNeededPoints: tsgGetNumNeeded" - stop 1 -endif -deallocate(points, pointsb, double_2d_a) - - - - - - - - -!======================================================================= -! tsgEvaluate(), tsgEvaluateBatch() -!======================================================================= -allocate(double_2d_a(4,2), double_2d_c(4,2), double_1d_a(4), double_1d_b(4), pointsb(2,2)) -pointsb = reshape( (/ 1./3.d0, 1./3.d0, pi/6.d0, -sqrt(2.d0)/2.d0 /), shape(pointsb) ) -double_2d_a(1,:) = 0.3d0 -double_2d_a(2,:) = pointsb(1,:) + pointsb(2,:) -double_2d_a(3,:) = pointsb(1,:)**2 + pointsb(2,:)**2 + pointsb(1,:)*pointsb(2,:) -double_2d_a(4,:) = pointsb(1,:)**3 + pointsb(2,:)**3 + pointsb(1,:)*(pointsb(2,:)**2) -call tsgMakeLocalPolynomialGrid(grid, 2, 4, 1, 1, tsg_localp) -points => tsgGetPoints(grid) -allocate(double_2d_b(4,tsgGetNumPoints(grid))) -double_2d_b(1,:) = 0.3d0 -double_2d_b(2,:) = points(1,:) + points(2,:) -double_2d_b(3,:) = points(1,:)**2 + points(2,:)**2 + points(1,:)*points(2,:) -double_2d_b(4,:) = points(1,:)**3 + points(2,:)**3 + points(1,:)*(points(2,:)**2) -call tsgLoadNeededPoints(grid, double_2d_b) -call tsgEvaluate(grid,reshape(pointsb(:,1),(/2/)),double_1d_a) -call tsgEvaluateFast(grid,reshape(pointsb(:,2),(/2/)),double_1d_b) -call tsgEvaluateBatch(grid,pointsb,2,double_2d_c) -do i=1,2 - if ( ( norm1d(double_2d_a(i,:)-double_2d_c(i,:)) > 1.d-11 ) .or. & - ( abs(double_2d_a(i,1)-double_1d_a(i)) + abs(double_2d_a(i,2)-double_1d_b(i)) ) > 1.d-11 ) then - write(*,*) "Mismatch in tsgEvaluateBatch: case 1, output ", i - stop 1 - endif -enddo -do i=3,4 - if ( ( norm1d(double_2d_a(i,:)-double_2d_c(i,:)) < 1.d-8 ) .or. & - ( abs(double_2d_a(i,1)-double_1d_a(i)) + abs(double_2d_a(i,2)-double_1d_b(i)) ) < 1.d-8 ) then - write(*,*) "Mismatch in tsgEvaluateBatch: case 1, output ", i - stop 1 - endif -enddo -deallocate(double_2d_a, double_2d_b, double_2d_c, double_1d_a, double_1d_b, pointsb, points) - - -allocate(double_2d_a(4,2), double_2d_c(4,2), double_1d_a(4), double_1d_b(4), pointsb(2,2)) -pointsb = reshape( (/ 1./3.d0, 1./3.d0, pi/6.d0, -sqrt(2.d0)/2.d0 /), shape(pointsb) ) -double_2d_a(1,:) = 0.3d0 -double_2d_a(2,:) = pointsb(1,:) + pointsb(2,:) -double_2d_a(3,:) = pointsb(1,:)**2 + pointsb(2,:)**2 + pointsb(1,:)*pointsb(2,:) -double_2d_a(4,:) = pointsb(1,:)**3 + pointsb(2,:)**3 + pointsb(1,:)*(pointsb(2,:)**2) -call tsgMakeLocalPolynomialGrid(grid, 2, 4, 1, 2, tsg_localp) -points => tsgGetPoints(grid) -allocate(double_2d_b(4,tsgGetNumPoints(grid))) -double_2d_b(1,:) = 0.3d0 -double_2d_b(2,:) = points(1,:) + points(2,:) -double_2d_b(3,:) = points(1,:)**2 + points(2,:)**2 + points(1,:)*points(2,:) -double_2d_b(4,:) = points(1,:)**3 + points(2,:)**3 + points(1,:)*(points(2,:)**2) -call tsgLoadNeededPoints(grid, double_2d_b) -call tsgEvaluate(grid,reshape(pointsb(:,1),(/2/)),double_1d_a) -call tsgEvaluate(grid,reshape(pointsb(:,2),(/2/)),double_1d_b) -call tsgEvaluateBatch(grid,pointsb,2,double_2d_c) -do i=1,2 - if ( norm1d(double_2d_a(i,:)-double_2d_c(i,:)) > 1.d-11 .or. & - ( abs(double_2d_a(i,1)-double_1d_a(i)) + abs(double_2d_a(i,2)-double_1d_b(i)) ) > 1.d-11 ) then - write(*,*) "Mismatch in tsgEvaluateBatch: case 2, output ", i - stop 1 - endif -enddo -do i=3,4 - if ( ( norm1d(double_2d_a(i,:)-double_2d_c(i,:)) < 1.d-8 ) .or. & - ( abs(double_2d_a(i,1)-double_1d_a(i)) + abs(double_2d_a(i,2)-double_1d_b(i)) ) < 1.d-8 ) then - write(*,*) "Mismatch in tsgEvaluateBatch: case 2, output ", i - stop 1 - endif -enddo -deallocate(double_2d_a, double_2d_b, double_2d_c, double_1d_a, double_1d_b, pointsb, points) - - -allocate(double_2d_a(4,2), double_2d_c(4,2), double_1d_a(4), double_1d_b(4), pointsb(2,2)) -pointsb = reshape( (/ 1./3.d0, 1./3.d0, pi/6.d0, -sqrt(2.d0)/2.d0 /), shape(pointsb) ) -double_2d_a(1,:) = 0.3d0 -double_2d_a(2,:) = pointsb(1,:) + pointsb(2,:) -double_2d_a(3,:) = pointsb(1,:)**2 + pointsb(2,:)**2 -double_2d_a(4,:) = pointsb(1,:)**3 + pointsb(2,:)**3 + pointsb(1,:)*(pointsb(2,:)**2) -call tsgMakeLocalPolynomialGrid(grid, 2, 4, 1, 2, tsg_semi_localp) -points => tsgGetPoints(grid) -allocate(double_2d_b(4,tsgGetNumPoints(grid))) -double_2d_b(1,:) = 0.3d0 -double_2d_b(2,:) = points(1,:) + points(2,:) -double_2d_b(3,:) = points(1,:)**2 + points(2,:)**2 -double_2d_b(4,:) = points(1,:)**3 + points(2,:)**3 + points(1,:)*(points(2,:)**2) -call tsgLoadNeededPoints(grid, double_2d_b) -call tsgEvaluate(grid,reshape(pointsb(:,1),(/2/)),double_1d_a) -call tsgEvaluate(grid,reshape(pointsb(:,2),(/2/)),double_1d_b) -call tsgEvaluateBatch(grid,pointsb,2,double_2d_c) -do i=1,3 - if ( norm1d(double_2d_a(i,:)-double_2d_c(i,:)) > 1.d-11 .or. & - ( abs(double_2d_a(i,1)-double_1d_a(i)) + abs(double_2d_a(i,2)-double_1d_b(i)) ) > 1.d-11 ) then - write(*,*) "Mismatch in tsgEvaluateBatch: case 3, output ", i - stop 1 - endif -enddo -do i=4,4 - if ( ( norm1d(double_2d_a(i,:)-double_2d_c(i,:)) < 1.d-8 ) .or. & - ( abs(double_2d_a(i,1)-double_1d_a(i)) + abs(double_2d_a(i,2)-double_1d_b(i)) ) < 1.d-8 ) then - write(*,*) "Mismatch in tsgEvaluateBatch: case 3, output ", i - stop 1 - endif -enddo -deallocate(double_2d_a, double_2d_b, double_2d_c, double_1d_a, double_1d_b, pointsb, points) - - -allocate(double_2d_a(4,2), double_2d_c(4,2), double_1d_a(4), double_1d_b(4), pointsb(2,2)) -pointsb = reshape( (/ 1./3.d0, 1./3.d0, pi/6.d0, -sqrt(2.d0)/2.d0 /), shape(pointsb) ) -double_2d_a(1,:) = 0.3d0 -double_2d_a(2,:) = pointsb(1,:) + pointsb(2,:) -double_2d_a(3,:) = pointsb(1,:)**2 + pointsb(2,:)**2 + pointsb(1,:)*pointsb(2,:) -double_2d_a(4,:) = pointsb(1,:)**3 + pointsb(2,:)**3 + pointsb(1,:)*(pointsb(2,:)**2) -call tsgMakeLocalPolynomialGrid(grid, 2, 4, 2, 2, tsg_localp) -points => tsgGetPoints(grid) -allocate(double_2d_b(4,tsgGetNumPoints(grid))) -double_2d_b(1,:) = 0.3d0 -double_2d_b(2,:) = points(1,:) + points(2,:) -double_2d_b(3,:) = points(1,:)**2 + points(2,:)**2 + points(1,:)*points(2,:) -double_2d_b(4,:) = points(1,:)**3 + points(2,:)**3 + points(1,:)*(points(2,:)**2) -call tsgLoadNeededPoints(grid, double_2d_b) -call tsgEvaluate(grid,reshape(pointsb(:,1),(/2/)),double_1d_a) -call tsgEvaluate(grid,reshape(pointsb(:,2),(/2/)),double_1d_b) -call tsgEvaluateBatch(grid,pointsb,2,double_2d_c) -do i=1,3 - if ( norm1d(double_2d_a(i,:)-double_2d_c(i,:)) > 1.d-11 .or. & - ( abs(double_2d_a(i,1)-double_1d_a(i)) + abs(double_2d_a(i,2)-double_1d_b(i)) ) > 1.d-11 ) then - write(*,*) "Mismatch in tsgEvaluateBatch: case 4, output ", i - stop 1 - endif -enddo -do i=4,4 - if ( ( norm1d(double_2d_a(i,:)-double_2d_c(i,:)) < 1.d-8 ) .or. & - ( abs(double_2d_a(i,1)-double_1d_a(i)) + abs(double_2d_a(i,2)-double_1d_b(i)) ) < 1.d-8 ) then - write(*,*) "Mismatch in tsgEvaluateBatch: case 4, output ", i - stop 1 - endif -enddo -deallocate(double_2d_a, double_2d_b, double_2d_c, double_1d_a, double_1d_b, pointsb, points) - - -allocate(double_2d_a(4,2), double_2d_c(4,2), double_1d_a(4), double_1d_b(4), pointsb(2,2)) -pointsb = reshape( (/ 1./3.d0, 1./3.d0, pi/6.d0, -sqrt(2.d0)/2.d0 /), shape(pointsb) ) -double_2d_a(1,:) = 0.3d0 -double_2d_a(2,:) = pointsb(1,:) + pointsb(2,:) -double_2d_a(3,:) = pointsb(1,:)**2 + pointsb(2,:)**2 + pointsb(1,:)*pointsb(2,:) -double_2d_a(4,:) = pointsb(1,:)**3 + pointsb(2,:)**3 + pointsb(1,:)*(pointsb(2,:)**2) -call tsgMakeLocalPolynomialGrid(grid, 2, 4, 3, 3, tsg_localp) -points => tsgGetPoints(grid) -allocate(double_2d_b(4,tsgGetNumPoints(grid))) -double_2d_b(1,:) = 0.3d0 -double_2d_b(2,:) = points(1,:) + points(2,:) -double_2d_b(3,:) = points(1,:)**2 + points(2,:)**2 + points(1,:)*points(2,:) -double_2d_b(4,:) = points(1,:)**3 + points(2,:)**3 + points(1,:)*(points(2,:)**2) -call tsgLoadNeededPoints(grid, double_2d_b) -call tsgEvaluate(grid,reshape(pointsb(:,1),(/2/)),double_1d_a) -call tsgEvaluate(grid,reshape(pointsb(:,2),(/2/)),double_1d_b) -call tsgEvaluateBatch(grid,pointsb,2,double_2d_c) -do i=1,4 - if ( ( norm1d(double_2d_a(i,:)-double_2d_c(i,:)) > 1.d-11 ) .or. & - ( abs(double_2d_a(i,1)-double_1d_a(i)) + abs(double_2d_a(i,2)-double_1d_b(i)) ) > 1.d-11 ) then - write(*,*) "Mismatch in tsgEvaluateBatch: case 5, output ", i - stop 1 - endif -enddo -deallocate(double_2d_a, double_2d_b, double_2d_c, double_1d_a, double_1d_b, pointsb, points) - - -call tsgMakeGlobalGrid(grid, 2, 1, 22, tsg_iptotal, tsg_chebyshev) -points => tsgGetPoints(grid) -allocate(double_2d_a(1,tsgGetNumPoints(grid))) -double_2d_a(1,:) = exp( -points(1,:)**2 - points(2,:)**2 ) -call tsgLoadNeededPoints(grid, double_2d_a) -allocate(pointsb(2,1000),double_2d_b(1,1000),double_2d_c(1,1000)) -rnd => random(2,1000) -pointsb = -1.d0 + 2.d0 * rnd -double_2d_b(1,:) = exp( -pointsb(1,:)**2 - pointsb(2,:)**2 ) -call tsgEvaluateBatch(grid,pointsb,1000,double_2d_c) -if ( norm2d(double_2d_b-double_2d_c) > 1.d-8 ) then - write(*,*) "Mismatch in tsgEvaluateBatch: global grid with chebyshev points, output ", norm2d(double_2d_b-double_2d_c) - stop 1 -endif -deallocate(points,pointsb,double_2d_a,double_2d_b,double_2d_c,rnd) - - - - - - - - -!======================================================================= -! tsgEvaluateHierarchicalFunctions() -!======================================================================= -call tsgMakeGlobalGrid(grid, 3, 1, 4, tsg_level, tsg_fejer2) -points => tsgGetPoints(grid); i_a = tsgGetNumPoints(grid) -allocate(double_2d_a(i_a,i_a)) -call tsgEvaluateHierarchicalFunctions(grid,points,i_a,double_2d_a) -forall(i = 1:i_a) double_2d_a(i,i) = double_2d_a(i,i) - 1 -if ( norm2d(double_2d_a) > 1.d-11 ) then - write(*,*) "Mismatch in tsgEvaluateHierarchy: lagrange polynomials do not form identity" - stop 1 -endif -deallocate(points, double_2d_a) - - -call tsgMakeSequenceGrid(grid, 2, 1, 2, tsg_level, tsg_leja) -allocate(points(2,6), double_2d_a(6,6), double_2d_b(6,6)) -points = reshape( (/ 0.33d0, 0.25d0, -0.270d0, 0.39d0, 0.970d0, -0.760d0, & - -0.44d0, 0.21d0, -0.813d0, 0.03d0, -0.666d0, 0.666d0 /), shape(points) ) -double_2d_a(1,:) = 1 -double_2d_a(2,:) = points(2,:) -double_2d_a(3,:) = 0.5d0 * points(2,:) * ( points(2,:) - 1.d0 ) -double_2d_a(4,:) = points(1,:) -double_2d_a(5,:) = points(1,:) * points(2,:) -double_2d_a(6,:) = 0.5d0 * points(1,:) * ( points(1,:) - 1.d0 ) -call tsgEvaluateHierarchicalFunctions(grid,points,6,double_2d_b) -if ( norm2d(double_2d_a-double_2d_b) > 1.d-11 ) then - write(*,*) "Mismatch in tsgEvaluateHierarchicalFunctions: sequence grid test" - stop 1 -endif -deallocate(points, double_2d_a, double_2d_b) - - -call tsgMakeLocalPolynomialGrid(grid, 2, 1, 4, 1, tsg_localp) -points => tsgGetPoints(grid) -i_a = tsgGetNumPoints(grid) -i_b = 13 -allocate(double_2d_a(1,i_a), double_2d_b(1,i_b), double_2d_c(i_a,i_b), double_1d_b(i_b), pointsb(2,i_b)) -double_2d_a(1,:) = exp( -points(1,:)**2 - 2.d0 * points(2,:)**2 ) -call tsgLoadNeededPoints(grid, double_2d_a) -rnd => random(2,i_b) -pointsb = -1.d0 + 2.d0 * rnd -call tsgEvaluateBatch(grid,pointsb,i_b,double_2d_b) -call tsgEvaluateHierarchicalFunctions(grid,pointsb,i_b,double_2d_c) -weights => tsgGetHierarchicalCoefficients(grid) -double_1d_b = matmul( weights, double_2d_c ) -if ( norm1d(double_1d_b-double_2d_b(1,:)) > 1.d-11 ) then - write(*,*) "Mismatch in tsgEvaluateHierarchicalFunctions: local grid test" - stop 1 -endif -allocate( double_1d_a(tsgGetNumPoints(grid)) ) -call tsgGetHierarchicalCoefficientsStatic(grid, double_1d_a) -double_1d_b = matmul( double_1d_a, double_2d_c ) -if ( norm1d(double_1d_b-double_2d_b(1,:)) > 1.d-11 ) then - write(*,*) "Mismatch in tsgGetHierarchicalCoefficientsStatic: local grid test" - stop 1 -endif -deallocate(points,double_1d_a,double_2d_a,double_2d_b,double_2d_c,weights,double_1d_b,rnd,pointsb) - - -call tsgMakeLocalPolynomialGrid(grid, 2, 1, 4, 1, tsg_localp) -points => tsgGetPoints(grid) -i_a = tsgGetNumPoints(grid) -i_b = 13 -allocate(double_2d_a(1,i_a), double_2d_b(1,i_b), double_1d_b(i_b), pointsb(2,i_b)) -double_2d_a(1,:) = exp( -points(1,:)**2 - 2.d0 * points(2,:)**2 ) -call tsgLoadNeededPoints(grid, double_2d_a) -rnd => random(2,i_b) -pointsb = -1.d0 + 2.d0 * rnd -call tsgEvaluateBatch(grid,pointsb,i_b,double_2d_b) -call tsgEvaluateSparseHierarchicalFunctions(grid,pointsb, i_b, int_pnt_1d_a, int_pnt_1d_b, double_pnt_1d_a) -weights => tsgGetHierarchicalCoefficients(grid) -double_1d_b = sparse_matmul( int_pnt_1d_a, int_pnt_1d_b, double_pnt_1d_a, weights ) -if ( norm1d(double_1d_b-double_2d_b(1,:)) > 1.d-11 ) then - write(*,*) "Mismatch in tsgEvaluateSparseHierarchicalFunctions: local grid test" - stop 1 -endif -deallocate(points,double_2d_a,double_2d_b,weights,double_1d_b,rnd,pointsb,int_pnt_1d_a,int_pnt_1d_b,double_pnt_1d_a) - - -! test reading a complex matrix -call tsgMakeFourierGrid(grid, 2, 1, 4, tsg_level) -points => tsgGetPoints(grid) -i_a = tsgGetNumPoints(grid) -i_b = 13 -allocate(double_2d_a(1,i_a), double_2d_b(1,i_b), double_2d_c(i_b,2*i_a), & - double_1d_a(i_b), double_1d_b(i_b), double_1d_c(i_b), pointsb(2,i_b), double_1d_d(2*i_a)) -double_2d_a(1,:) = exp( -points(1,:)**2 - 2.d0 * points(2,:)**2 ) -call tsgLoadNeededPoints(grid, double_2d_a) -rnd => random(2,i_b) -pointsb = -1.d0 + 2.d0 * rnd -call tsgEvaluateBatch(grid,pointsb,i_b,double_2d_b) - -allocate(dcmplx_1d_a(i_a),dcmplx_2d_c(i_a,i_b)) -call tsgEvaluateComplexHierarchicalFunctions(grid,pointsb,i_b,dcmplx_2d_c) -call tsgGetComplexHierarchicalCoefficientsStatic(grid, dcmplx_1d_a) -double_1d_b = real(matmul(dcmplx_1d_a,dcmplx_2d_c)) - -dcmplx_pnt_1d_a => tsgGetComplexHierarchicalCoefficients(grid) -double_1d_c = real(matmul(dcmplx_pnt_1d_a,dcmplx_2d_c)) - -if ( norm1d(double_1d_b - double_2d_b(1,:)) > 1.d-11 ) then - write(*,*) "Mismatch in tsgGetComplexHierarchicalCoefficientsStatic: fourier grid test 1" - stop 1 -endif -if ( norm1d(double_1d_c - double_2d_b(1,:)) > 1.d-11 ) then - write(*,*) "Mismatch in tsgGetComplexHierarchicalCoefficientsStatic: fourier grid test 2" - stop 1 -endif -deallocate(points,pointsb,double_2d_a,double_2d_b,double_2d_c,double_1d_a,double_1d_b,double_1d_c,double_1d_d,rnd) -deallocate(dcmplx_pnt_1d_a,dcmplx_1d_a,dcmplx_2d_c) - -! load coefficients -i_a = 13 -allocate( double_2d_b(1,i_a), double_2d_c(1,i_a), pointsb(2,i_a) ) -call tsgMakeLocalPolynomialGrid(grid, 2, 1, 5, 1, tsg_semi_localp) -call tsgMakeLocalPolynomialGrid(grid_II, 2, 1, 5, 1, tsg_semi_localp) -points => tsgGetPoints(grid); i_b = tsgGetNumPoints(grid) -rnd => random(2,i_a) -pointsb = -1.d0 + 2.d0 * rnd -allocate( double_2d_a(1,i_b), double_2d_d(i_b,i_b) ) -double_2d_a(1,:) = exp( -points(1,:)**2 - 2.d0 * points(2,:)**2 ) -call tsgLoadNeededPoints(grid, double_2d_a) -double_pnt_1d_b => tsgGetHierarchicalCoefficients(grid) -call tsgSetHierarchicalCoefficients(grid_II,double_pnt_1d_b) -call tsgEvaluateBatch(grid,pointsb,i_a,double_2d_b) -call tsgEvaluateBatch(grid_II,pointsb,i_a,double_2d_c) -if ( norm2d(double_2d_b - double_2d_c) > 1.d-11 ) then - write(*,*) "Mismatch in setHierarchicalCoefficients: localp grid solve ", norm2d(double_2d_b - double_2d_c) - stop 1 -endif -deallocate(points,pointsb,double_2d_a,double_2d_b,double_2d_c,double_2d_d,double_pnt_1d_b,rnd) - - -call tsgMakeLocalPolynomialGrid(grid, 2, 1, 5, 1, tsg_semi_localp) -call tsgMakeLocalPolynomialGrid(grid_II, 2, 1, 5, 1, tsg_semi_localp) -call tsgSetDomainTransform(grid, (/-1.d0, 7.d0/), (/2.d0, 9.d0/)) -call tsgSetDomainTransform(grid_II, (/-1.d0, 7.d0/), (/2.d0, 9.d0/)) -points => tsgGetPoints(grid); i_b = tsgGetNumPoints(grid) -i_a = 13 -allocate(double_2d_a(1,i_b), double_2d_b(1,i_a), double_2d_c(1,i_a), double_2d_d(i_b,i_b), pointsb(2,i_a)) -double_2d_a(1,:) = exp( -points(1,:)**2 - 2.d0 * points(2,:)**2 ) -call tsgLoadNeededPoints(grid, double_2d_a) -double_pnt_1d_b => tsgGetHierarchicalCoefficients(grid) -call tsgSetHierarchicalCoefficients(grid_II,double_pnt_1d_b) -rnd => random(2,i_a) -pointsb = -1.d0 + 2.d0 * rnd -call tsgEvaluateBatch(grid, pointsb, i_a, double_2d_b) -call tsgEvaluateBatch(grid_II, pointsb, i_a, double_2d_c) -if ( norm2d(double_2d_b - double_2d_c) > 1.d-11 ) then - write(*,*) "Mismatch in setHierarchicalCoefficients: local grid solve, transform" - stop 1 -endif -deallocate(points,pointsb,double_2d_a,double_2d_b,double_2d_c,double_2d_d,double_pnt_1d_b,rnd) - - -call tsgMakeSequenceGrid(grid, 2, 1, 5, tsg_level, tsg_rleja) -call tsgMakeSequenceGrid(grid_II, 2, 1, 5, tsg_level, tsg_rleja) -call tsgSetDomainTransform(grid, (/1.d0, 1.d0/), (/2.d0, 2.d0/)) -call tsgSetDomainTransform(grid_II, (/1.d0, 1.d0/), (/2.d0, 2.d0/)) -points => tsgGetPoints(grid); i_b = tsgGetNumPoints(grid) -i_a = 13 -allocate(double_2d_a(1,i_b), double_2d_b(1,i_a), double_2d_c(1,i_a), double_2d_d(i_b,i_b), pointsb(2,i_a)) -double_2d_a(1,:) = exp( -points(1,:)**2 - 2.d0 * points(2,:)**2 ) -call tsgLoadNeededPoints(grid, double_2d_a) -double_pnt_1d_b => tsgGetHierarchicalCoefficients(grid) -call tsgSetHierarchicalCoefficients(grid_II,double_pnt_1d_b) -rnd => random(2,i_a) -pointsb = -1.d0 + 2.d0 * rnd -call tsgEvaluateBatch(grid, pointsb, i_a, double_2d_b) -call tsgEvaluateBatch(grid_II, pointsb, i_a, double_2d_c) -if ( norm2d(double_2d_b - double_2d_c) > 1.d-11 ) then - write(*,*) "Mismatch in setHierarchicalCoefficients: sequence grid solve, transform" - stop 1 -endif -deallocate(points,pointsb,double_2d_a,double_2d_b,double_2d_c,double_2d_d,double_pnt_1d_b,rnd) - - -call tsgMakeWaveletGrid(grid, 1, 1, 1, 3) -points => tsgGetHierarchicalSupport(grid) -double_2d_b = reshape( (/ 2.d0, 2.d0, 2.d0, 2.d0, 2.d0, 2.d0, 2.d0, 2.d0, 2.d0, & - 1.d4, 1.d4, 1.d4, 1.d4, 1.d4, 1.d4, 1.d4, 1.d4 /), shape(points) ) -if ( norm2d(points - double_2d_b) > 1.d-11 ) then - write(*,*) "Mismatch in tsgGetHierarchicalSupport: wavelet grid" - stop 1 -endif -deallocate(points, double_2d_b) - - - - - -!======================================================================= -! tsgIntegrate() -!======================================================================= -call tsgMakeGlobalGrid(grid, 1, 1, 2, tsg_level, tsg_gauss_hermite, alpha=0.d0, beta=0.d0) -points => tsgGetPoints(grid) -i_a = tsgGetNumPoints(grid) -allocate(double_2d_a(1,i_a),double_1d_a(tsgGetNumOutputs(grid))) -double_2d_a(1,:) = points(1,:)**2 -call tsgLoadNeededPoints(grid, double_2d_a) -call tsgIntegrate(grid, double_1d_a) -if ( abs(double_1d_a(1) - sqrt(pi)/2.d0) > 1.d-11 ) then - write(*,*) "Mismatch in tsgIntegrate: case 1" - stop 1 -endif -deallocate(points, double_2d_a, double_1d_a) - - -call tsgMakeGlobalGrid(grid, 1, 1, 2, tsg_level, tsg_gauss_hermite, alpha=2.d0, beta=0.d0) -points => tsgGetPoints(grid) -i_a = tsgGetNumPoints(grid) -allocate(double_2d_a(1,i_a),double_1d_a(tsgGetNumOutputs(grid))) -double_2d_a(1,:) = sqrt(2.d0) -call tsgLoadNeededPoints(grid, double_2d_a) -call tsgIntegrate(grid, double_1d_a) -if ( abs(double_1d_a(1) - sqrt(2.d0*pi)/2.d0) > 1.d-11 ) then - write(*,*) "Mismatch in tsgIntegrate: case 2" - stop 1 -endif -deallocate(points, double_2d_a, double_1d_a) - - - - -write(*,*) "Core I/O and evaluate: PASS" - - - -!======================================================================= -! tsgGetInterpolationWeights() -!======================================================================= -i_a = 32 -call tsgMakeGlobalGrid(grid, 2, 1, 4, tsg_level, tsg_fejer2) -call tsgMakeGlobalGrid(grid_II, 2, 1, 4, tsg_level, tsg_fejer2) -points => tsgGetPoints(grid); i_b = tsgGetNumPoints(grid) -allocate( pointsb(2,i_a), double_2d_a(1,i_b), double_2d_b(1,i_a), double_2d_c(1,i_a) ) -rnd => random(2,i_a) -pointsb = -1.d0 + 2.d0 * rnd -double_2d_a(1,:) = exp( -points(1,:)**2 - 2.d0 * points(2,:)**2 ) -call tsgLoadNeededPoints(grid, double_2d_a) -call tsgEvaluateBatch(grid, pointsb, i_a, double_2d_b) -do i = 1,i_a - double_pnt_1d_a => tsgGetInterpolationWeights(grid_II, pointsb(:,i)) - double_2d_c(:,i) = matmul(double_2d_a,double_pnt_1d_a) - deallocate(double_pnt_1d_a) -enddo -if ( norm2d(double_2d_b - double_2d_c) > 1.d-11 ) then - write(*,*) "Mismatch in tsgGetInterpolationWeights: global case" - stop 1 -endif -do i = 1,i_a - allocate( double_1d_a(tsgGetNumPoints(grid)) ) - call tsgGetInterpolationWeightsStatic(grid, pointsb(:,i), double_1d_a) - double_2d_c(:,i) = matmul(double_2d_a,double_1d_a) - deallocate(double_1d_a) -enddo -if ( norm2d(double_2d_b - double_2d_c) > 1.d-11 ) then - write(*,*) "Mismatch in tsgGetInterpolationWeights: global case" - stop 1 -endif -deallocate(points, pointsb, double_2d_a, double_2d_b, double_2d_c, rnd) - - -i_a = 32 -call tsgMakeSequenceGrid(grid, 2, 1, 7, tsg_level, tsg_min_delta) -call tsgMakeSequenceGrid(grid_II, 2, 1, 7, tsg_level, tsg_min_delta) -points => tsgGetPoints(grid); i_b = tsgGetNumPoints(grid) -allocate( pointsb(2,i_a), double_2d_a(1,i_b), double_2d_b(1,i_a), double_2d_c(1,i_a) ) -rnd => random(2,i_a) -pointsb = -1.d0 + 2.d0 * rnd -double_2d_a(1,:) = exp( -points(1,:)**2 - 2.d0 * points(2,:)**2 ) -call tsgLoadNeededPoints(grid, double_2d_a) -call tsgEvaluateBatch(grid, pointsb, i_a, double_2d_b) -do i = 1,i_a - double_pnt_1d_a => tsgGetInterpolationWeights(grid_II, pointsb(:,i)) - double_2d_c(:,i) = matmul(double_2d_a,double_pnt_1d_a) - deallocate(double_pnt_1d_a) -enddo -if ( norm2d(double_2d_b - double_2d_c) > 1.d-11 ) then - write(*,*) "Mismatch in tsgGetInterpolationWeights: global case" - stop 1 -endif -deallocate(points, pointsb, double_2d_a, double_2d_b, double_2d_c, rnd) - - - - - - - - -!======================================================================= -! tsgEstimateAnisotropicCoefficients() -!======================================================================= -call tsgMakeGlobalGrid(grid, 2, 1, 9, tsg_level, tsg_rleja) -points => tsgGetPoints(grid); i_b = tsgGetNumPoints(grid) -allocate( double_2d_a(1,i_b) ) -double_2d_a(1,:) = exp( points(1,:) + points(2,:)**2 ) -call tsgLoadNeededPoints(grid, double_2d_a) -int_pnt_1d_a => tsgEstimateAnisotropicCoefficients(grid,tsg_iptotal,1) -if ( abs( int_pnt_1d_a(1)/dble(int_pnt_1d_a(2)) - 2.0 ) > 0.2 ) then - write(*,*) "Mismatch in tsgEstimateAnisotropicCoefficients: total degree" - stop 1 -endif -deallocate(int_pnt_1d_a) -int_pnt_1d_a => tsgEstimateAnisotropicCoefficients(grid,tsg_ipcurved,1) -if ( size(int_pnt_1d_a) .ne. 4 ) then - write(*,*) "Mismatch in tsgEstimateAnisotropicCoefficients: curved dimensions" - stop 1 -endif -if ( (abs( int_pnt_1d_a(1)/dble(int_pnt_1d_a(2)) - 2.0 ) > 0.2) & - .or. (int_pnt_1d_a(3)>0) .or. (int_pnt_1d_a(4)>0) ) then - write(*,*) "Mismatch in tsgEstimateAnisotropicCoefficients: curved" - stop 1 -endif -deallocate(int_pnt_1d_a) -deallocate(points, double_2d_a) - - - - - - - - -!======================================================================= -! tsgSetAnisotropicRefinement() -!======================================================================= -call tsgMakeSequenceGrid(grid, 3, 1, 3, tsg_level, tsg_leja, levelLimits=(/3,2,1/)) -points => tsgGetPoints(grid); i_b = tsgGetNumPoints(grid) -if ( check_points( points(1,:), (/0.d0,-1.d0,1.d0,1.d0/sqrt(3.d0)/), (/12345.d0/) ) .or. & - check_points( points(2,:), (/0.d0,-1.d0,1.d0/), (/1.d0/sqrt(3.d0)/) ) .or. & - check_points( points(3,:), (/0.d0,1.d0/), (/-1.d0,1.d0/sqrt(3.d0)/) ) ) then - write(*,*) "Mismatch in tsgSetAnisotropicRefinement: limits in make" - stop 1 -endif -allocate( double_2d_a(1,i_b) ) -double_2d_a(1,:) = exp( -points(1,:)**2 - points(2,:)**2 ) -call tsgLoadNeededPoints(grid, double_2d_a) -call tsgSetAnisotropicRefinement( grid, tsg_iptotal, 5, 1 ) -pointsb => tsgGetNeededPoints(grid) -if ( size(pointsb,2) .eq. 0 ) then - write(*,*) "Mismatch in tsgSetAnisotropicRefinement: did not refine" - stop 1 -endif -if ( check_points( points(2,:), (/12345.d0/), (/1.d0/sqrt(3.d0)/) ) .or. & - check_points( points(3,:), (/12345.d0/), (/-1.d0,1.d0/sqrt(3.d0)/) ) ) then - write(*,*) "Mismatch in tsgSetAnisotropicRefinement: limits refine using existing limits" - stop 1 -endif -deallocate( pointsb ) -call tsgSetAnisotropicRefinement(grid,tsg_iptotal,10,1,levelLimits=(/3,2,2/)) -pointsb => tsgGetNeededPoints(grid) -if ( size(pointsb) .eq. 0 ) then - write(*,*) "Mismatch in tsgSetAnisotropicRefinement: did not refine" - stop 1 -endif -if ( check_points( points(2,:), (/12345.d0/), (/1.d0/sqrt(3.d0)/) ) .or. & - check_points( points(3,:), (/-1.d0/), (/1.d0/sqrt(3.d0)/) ) ) then - write(*,*) "Mismatch in tsgSetAnisotropicRefinement: limits refine using new limits" - stop 1 -endif -deallocate( points, pointsb, double_2d_a ) - - - - - - - - -!======================================================================= -! tsgSetLocalSurplusRefinement() -!======================================================================= -call tsgMakeLocalPolynomialGrid(grid, 3, 1, 3, 1, tsg_localp, levelLimits=(/1,2,3/)) -points => tsgGetPoints(grid); i_b = tsgGetNumPoints(grid) -if ( check_points( points(1,:), (/0.d0,-1.d0,1.d0/), (/0.5d0,-0.5d0,-0.75d0,-0.25d0,0.25d0,0.75d0/) ) .or. & - check_points( points(2,:), (/0.d0,-1.d0,1.d0,0.5d0,-0.5d0/), (/-0.75d0,-0.25d0,0.25d0,0.75d0/) ) .or. & - check_points( points(3,:), (/0.d0,-1.d0,1.d0,0.5d0,-0.5d0,-0.75d0,-0.25d0,0.25d0,0.75d0/), (/12345.d0/) ) ) then - write(*,*) "Mismatch in tsgSetLocalSurplusRefinement: limits in make" - stop 1 -endif -allocate( double_2d_a(1,i_b) ) -double_2d_a(1,:) = exp( -points(1,:)**2 - points(2,:)**2 ) -call tsgLoadNeededPoints(grid, double_2d_a) -call tsgSetLocalSurplusRefinement(grid,1.d-8,tsg_classic,1) -if ( tsgGetNumNeeded(grid) .eq. 0 ) then - write(*,*) "Mismatch in tsgSetLocalSurplusRefinement: did not refine local polynomial" - stop 1 -endif -if ( check_points( points(1,:), (/12345.d0/), (/0.5d0,-0.5d0,-0.75d0,-0.25d0,0.25d0,0.75d0/) ) .or. & - check_points( points(2,:), (/12345.d0/), (/-0.75d0,-0.25d0,0.25d0,0.75d0/) ) .or. & - check_points( points(3,:), (/12345.d0/), (/12345.d0/) ) ) then - write(*,*) "Mismatch in tsgSetSurplusRefinement: limits refine using existing limits" - stop 1 -endif -call tsgSetLocalSurplusRefinement(grid,1.d-8,tsg_classic,1,(/2,2,3/)) -if ( tsgGetNumNeeded(grid) .eq. 0 ) then - write(*,*) "Mismatch in tsgSetLocalSurplusRefinement: did not refine on second pass" - stop 1 -endif -if ( check_points( points(1,:), (/0.5d0,-0.5d0/), (/-0.75d0,-0.25d0,0.25d0,0.75d0/) ) .or. & - check_points( points(2,:), (/12345.d0/), (/-0.75d0,-0.25d0,0.25d0,0.75d0/) ) ) then - write(*,*) "Mismatch in tsgSetSurplusRefinement: limits refine using new limits" - stop 1 -endif -deallocate(points, double_2d_a) - - - - - - - - -!======================================================================= -! tsgClearRefinement() -!======================================================================= -call tsgMakeLocalPolynomialGrid(grid, 3, 1, 4, 2, tsg_localp, levelLimits=(/300,300,300/)) -points => tsgGetPoints(grid); i_b = tsgGetNumPoints(grid) -allocate( double_2d_a(1,i_b) ) -double_2d_a(1,:) = exp( -points(1,:)**2 - 0.5*points(2,:)**2 - 2.d0*points(3,:)**2 ) -call tsgLoadNeededPoints(grid, double_2d_a) -if ( tsgGetNumNeeded(grid) > 0 ) then - write(*,*) "Mismatch in cancel refine: did not load values" - stop 1 -endif -call tsgSetLocalSurplusRefinement(grid,1.d-4,tsg_directional) -if ( tsgGetNumNeeded(grid) .eq. 0 ) then - write(*,*) "Mismatch in cancel refine: did not set refinement at output -1" - stop 1 -endif -call tsgClearRefinement(grid) -if ( tsgGetNumNeeded(grid) > 0 ) then - write(*,*) "Mismatch in cancel refine: did not cancel the refinement" - stop 1 -endif -deallocate(points, double_2d_a) - - - - - - - - -!======================================================================= -! tsgMergeRefinement() -!======================================================================= -call tsgMakeGlobalGrid(grid, 2, 1, 4, tsg_level, tsg_fejer2) -call tsgMakeGlobalGrid(grid_II, 2, 1, 4, tsg_level, tsg_fejer2) - -points => tsgGetPoints(grid); i_b = tsgGetNumPoints(grid) -allocate( double_2d_a(1,i_b) ) -double_2d_a(1,:) = exp( -points(1,:)**2 - points(2,:)**2 ) -call tsgLoadNeededPoints(grid, double_2d_a) -call tsgLoadNeededPoints(grid_II, double_2d_a) -call tsgSetAnisotropicRefinement( grid, tsg_iptotal, 10, 1 ) -call tsgSetAnisotropicRefinement( grid_II, tsg_iptotal, 10, 1 ) -deallocate(points, double_2d_a) - -points => tsgGetNeededPoints(grid); i_b = tsgGetNumNeeded(grid) -allocate( double_2d_a(1,i_b) ) -double_2d_a(1,:) = exp( -points(1,:)**2 - points(2,:)**2 ) -call tsgLoadNeededPoints(grid, double_2d_a) -call tsgMergeRefinement(grid_II) -deallocate(points, double_2d_a) - -i_a = 20 -points => tsgGetPoints(grid); i_b = tsgGetNumPoints(grid) -pointsb => tsgGetPoints(grid_II); i_c = tsgGetNumPoints(grid_II) -if ( norm2d(points-pointsb) > 1.d-11 )then - write(*,*) "Mismatch in tsgMergeRefine(): case 2, tsgGetPoints()" - stop 1 -endif -deallocate(points, pointsb) - -allocate(points(2,i_a), double_2d_a(1,i_a)) -rnd => random(2,i_a) -points = -1.d0 + 2.d0 * rnd -call tsgEvaluateBatch(grid_II,points,i_a,double_2d_a) -if ( norm2d(double_2d_a) > 1.d-11 ) then - write(*,*) "Mismatch in tsgMergeRefine(): case 3, tsgEvaluate() not zero" - stop 1 -endif -deallocate(rnd, points, double_2d_a) - -points => tsgGetPoints(grid_II); i_b = tsgGetNumPoints(grid_II) -allocate( double_2d_a(1,i_b) ) -double_2d_a(1,:) = exp( -points(1,:)**2 - points(2,:)**2 ) -call tsgLoadNeededPoints(grid_II, double_2d_a) -deallocate(points, double_2d_a) - -i_a = 30 -allocate( points(2,i_a), double_2d_a(1,i_a), double_2d_b(1,i_a) ) -rnd => random(2,i_a) -points = -1.d0 + 2.d0 * rnd -call tsgEvaluateBatch(grid, points, i_a, double_2d_a) -call tsgEvaluateBatch(grid_II, points, i_a, double_2d_b) -if ( norm2d(double_2d_a-double_2d_b) > 1.d-11 )then - write(*,*) "Mismatch in tsgMergeRefine(): case 3, tsgEvaluate() not equal" - stop 1 -endif -deallocate(points, double_2d_a, double_2d_b, rnd) - - -write(*,*) "Refinement functions: PASS" - - - - -call tsgDeallocateGrid(grid) -call tsgDeallocateGrid(grid_II) - - - - - -contains - -function random(n,m) result(res) - integer :: n, m - double precision, pointer :: res(:,:) - - allocate(res(n,m)) - call random_number(res) -end function - -function norm(x,n) result(res) - integer :: n - double precision :: x(n) - double precision :: res - - res = sqrt(sum(x**2)) -end function - -function norm1d(x) result(res) - double precision :: x(:) - double precision :: res - - res = sqrt(sum(x**2)) -end function - -function norm2d(x) result(res) - double precision :: x(:,:) - double precision :: res - - res = sqrt(sum(x**2)) -end function - -function check_points(points,mustHave,mustNotHave) result(res) - double precision :: points(:), mustHave(:), mustNotHave(:) - logical :: res - integer :: i, j - - res = .false. - do i = 1,size(points) - if ( any(points(i).ne.mustHave) .and. any(points(i).eq.mustNotHave) ) then - res = .true. - exit - endif - enddo -end function - -subroutine print_vector(vec) - double precision :: vec(:) - integer :: i - - do i = 1,size(vec) - write(*,'(f6.3 a)',advance='no') vec(i), " " - enddo - write(*,*) -end subroutine - -subroutine print_array(arr) - double precision :: arr(:,:) - integer :: shp(2), i, j - - shp = shape(arr) - - do i = 1,shp(1) - do j = 1,shp(2) - write(*,'(f6.3 a)',advance='no') arr(i,j), " " - enddo - write(*,*) - enddo -end subroutine - -function sparse_matmul(pntr,indx,val,vec) result(res) - integer :: pntr(:), indx(:) - double precision :: val(:), vec(:) - integer :: i, j - double precision :: res(size(pntr)-1) - - res = 0.d0 - do i = 1,size(pntr)-1 - do j = pntr(i)+1,pntr(i+1) - res(i) = res(i) + val(j) * vec(indx(j)+1) - enddo - enddo -end function - -function Re_complex_matmul(mat,vec) result(res) - double precision :: mat(:,:), vec(:) - integer :: i, j - integer :: shp(2) - double precision :: res(size(mat,1)) - - shp = shape(mat) - - res = 0.d0 - do i = 1,shp(1) - do j = 1,shp(2)/2 - res(i) = res(i) + mat(i,2*j-1) * vec(2*j-1) - mat(i,2*j) * vec(2*j) - enddo - enddo -end function - -end program FORTESTER diff --git a/InterfaceFortran/mpitester.f90 b/InterfaceFortran/mpitester.f90 deleted file mode 100644 index aa21d3ed1..000000000 --- a/InterfaceFortran/mpitester.f90 +++ /dev/null @@ -1,103 +0,0 @@ -!================================================================================================================================================================================== -! Copyright (c) 2018, Miroslav Stoyanov -! -! This file is part of -! Toolkit for Adaptive Stochastic Modeling And Non-Intrusive ApproximatioN: TASMANIAN -! -! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: -! -! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. -! -! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions -! and the following disclaimer in the documentation and/or other materials provided with the distribution. -! -! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse -! or promote products derived from this software without specific prior written permission. -! -! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, -! INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. -! IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, -! OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, -! OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -! -! UT-BATTELLE, LLC AND THE UNITED STATES GOVERNMENT MAKE NO REPRESENTATIONS AND DISCLAIM ALL WARRANTIES, BOTH EXPRESSED AND IMPLIED. -! THERE ARE NO EXPRESS OR IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, OR THAT THE USE OF THE SOFTWARE WILL NOT INFRINGE ANY PATENT, -! COPYRIGHT, TRADEMARK, OR OTHER PROPRIETARY RIGHTS, OR THAT THE SOFTWARE WILL ACCOMPLISH THE INTENDED RESULTS OR THAT THE SOFTWARE OR ITS USE WILL NOT RESULT IN INJURY OR DAMAGE. -! THE USER ASSUMES RESPONSIBILITY FOR ALL LIABILITIES, PENALTIES, FINES, CLAIMS, CAUSES OF ACTION, AND COSTS AND EXPENSES, CAUSED BY, RESULTING FROM OR ARISING OUT OF, -! IN WHOLE OR IN PART THE USE, STORAGE OR DISPOSAL OF THE SOFTWARE. -!================================================================================================================================================================================== -program MPITESTER - use TasmanianSG - use mpi - implicit none - -integer :: ierr, me -integer :: num_points - -type(TasmanianSparseGrid) :: grid, reference_grid - -call MPI_Init(ierr) - -call MPI_Comm_rank(MPI_COMM_WORLD, me, ierr) - -if (ierr /= MPI_SUCCESS) then - write(*,*) "ERROR: failed to initialize MPI" - stop 1 -endif - -!write(*,*) "HELLO CRUEL WORLD from me: ", me, MPI_STATUS_SIZE - -call tsgAllocateGrid(grid) -call tsgAllocateGrid(reference_grid) - -call tsgMakeSequenceGrid(reference_grid, 3, 1, 3, tsg_level, tsg_rleja) - -if (me == 0) then - call tsgMPIGridSend(reference_grid, 1, 11, MPI_COMM_WORLD, ierr) -endif -if (me == 1) then - call tsgMPIGridRecv(grid, 0, 11, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) - if (tsgGetNumPoints(grid) /= tsgGetNumPoints(reference_grid)) then - write(*,*) "ERROR: failed to receive the grid" - stop 1 - endif - call tsgDeallocateGrid(grid) ! reset the grid - call tsgAllocateGrid(grid) -endif - -if (ierr /= MPI_SUCCESS) then - write(*,*) "ERROR: failed at send/recv stage" - stop 1 -endif - -if (me == 2) then - call tsgMPIGridBcast(reference_grid, 2, MPI_COMM_WORLD, ierr) - call tsgCopyGrid(grid, reference_grid) -endif -if (me /= 2) then - call tsgMPIGridBcast(grid, 2, MPI_COMM_WORLD, ierr) -endif - -if (tsgGetNumPoints(grid) /= tsgGetNumPoints(reference_grid)) then - write(*,*) "ERROR: failed to receive the grid" - stop 1 -endif - -if (ierr /= MPI_SUCCESS) then - write(*,*) "ERROR: failed at bcast stage" - stop 1 -endif - -call tsgDeallocateGrid(grid) -call tsgDeallocateGrid(reference_grid) - -if (me == 0) then - write(*,*) "" - write(*,*) " MPI Send/Recv/Bcast Pass" - write(*,*) "" -endif - -call MPI_Finalize(ierr) - -end program MPITESTER diff --git a/InterfaceFortran/tsgC2Fortran.cpp b/InterfaceFortran/tsgC2Fortran.cpp deleted file mode 100644 index 43dbb871b..000000000 --- a/InterfaceFortran/tsgC2Fortran.cpp +++ /dev/null @@ -1,263 +0,0 @@ -/* - * Copyright (c) 2017, Miroslav Stoyanov - * - * This file is part of - * Toolkit for Adaptive Stochastic Modeling And Non-Intrusive ApproximatioN: TASMANIAN - * - * Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: - * - * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. - * - * 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions - * and the following disclaimer in the documentation and/or other materials provided with the distribution. - * - * 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse - * or promote products derived from this software without specific prior written permission. - * - * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, - * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. - * IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, - * OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, - * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, - * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - * - * UT-BATTELLE, LLC AND THE UNITED STATES GOVERNMENT MAKE NO REPRESENTATIONS AND DISCLAIM ALL WARRANTIES, BOTH EXPRESSED AND IMPLIED. - * THERE ARE NO EXPRESS OR IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, OR THAT THE USE OF THE SOFTWARE WILL NOT INFRINGE ANY PATENT, - * COPYRIGHT, TRADEMARK, OR OTHER PROPRIETARY RIGHTS, OR THAT THE SOFTWARE WILL ACCOMPLISH THE INTENDED RESULTS OR THAT THE SOFTWARE OR ITS USE WILL NOT RESULT IN INJURY OR DAMAGE. - * THE USER ASSUMES RESPONSIBILITY FOR ALL LIABILITIES, PENALTIES, FINES, CLAIMS, CAUSES OF ACTION, AND COSTS AND EXPENSES, CAUSED BY, RESULTING FROM OR ARISING OUT OF, - * IN WHOLE OR IN PART THE USE, STORAGE OR DISPOSAL OF THE SOFTWARE. - */ - -#ifndef __TSG_C2FORT_CPP -#define __TSG_C2FORT_CPP - -#include "TasmanianAddons.hpp" - -using std::cerr; -using std::endl; - -using namespace TasGrid; - -extern "C" void tsgc2fstr_(int *length, const char *str); -extern "C" void tsgc2fvec_(int *length, double *vect); -extern "C" void tsgc2fmat_(int *rows, int *cols, double *mat); - -extern "C"{ - -//////////////////////////////////////////////////////////////////////// -// MAIN INTERFACE -//////////////////////////////////////////////////////////////////////// - -void tsggvm_(int *ver){ *ver = TasmanianSparseGrid::getVersionMajor(); } -void tsggvn_(int *ver){ *ver = TasmanianSparseGrid::getVersionMinor(); } -void tsggli_(){ - // this is clumsy, may have to think of something else - const char *lic = TasmanianSparseGrid::getLicense(); - int l = 0; - while(lic[l] != '\0') l++; - tsgc2fstr_(&l, lic); -} - -// read/write -void tsgwri_(TasmanianSparseGrid **grid, const char *filename, int *binary){ (*grid)->write(filename, (*binary != 0)); } -void tsgrea_(TasmanianSparseGrid **grid, const char *filename, int *status){ - try{ - (*grid)->read(filename); - *status = 1; - }catch(std::runtime_error &e){ - cerr << e.what() << endl; - *status = 0; - } -} - -void tsgalloc_(TasmanianSparseGrid **grid){ - (*grid) = new TasmanianSparseGrid(); -} -void tsgfree_(TasmanianSparseGrid **grid){ - delete (*grid); -} - -// create -void tsgmg_(TasmanianSparseGrid **grid, int *dimensions, int *outputs, int *depth, int *type, int *rule, int *opt_flags, - const int *aniso_weights, double *alpha, double *beta, const char *customRuleFilename, const int *llimits){ - const int *aw = (opt_flags[0] != 0 ) ? aniso_weights : nullptr; - double al = (opt_flags[1] != 0 ) ? *alpha : 0.0; - double be = (opt_flags[2] != 0 ) ? *beta : 0.0; - const char *cfn = (opt_flags[3] != 0 ) ? customRuleFilename : nullptr; - const int *ll = (opt_flags[4] != 0 ) ? llimits : nullptr; - - (*grid)->makeGlobalGrid(*dimensions, *outputs, *depth, IO::getDepthTypeInt(*type), IO::getRuleInt(*rule), aw, al, be, cfn, ll); -} -void tsgms_(TasmanianSparseGrid **grid, int *dimensions, int *outputs, int *depth, int *type, int *rule, int *opt_flags, - const int *aniso_weights, const int *llimits){ - const int *aw = (opt_flags[0] != 0 ) ? aniso_weights : nullptr; - const int *ll = (opt_flags[1] != 0 ) ? llimits : nullptr; - - (*grid)->makeSequenceGrid(*dimensions, *outputs, *depth, IO::getDepthTypeInt(*type), IO::getRuleInt(*rule), aw, ll); -} -void tsgml_(TasmanianSparseGrid **grid, int *dimensions, int *outputs, int *depth, int *opt_flags, - int *order, int *rule, const int *llimits){ - int ru = (opt_flags[0] != 0) ? *rule : 1; - int ord = (opt_flags[1] != 0) ? *order : 1; - const int *ll = (opt_flags[2] != 0) ? llimits : nullptr; - - (*grid)->makeLocalPolynomialGrid(*dimensions, *outputs, *depth, ord, IO::getRuleInt(ru), ll); -} -void tsgmw_(TasmanianSparseGrid **grid, int *dimensions, int *outputs, int *depth, int *opt_flags, int *order, const int *llimits){ - int ord = (opt_flags[0] != 0) ? *order : 1; - const int *ll = (opt_flags[1] != 0) ? llimits : nullptr; - - (*grid)->makeWaveletGrid(*dimensions, *outputs, *depth, ord, ll); -} -void tsgmf_(TasmanianSparseGrid **grid, int *dimensions, int *outputs, int *depth, int *type, int *opt_flags, const int *aniso_weights, const int *llimits){ - const int *aw = (opt_flags[0] != 0 ) ? aniso_weights : nullptr; - const int *ll = (opt_flags[1] != 0 ) ? llimits : nullptr; - - (*grid)->makeFourierGrid(*dimensions, *outputs, *depth, IO::getDepthTypeInt(*type), aw, ll); -} - -// copy/updare -void tsgcp_(TasmanianSparseGrid **grid, TasmanianSparseGrid **source){ (*grid)->copyGrid(*source); } - -void tsgug_(TasmanianSparseGrid **grid, int *depth, int *type, int *opt_flags, const int *anisotropic_weights){ - const int *aw = (opt_flags[0] != 0) ? anisotropic_weights : nullptr; - (*grid)->updateGlobalGrid(*depth, IO::getDepthTypeInt(*type), aw); -} -void tsgus_(TasmanianSparseGrid **grid, int *depth, int *type, int *opt_flags, const int *anisotropic_weights){ - const int *aw = (opt_flags[0] != 0) ? anisotropic_weights : nullptr; - (*grid)->updateSequenceGrid(*depth, IO::getDepthTypeInt(*type), aw); -} - -// getAlpha/Beta/Order/Dims/Outs/Rule // -void tsggal_(TasmanianSparseGrid **grid, double *alpha){ *alpha = (*grid)->getAlpha(); } -void tsggbe_(TasmanianSparseGrid **grid, double *beta){ *beta = (*grid)->getBeta(); } -void tsggor_(TasmanianSparseGrid **grid, int *order){ *order = (*grid)->getOrder(); } -void tsggnd_(TasmanianSparseGrid **grid, int *dims){ *dims = (*grid)->getNumDimensions(); } -void tsggno_(TasmanianSparseGrid **grid, int *outs){ *outs = (*grid)->getNumOutputs(); } -void tsggru_(TasmanianSparseGrid **grid, int *rule){ *rule = IO::getRuleInt((*grid)->getRule()); } - -// getNumNeeded/Loaded/Points -void tsggnn_(TasmanianSparseGrid **grid, int *num_points){ *num_points = (*grid)->getNumNeeded(); } -void tsggnl_(TasmanianSparseGrid **grid, int *num_points){ *num_points = (*grid)->getNumLoaded(); } -void tsggnp_(TasmanianSparseGrid **grid, int *num_points){ *num_points = (*grid)->getNumPoints(); } - -// getLoaded/Needed/Points -void tsgglp_(TasmanianSparseGrid **grid, double *points){ (*grid)->getLoadedPoints(points); } -void tsggdp_(TasmanianSparseGrid **grid, double *points){ (*grid)->getNeededPoints(points); } -void tsggpp_(TasmanianSparseGrid **grid, double *points){ (*grid)->getPoints(points); } - -// getQuadrature/InterpolationWeights -void tsggqw_(TasmanianSparseGrid **grid, double *weights){ (*grid)->getQuadratureWeights(weights); } -void tsggiw_(TasmanianSparseGrid **grid, const double *x, double *weights){ (*grid)->getInterpolationWeights(x,weights); } - -// set/is/clear/getDomainTransform -void tsgsdt_(TasmanianSparseGrid **grid, const double *transform_a, const double *transform_b){ (*grid)->setDomainTransform(transform_a, transform_b); } -void tsgidt_(TasmanianSparseGrid **grid, int *result){ *result = ((*grid)->isSetDomainTransfrom()) ? 1 : 0; } -void tsgcdt_(TasmanianSparseGrid **grid){ (*grid)->clearDomainTransform(); } -void tsggdt_(TasmanianSparseGrid **grid, double *transform_a, double *transform_b){ (*grid)->getDomainTransform(transform_a, transform_b); } - -// loadNeededPoints -void tsglnp_(TasmanianSparseGrid **grid, const double *vals){ (*grid)->loadNeededPoints(vals); } - -// evaluate/Fast/Batch/integrate -void tsgeva_(TasmanianSparseGrid **grid, const double *x, double *y){ (*grid)->evaluate(x, y); } -void tsgevf_(TasmanianSparseGrid **grid, const double *x, double *y){ (*grid)->evaluateFast(x, y); } -void tsgevb_(TasmanianSparseGrid **grid, const double *x, int *num_x, double *y){ (*grid)->evaluateBatch(x, (*num_x), y); } -void tsgint_(TasmanianSparseGrid **grid, double *q){ (*grid)->integrate(q); } - -// hierarchical functions/coefficients -void tsgehf_(TasmanianSparseGrid **grid, const double *x, int *num_x, double *y){ - (*grid)->evaluateHierarchicalFunctions(x, (*num_x), y);} -void tsgehs_(TasmanianSparseGrid **grid, const double *x, int *num_x, int *pntr, int *indx, double *vals){ - (*grid)->evaluateSparseHierarchicalFunctionsStatic(x, *num_x, pntr, indx, vals);} -void tsgehz_(TasmanianSparseGrid **grid, const double *x, int *num_x, int *num_nz){ - *num_nz = (*grid)->evaluateSparseHierarchicalFunctionsGetNZ(x, *num_x);} -void tsgghc_(TasmanianSparseGrid **grid, double *c){ - const double *cc = (*grid)->getHierarchicalCoefficients(); - (*grid)->isFourier() ? std::copy(cc, cc + 2 * (*grid)->getNumPoints() * (*grid)->getNumOutputs(), c) - : std::copy(cc, cc + (*grid)->getNumPoints() * (*grid)->getNumOutputs(), c); -} -void tsgshc_(TasmanianSparseGrid **grid, double *c){(*grid)->setHierarchicalCoefficients(c);} -void tsghsu_(TasmanianSparseGrid **grid, double *c){ - std::vector sup = (*grid)->getHierarchicalSupport(); - std::copy(sup.begin(), sup.end(), c); -} - -// setAnisotropic/Surplus/Refinement -void tsgsar_(TasmanianSparseGrid **grid, int *type, int *min_growth, int *output, int *opt_flags, const int *llimits){ - const int *ll = (opt_flags[0] != 0) ? llimits : nullptr; - (*grid)->setAnisotropicRefinement(IO::getDepthTypeInt(*type), *min_growth, *output, ll); -} -void tsgeac_(TasmanianSparseGrid **grid, int *type, int *output, int *result){ - auto coeff = (*grid)->estimateAnisotropicCoefficients(IO::getDepthTypeInt(*type), *output); - int num_coeff = (*grid)->getNumDimensions(); - if ((*type == 2) || (*type == 4) || (*type == 6)) num_coeff *= 2; - for(int i=0; isetSurplusRefinement(*tol, *output, ll); } -void tsgshr_(TasmanianSparseGrid **grid, double *tol, int *type, int *opt_flags, int *output, const int *llimits){ - int theout = (opt_flags[0] != 0) ? *output : -1; - const int *ll = (opt_flags[1] != 0) ? llimits : nullptr; - - (*grid)->setSurplusRefinement(*tol, IO::getTypeRefinementInt(*type), theout, ll); -} -void tsgcre_(TasmanianSparseGrid **grid){ (*grid)->clearRefinement(); } -void tsgmre_(TasmanianSparseGrid **grid){ (*grid)->mergeRefinement(); } - -// set/is/clear/getConformalTransform -void tsgsca_(TasmanianSparseGrid **grid, int *trunc){ (*grid)->setConformalTransformASIN(Utils::copyArray(trunc, (*grid)->getNumDimensions())); } -void tsgica_(TasmanianSparseGrid **grid, int *result){ *result = ((*grid)->isSetConformalTransformASIN()) ? 1 : 0; } -void tsgcct_(TasmanianSparseGrid **grid){ (*grid)->clearConformalTransform(); } -void tsggca_(TasmanianSparseGrid **grid, int *trunc){ - auto truncation_vector = (*grid)->getConformalTransformASIN(); - std::copy(truncation_vector.begin(), truncation_vector.end(), trunc); -} - -// isGlobal/Sequence/LocalPolynomial/Wavelet/Fourier -void tsgisg_(TasmanianSparseGrid **grid, int *status){*status = ((*grid)->isGlobal()) ? 1 : 0;} -void tsgiss_(TasmanianSparseGrid **grid, int *status){*status = ((*grid)->isSequence()) ? 1 : 0;} -void tsgisl_(TasmanianSparseGrid **grid, int *status){*status = ((*grid)->isLocalPolynomial()) ? 1 : 0;} -void tsgisw_(TasmanianSparseGrid **grid, int *status){*status = ((*grid)->isWavelet()) ? 1 : 0;} -void tsgisf_(TasmanianSparseGrid **grid, int *status){*status = ((*grid)->isFourier()) ? 1 : 0;} - -// print stats -void tsgpri_(TasmanianSparseGrid **grid){ (*grid)->printStats(); } - -// get/enableAcceleration -void tsgacc_(TasmanianSparseGrid **grid, int *acc){ - (*grid)->enableAcceleration(AccelerationMeta::getIOIntAcceleration(*acc)); -} -void tsggac_(TasmanianSparseGrid **grid, int *acc){ - TypeAcceleration accel = (*grid)->getAccelerationType(); - *acc = AccelerationMeta::getIOAccelerationInt(accel); -} -void tsgsgi_(TasmanianSparseGrid **grid, int *gpuID){ (*grid)->setGPUID(*gpuID); } -void tsgggi_(TasmanianSparseGrid **grid, int *gpuID){ *gpuID = (*grid)->getGPUID(); } -void tsggng_(int *gpus){ *gpus = TasmanianSparseGrid::getNumGPUs(); } -void tsgggm_(int *gpuID, int *mem){ *mem = TasmanianSparseGrid::getGPUMemory(*gpuID); } - -// MPI send/recv/bcast -#ifdef Tasmanian_ENABLE_MPI -void tsgmpigsend_(TasmanianSparseGrid **grid, int *dest, int *tag, MPI_Fint *comm, int *ierr){ - *ierr = MPIGridSend(**grid, *dest, *tag, MPI_Comm_f2c(*comm)); -} -void tsgmpigrecv_(TasmanianSparseGrid **grid, int *src, int *tag, MPI_Fint *comm, MPI_Fint *stat, int *ierr){ - MPI_Status status; - *ierr = MPI_Status_f2c(stat, &status); - if (*ierr != MPI_SUCCESS) return; - *ierr = MPIGridRecv(**grid, *src, *tag, MPI_Comm_f2c(*comm), &status); -} -void tsgmpigcast_(TasmanianSparseGrid **grid, int *root, MPI_Fint *comm, int *ierr){ - *ierr = MPIGridBcast(**grid, *root, MPI_Comm_f2c(*comm)); -} -#else -void tsgmpigsend_(TasmanianSparseGrid**, int*, int*, int*, int*){} -void tsgmpigrecv_(TasmanianSparseGrid**, int*, int*, int*, int*, int*){} -void tsgmpigcast_(TasmanianSparseGrid**, int*, int*, int*){} -#endif - -} -#endif diff --git a/InterfaceFortran/tsgC2FortranBridge.f90 b/InterfaceFortran/tsgC2FortranBridge.f90 deleted file mode 100644 index d1b7d1025..000000000 --- a/InterfaceFortran/tsgC2FortranBridge.f90 +++ /dev/null @@ -1,65 +0,0 @@ -!================================================================================================================================================================================== -! Copyright (c) 2017, Miroslav Stoyanov -! -! This file is part of -! Toolkit for Adaptive Stochastic Modeling And Non-Intrusive ApproximatioN: TASMANIAN -! -! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: -! -! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. -! -! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions -! and the following disclaimer in the documentation and/or other materials provided with the distribution. -! -! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse -! or promote products derived from this software without specific prior written permission. -! -! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, -! INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. -! IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, -! OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, -! OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, -! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -! -! UT-BATTELLE, LLC AND THE UNITED STATES GOVERNMENT MAKE NO REPRESENTATIONS AND DISCLAIM ALL WARRANTIES, BOTH EXPRESSED AND IMPLIED. -! THERE ARE NO EXPRESS OR IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE, OR THAT THE USE OF THE SOFTWARE WILL NOT INFRINGE ANY PATENT, -! COPYRIGHT, TRADEMARK, OR OTHER PROPRIETARY RIGHTS, OR THAT THE SOFTWARE WILL ACCOMPLISH THE INTENDED RESULTS OR THAT THE SOFTWARE OR ITS USE WILL NOT RESULT IN INJURY OR DAMAGE. -! THE USER ASSUMES RESPONSIBILITY FOR ALL LIABILITIES, PENALTIES, FINES, CLAIMS, CAUSES OF ACTION, AND COSTS AND EXPENSES, CAUSED BY, RESULTING FROM OR ARISING OUT OF, -! IN WHOLE OR IN PART THE USE, STORAGE OR DISPOSAL OF THE SOFTWARE. -!================================================================================================================================================================================== -!SUBROUTINE tsgc2fmat(r, c, M) -! USE TasmanianSG, only: tsgReceiveMatrix -! INTEGER, INTENT(IN) :: r, c -! DOUBLE PRECISION, INTENT(IN) :: M(r*c) -! CALL tsgReceiveMatrix(r, c, M) -!END SUBROUTINE tsgc2fmat -!! ========================================= -!SUBROUTINE tsgc2fvec(s, V) -! USE TasmanianSG, only: tsgReceiveVector -! INTEGER, INTENT(IN) :: s -! DOUBLE PRECISION, INTENT(IN) :: V(s) -! CALL tsgReceiveVector(s, V) -!END SUBROUTINE tsgc2fvec -!! ========================================= -!SUBROUTINE tsgc2fdouble(d) -! USE TasmanianSG, only: tsgReceiveScalar -! DOUBLE PRECISION, INTENT(IN) :: d -! CALL tsgReceiveScalar(d) -!END SUBROUTINE tsgc2fdouble -!! ========================================= -!SUBROUTINE tsgc2fint(i) -! USE TasmanianSG, only: tsgReceiveInt -! INTEGER, INTENT(IN) :: i -! CALL tsgReceiveInt(i) -!END SUBROUTINE tsgc2fint -!! ========================================= - -! ========================================= -SUBROUTINE tsgc2fstr(l, S) - USE TasmanianSG, only: tsgReceiveString - IMPLICIT NONE - INTEGER, INTENT(IN) :: l - CHARACTER, INTENT(IN) :: S(l) - CALL tsgReceiveString(l, S) -END SUBROUTINE tsgc2fstr -! =========================================