From 3a22bf7e27f0011922909af24f676a1c21997d0e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 13 Jan 2012 16:18:16 +0000 Subject: [PATCH] changed fftw from legacy fortran to new (2003) fortran (calling c routines directly) renamed "steps" consequently to "incs" moved kdtree2 to math.f90, put original source to private folder --- code/DAMASK_spectral.f90 | 390 +++---- code/FEsolving.f90 | 8 +- code/IO.f90 | 4 +- code/core_modules.f90 | 2 +- code/damask.core.pyf | 17 +- code/fftw3.f03 | 1819 +++++++++++++++++++++++++++++++ code/kdtree2.f90 | 1885 -------------------------------- code/makefile | 21 +- code/math.f90 | 2186 ++++++++++++++++++++++++++++++++++++-- code/numerics.f90 | 2 +- 10 files changed, 4142 insertions(+), 2192 deletions(-) create mode 100644 code/fftw3.f03 delete mode 100644 code/kdtree2.f90 diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index 7d0a8acb4..9f18a1c29 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -38,20 +38,21 @@ program DAMASK_spectral !******************************************************************** use DAMASK_interface - use prec, only: pInt, pReal, DAMASK_NaN + use prec, only: pInt, pReal, DAMASK_NaN use IO - use debug, only: debug_spectral, & - debug_spectralGeneral, & - debug_spectralDivergence, & - debug_spectralRestart + use debug, only: debug_spectral, & + debug_spectralGeneral, & + debug_spectralDivergence, & + debug_spectralRestart, & + debug_spectralFFTW use math use kdtree2_module - use mesh, only: mesh_ipCenterOfGravity - use CPFEM, only: CPFEM_general, CPFEM_initAll - use FEsolving, only: restartWrite, restartReadStep - use numerics, only: err_div_tol, err_stress_tolrel , rotation_tol,& - itmax, memory_efficient, DAMASK_NumThreadsInt, divergence_correction, & - fftw_planner_flag, fftw_timelimit + use mesh, only: mesh_ipCenterOfGravity + use CPFEM, only: CPFEM_general, CPFEM_initAll + use FEsolving, only: restartWrite, restartReadInc + use numerics, only: err_div_tol, err_stress_tolrel , rotation_tol,& + itmax, memory_efficient, DAMASK_NumThreadsInt, divergence_correction, & + fftw_planner_flag, fftw_timelimit use homogenization, only: materialpoint_sizeResults, materialpoint_results !$ use OMP_LIB ! the openMP function library @@ -70,23 +71,23 @@ program DAMASK_spectral character(len=1024) :: path, line, keyword logical :: gotResolution =.false., gotDimension =.false., gotHomogenization = .false. type bc_type - real(pReal), dimension (3,3) :: deformation, & ! applied velocity gradient or time derivative of deformation gradient - stress, & ! stress BC (if applicable) - rotation ! rotation of BC (if applicable) - real(pReal) :: timeIncrement, & ! length of increment - temperature ! isothermal starting conditions - integer(pInt) :: steps, & ! number of steps - outputfrequency, & ! frequency of result writes - restartfrequency, & ! frequency of restart writes - logscale ! linear/logaritmic time step flag - logical :: followFormerTrajectory,& ! follow trajectory of former loadcase - velGradApplied ! decide wether velocity gradient or fdot is given - logical, dimension(3,3) :: maskDeformation, & ! mask of boundary conditions - maskStress - logical, dimension(9) :: maskStressVector ! linear mask of boundary conditions + real(pReal), dimension (3,3) :: deformation = 0.0_pReal, & ! applied velocity gradient or time derivative of deformation gradient + stress = 0.0_pReal, & ! stress BC (if applicable) + rotation = math_I3 ! rotation of BC (if applicable) + real(pReal) :: time = 0.0_pReal, & ! length of increment + temperature = 300_pReal ! isothermal starting conditions + integer(pInt) :: incs = 0_pInt, & ! number of increments + outputfrequency = 1_pInt, & ! frequency of result writes + restartfrequency = 0_pInt, & ! frequency of restart writes + logscale = 0_pInt ! linear/logaritmic time inc flag + logical :: followFormerTrajectory = .true., & ! follow trajectory of former loadcase + velGradApplied = .false. ! decide wether velocity gradient or fdot is given + logical, dimension(3,3) :: maskDeformation = .false., & ! mask of deformation boundary conditions + maskStress = .false. ! mask of stress boundary conditions + logical, dimension(9) :: maskStressVector = .false. ! linear mask of boundary conditions end type + type(bc_type), allocatable, dimension(:) :: bc - type(bc_type) :: bc_default character(len=6) :: loadcase_string ! variables storing information from geom file @@ -95,6 +96,7 @@ program DAMASK_spectral integer(pInt) :: Npoints,& ! number of Fourier points homog ! homogenization scheme used integer(pInt), dimension(3) :: res = 1_pInt ! resolution (number of Fourier points) in each direction + integer(pInt) :: res1_red ! stress, stiffness and compliance average etc. real(pReal), dimension(3,3) :: pstress, pstress_av, defgrad_av_lab, & @@ -109,7 +111,7 @@ program DAMASK_spectral integer(pInt) :: size_reduced = 0.0_pReal ! number of stress BCs ! pointwise data - real(pReal), dimension(:,:,:,:,:), allocatable :: workfft, defgrad, defgradold + real(pReal), dimension(:,:,:,:,:), allocatable :: defgrad, defgradold real(pReal), dimension(:,:,:,:), allocatable :: coordinates real(pReal), dimension(:,:,:), allocatable :: temperature @@ -118,8 +120,14 @@ program DAMASK_spectral real(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat ! gamma operator (field) for spectral method real(pReal), dimension(:,:,:,:), allocatable :: xi ! wave vector field for divergence and for gamma operator integer(pInt), dimension(3) :: k_s - integer*8, dimension(3) :: fftw_plan ! plans for fftw (forward and backward) - + type(C_PTR) :: data_fftw, fftw_stress, fftw_fluctuation + real(pReal), dimension(:,:,:,:,:), pointer :: data_real + complex(pReal), dimension(:,:,:,:,:), pointer :: data_complex +!debugging (proof of correct transformation) + type(C_PTR) :: fftw_debug, fftw_debug_forward, fftw_debug_backward + real(pReal), dimension(:,:,:), pointer :: fftw_debug_real + complex(pReal), dimension(:,:,:), pointer :: fftw_debug_complex + ! variables for regriding real(pReal), dimension(:,:,:,:) ,allocatable :: deformed_small real(pReal), dimension(:,:) ,allocatable :: deformed_large @@ -131,42 +139,36 @@ program DAMASK_spectral real(pReal) :: time = 0.0_pReal, time0 = 0.0_pReal, timeinc ! elapsed time, begin of interval, time interval real(pReal) :: guessmode, err_div, err_stress, err_stress_tol, p_hat_avg complex(pReal) :: err_div_avg_complex - complex(pReal), dimension(3) :: divergence_complex complex(pReal), parameter :: img = cmplx(0.0_pReal,1.0_pReal) real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal + complex(pReal), dimension(3) :: temp3_Complex complex(pReal), dimension(3,3) :: temp33_Complex - real(pReal), dimension(3,3) :: temp33_Real + real(pReal), dimension(3,3) :: temp33_Real integer(pInt) :: i, j, k, l, m, n, p, errorID - integer(pInt) :: N_Loadcases, loadcase, step, iter, ielem, CPFEM_mode, & - ierr, notConvergedCounter = 0_pInt, totalStepsCounter = 0_pInt - logical :: errmatinv, regrid = .false. + integer(pInt) :: N_Loadcases, loadcase, inc, iter, ielem, CPFEM_mode, & + ierr, notConvergedCounter = 0_pInt, totalIncsCounter = 0_pInt + logical :: errmatinv real(pReal) :: defgradDet, defgradDetMax, defgradDetMin real(pReal) :: correctionFactor integer(pInt), dimension(3) :: cutting_freq ! --- debugging variables - real(pReal), dimension(:,:,:,:), allocatable :: divergence + type(C_PTR) :: divergence + real(pReal), dimension(:,:,:,:), pointer :: divergence_real + complex(pReal), dimension(:,:,:,:), pointer :: divergence_complex real(pReal) :: p_real_avg, err_div_max, err_real_div_avg, err_real_div_max - logical :: debugGeneral, debugDivergence, debugRestart + logical :: debugGeneral, debugDivergence, debugRestart, debugFFTW + type(C_PTR) :: fftw_divergence ! plan for fftw backward transform of divergence + integer(pInt) :: row, column -! --- initialize default value for loadcase - bc_default%steps = 0_pInt; bc_default%timeIncrement = 0.0_pReal - bc_default%temperature = 300.0_pReal; bc_default%logscale = 0_pInt - bc_default%outputfrequency = 1_pInt; bc_default%restartfrequency = 0_pInt - bc_default%deformation = zeroes; bc_default%stress = zeroes - bc_default%maskDeformation = .false.; bc_default%maskStress = .false. - bc_default%velGradApplied = .false.; bc_default%maskStressVector = .false. - bc_default%followFormerTrajectory = .true. - bc_default%rotation = math_I3 ! assume no rotation - ! --- initializing model size independed parameters !$ call omp_set_num_threads(DAMASK_NumThreadsInt) ! set number of threads for parallel execution set by DAMASK_NUM_THREADS call DAMASK_interface_init() print '(a)', '' - print '(a,a)', ' <<<+- DAMASK_spectral init -+>>>' - print '(a,a)', ' $Id$' + print '(a)', ' <<<+- DAMASK_spectral init -+>>>' + print '(a)', ' $Id$' print '(a)', '' print '(a,a)', ' Working Directory: ',trim(getSolverWorkingDirectoryName()) print '(a,a)', ' Solver Job Name: ',trim(getSolverJobName()) @@ -182,13 +184,13 @@ program DAMASK_spectral positions = IO_stringPos(line,maxNchunksLoadcase) do i = 1_pInt, maxNchunksLoadcase, 1_pInt ! reading compulsory parameters for loadcase select case (IO_lc(IO_stringValue(line,positions,i))) - case('l', 'velocitygrad', 'velgrad','velocitygradient') + case('l','velocitygrad','velgrad','velocitygradient') N_l = N_l + 1_pInt case('fdot') N_Fdot = N_Fdot + 1_pInt - case('t', 'time', 'delta') + case('t','time','delta') N_t = N_t + 1_pInt - case('n', 'incs', 'increments', 'steps', 'logincs', 'logsteps') + case('n','incs','increments','steps','logincs','logsteps') N_n = N_n + 1_pInt end select enddo ! count all identifiers to allocate memory and do sanity check @@ -207,7 +209,6 @@ program DAMASK_spectral read(myUnit,'(a1024)',END = 101) line if (IO_isBlank(line)) cycle ! skip empty lines loadcase = loadcase + 1_pInt - bc(loadcase) = bc_default positions = IO_stringPos(line,maxNchunksLoadcase) do j = 1_pInt,maxNchunksLoadcase select case (IO_lc(IO_stringValue(line,positions,j))) @@ -224,7 +225,7 @@ program DAMASK_spectral enddo bc(loadcase)%maskDeformation = transpose(reshape(temp_maskVector,(/3,3/))) bc(loadcase)%deformation = math_plain9to33(temp_valueVector) - case('p', 'pk1', 'piolakirchhoff', 'stress') + case('p','pk1','piolakirchhoff','stress') temp_valueVector = 0.0_pReal forall (k = 1_pInt:9_pInt) bc(loadcase)%maskStressVector(k) = IO_stringValue(line,positions,j+k) /= '*' do k = 1_pInt,9_pInt @@ -233,13 +234,13 @@ program DAMASK_spectral bc(loadcase)%maskStress = transpose(reshape(bc(loadcase)%maskStressVector,(/3,3/))) bc(loadcase)%stress = math_plain9to33(temp_valueVector) case('t','time','delta') ! increment time - bc(loadcase)%timeIncrement = IO_floatValue(line,positions,j+1_pInt) + bc(loadcase)%time = IO_floatValue(line,positions,j+1_pInt) case('temp','temperature') ! starting temperature bc(loadcase)%temperature = IO_floatValue(line,positions,j+1_pInt) case('n','incs','increments','steps') ! number of increments - bc(loadcase)%steps = IO_intValue(line,positions,j+1_pInt) + bc(loadcase)%incs = IO_intValue(line,positions,j+1_pInt) case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling) - bc(loadcase)%steps = IO_intValue(line,positions,j+1_pInt) + bc(loadcase)%incs = IO_intValue(line,positions,j+1_pInt) bc(loadcase)%logscale = 1_pInt case('f','freq','frequency','outputfreq') ! frequency of result writings bc(loadcase)%outputfrequency = IO_intValue(line,positions,j+1_pInt) @@ -319,12 +320,13 @@ program DAMASK_spectral enddo close(myUnit) if (.not.(gotDimension .and. gotHomogenization .and. gotResolution)) call IO_error(error_ID = 45_pInt) - + if(mod(res(1),2_pInt)/=0_pInt .or.& mod(res(2),2_pInt)/=0_pInt .or.& (mod(res(3),2_pInt)/=0_pInt .and. res(3)/= 1_pInt)) call IO_error(error_ID = 103_pInt) - + res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) Npoints = res(1)*res(2)*res(3) + wgt = 1.0_pReal/real(Npoints, pReal) ! --- initialization of CPFEM_general (= constitutive law) call CPFEM_initAll(bc(1)%temperature,1_pInt,1_pInt) @@ -333,6 +335,7 @@ program DAMASK_spectral debugGeneral = iand(debug_spectral,debug_spectralGeneral) > 0_pInt debugDivergence = iand(debug_spectral,debug_spectralDivergence) > 0_pInt debugRestart = iand(debug_spectral,debug_spectralRestart) > 0_pInt + debugFFTW = iand(debug_spectral,debug_spectralFFTW) > 0_pInt ! --- output of geometry print '(a)', '' @@ -350,7 +353,7 @@ program DAMASK_spectral print '(a,a)', 'loadcase file: ',trim(getLoadcaseName()) if (bc(1)%followFormerTrajectory) then - call IO_warning(warning_ID = 33_pInt) ! cannot guess along trajectory for first step of first loadcase + call IO_warning(warning_ID = 33_pInt) ! cannot guess along trajectory for first inc of first loadcase bc(1)%followFormerTrajectory = .false. endif @@ -381,8 +384,8 @@ program DAMASK_spectral if (any(bc(loadcase)%rotation /= math_I3)) & print '(a,3(3(f12.6,x)/)$)','rotation of loadframe:',math_transpose3x3(bc(loadcase)%rotation) print '(a,f12.6)','temperature:',bc(loadcase)%temperature - print '(a,f12.6)','time: ',bc(loadcase)%timeIncrement - print '(a,i5)','steps: ',bc(loadcase)%steps + print '(a,f12.6)','time: ',bc(loadcase)%time + print '(a,i5)' ,'increments: ',bc(loadcase)%incs print '(a,i5)','output frequency: ',bc(loadcase)%outputfrequency print '(a,i5)','restart frequency: ',bc(loadcase)%restartfrequency @@ -395,21 +398,23 @@ program DAMASK_spectral > reshape(spread(rotation_tol,1,9),(/3,3/)))& .or. abs(math_det3x3(bc(loadcase)%rotation)) > 1.0_pReal + rotation_tol) & errorID = 46_pInt ! given rotation matrix contains strain - if (bc(loadcase)%timeIncrement < 0.0_pReal) errorID = 34_pInt ! negative time increment - if (bc(loadcase)%steps < 1_pInt) errorID = 35_pInt ! non-positive increment count + if (bc(loadcase)%time < 0.0_pReal) errorID = 34_pInt ! negative time increment + if (bc(loadcase)%incs < 1_pInt) errorID = 35_pInt ! non-positive incs count if (bc(loadcase)%outputfrequency < 1_pInt) errorID = 36_pInt ! non-positive result frequency if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) enddo ! Initialization of fftw (see manual on fftw.org for more details) + if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=102) #ifdef _OPENMP if(DAMASK_NumThreadsInt > 0_pInt) then - call dfftw_init_threads(ierr) + ierr = fftw_init_threads() if (ierr == 0_pInt) call IO_error(error_ID = 104_pInt) - call dfftw_plan_with_nthreads(DAMASK_NumThreadsInt) - endif + call fftw_plan_with_nthreads(DAMASK_NumThreadsInt) + endif #endif - call dfftw_set_timelimit(fftw_timelimit) + call fftw_set_timelimit(fftw_timelimit) + !************************************************************* ! Loop over loadcases defined in the loadcase file do loadcase = 1_pInt, N_Loadcases @@ -418,83 +423,92 @@ program DAMASK_spectral if (bc(loadcase)%followFormerTrajectory) then ! continue to guess along former trajectory where applicable guessmode = 1.0_pReal else - guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first step + guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first inc endif - mask_defgrad = merge(ones,zeroes,bc(loadcase)%maskDeformation) - mask_stress = merge(ones,zeroes,bc(loadcase)%maskStress) + mask_defgrad = merge(ones,zeroes,bc(loadcase)%maskDeformation) + mask_stress = merge(ones,zeroes,bc(loadcase)%maskStress) size_reduced = count(bc(loadcase)%maskStressVector) allocate (c_reduced(size_reduced,size_reduced)); c_reduced = 0.0_pReal allocate (s_reduced(size_reduced,size_reduced)); s_reduced = 0.0_pReal - timeinc = bc(loadcase)%timeIncrement/bc(loadcase)%steps ! only valid for given linear time scale. will be overwritten later in case logarithmic scale is used + timeinc = bc(loadcase)%time/bc(loadcase)%incs ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used + fDot = bc(loadcase)%deformation ! only valid for given fDot. will be overwritten later in case L is given !************************************************************* -! loop oper steps defined in input file for current loadcase - do step = 1_pInt, bc(loadcase)%steps +! loop oper incs defined in input file for current loadcase + do inc = 1_pInt, bc(loadcase)%incs !************************************************************* ! forwarding time - if (bc(loadcase)%logscale == 1_pInt) then ! logarithmic scale - if (loadcase == 1_pInt) then ! 1st loadcase of logarithmic scale - if (step == 1_pInt) then ! 1st step of 1st loadcase of logarithmic scale - timeinc = bc(1)%timeIncrement*(2.0_pReal**real( 1_pInt-bc(1)%steps ,pReal)) ! assume 1st step is equal to 2nd - else ! not-1st step of 1st loadcase of logarithmic scale - timeinc = bc(1)%timeIncrement*(2.0_pReal**real(step-1_pInt-bc(1)%steps ,pReal)) + if (bc(loadcase)%logscale == 1_pInt) then ! loglinear scale + if (loadcase == 1_pInt) then ! 1st loadcase of loglinear scale + if (inc == 1_pInt) then ! 1st inc of 1st loadcase of loglinear scale + timeinc = bc(1)%time*(2.0_pReal**real( 1_pInt-bc(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd + else ! not-1st inc of 1st loadcase of loglinear scale + timeinc = bc(1)%time*(2.0_pReal**real(inc-1_pInt-bc(1)%incs ,pReal)) endif else ! not-1st loadcase of logarithmic scale - timeinc = time0 *( (1.0_pReal + bc(loadcase)%timeIncrement/time0 )**real( step/bc(loadcase)%steps ,pReal) & - -(1.0_pReal + bc(loadcase)%timeIncrement/time0 )**real( (step-1_pInt)/bc(loadcase)%steps ,pReal) ) + timeinc = time0 *( (1.0_pReal + bc(loadcase)%time/time0 )**real( inc/bc(loadcase)%incs ,pReal) & + -(1.0_pReal + bc(loadcase)%time/time0 )**real( (inc-1_pInt)/bc(loadcase)%incs ,pReal) ) endif endif time = time + timeinc - totalStepsCounter = totalStepsCounter + 1_pInt + totalIncsCounter = totalIncsCounter + 1_pInt !************************************************************* ! Initialization Start !************************************************************* - if(totalStepsCounter == restartReadStep) then ! Initialize values - if (regrid .eqv. .true. ) then ! 'DeInitialize' the values changing in case of regridding - regrid = .false. - - call dfftw_destroy_plan(fftw_plan(1)); call dfftw_destroy_plan(fftw_plan(2)) - if(debugDivergence) call dfftw_destroy_plan(fftw_plan(3)) - deallocate (defgrad) - deallocate (defgradold) - deallocate (coordinates) - deallocate (temperature) - deallocate (xi) - deallocate (workfft) - !ToDo: here we have to create the new geometry and assign the values from the previous step + if(totalIncsCounter == restartReadInc) then ! Initialize values + guessmode = 0.0_pReal ! no old values + allocate (defgrad ( res(1), res(2),res(3),3,3)); defgrad = 0.0_pReal + allocate (defgradold ( res(1), res(2),res(3),3,3)); defgradold = 0.0_pReal + allocate (coordinates(3,res(1), res(2),res(3))); coordinates = 0.0_pReal + allocate (temperature( res(1), res(2),res(3))); temperature = bc(1)%temperature ! start out isothermally + allocate (xi (3,res1_red,res(2),res(3))); xi = 0.0_pReal + data_fftw = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8) + call c_f_pointer(data_fftw, data_real, [ res(1)+2_pInt,res(2),res(3),3,3]) + call c_f_pointer(data_fftw, data_complex, [ res1_red, res(2),res(3),3,3]) + + if (debugDivergence) then + divergence = fftw_alloc_complex(int(res1_red*res(2)*res(3)*3_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8) + call c_f_pointer(divergence, divergence_real, [ res(1)+2_pInt,res(2),res(3),3]) + call c_f_pointer(divergence, divergence_complex, [ res1_red, res(2),res(3),3]) + endif + if (debugFFTW) then + fftw_debug = fftw_alloc_complex(int(res1_red*res(2)*res(3),C_SIZE_T)) !C_SIZE_T is of type integer(8) + call c_f_pointer(fftw_debug, fftw_debug_real, [ res(1)+2_pInt,res(2),res(3)]) + call c_f_pointer(fftw_debug, fftw_debug_complex, [ res1_red, res(2),res(3)]) endif - guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first step - allocate (defgrad ( res(1),res(2),res(3),3,3)); defgrad = 0.0_pReal - allocate (defgradold ( res(1),res(2),res(3),3,3)); defgradold = 0.0_pReal - allocate (coordinates(3,res(1),res(2),res(3))); coordinates = 0.0_pReal - allocate (temperature( res(1),res(2),res(3))); temperature = bc(1)%temperature ! start out isothermally - allocate (xi (3,res(1)/2+1,res(2),res(3))); xi =0.0_pReal - allocate (workfft(res(1)+2,res(2),res(3),3,3)); workfft = 0.0_pReal - if (debugDivergence) allocate (divergence(res(1)+2,res(2),res(3),3)); divergence = 0.0_pReal - wgt = 1.0_pReal/real(Npoints, pReal) - - call dfftw_plan_many_dft_r2c(fftw_plan(1),3,(/res(1),res(2),res(3)/),9,& - workfft,(/res(1) +2_pInt,res(2),res(3)/),1,(res(1) +2_pInt)*res(2)*res(3),& - workfft,(/res(1)/2_pInt+1_pInt,res(2),res(3)/),1,(res(1)/2_pInt+1_pInt)*res(2)*res(3),fftw_planner_flag) - call dfftw_plan_many_dft_c2r(fftw_plan(2),3,(/res(1),res(2),res(3)/),9,& - workfft,(/res(1)/2_pInt+1_pInt,res(2),res(3)/),1,(res(1)/2_pInt+1_pInt)*res(2)*res(3),& - workfft,(/res(1) +2_pInt,res(2),res(3)/),1,(res(1) +2_pInt)*res(2)*res(3),fftw_planner_flag) + fftw_stress = fftw_plan_many_dft_r2c(3,(/res(3),res(2) ,res(1)/),9,& ! dimensions , length in each dimension in reversed order + data_real,(/res(3),res(2) ,res(1)+2_pInt/),& ! input data , physical length in each dimension in reversed order + 1, res(3)*res(2)*(res(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions + data_complex,(/res(3),res(2) ,res1_red/),& + 1, res(3)*res(2)* res1_red,fftw_planner_flag) + + fftw_fluctuation = fftw_plan_many_dft_c2r(3,(/res(3),res(2) ,res(1)/),9,& + data_complex,(/res(3),res(2) ,res1_red/),& + 1, res(3)*res(2)* res1_red,& + data_real,(/res(3),res(2) ,res(1)+2_pInt/),& + 1, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) if (debugDivergence) & - call dfftw_plan_many_dft_c2r(fftw_plan(3),3,(/res(1),res(2),res(3)/),3,& - divergence,(/res(1)/2_pInt+1_pInt,res(2),res(3)/),1,(res(1)/2_pInt+1_pInt)*res(2)*res(3),& - divergence,(/res(1) +2_pInt,res(2),res(3)/),1,(res(1) +2_pInt)*res(2)*res(3),fftw_planner_flag) + fftw_divergence = fftw_plan_many_dft_c2r(3,(/res(3),res(2) ,res(1)/),3,& + divergence_complex,(/res(3),res(2) ,res1_red/),& + 1, res(3)*res(2)* res1_red,& + divergence_real,(/res(3),res(2) ,res(1)+2_pInt/),& + 1, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) + + if (debugFFTW) fftw_debug_forward = fftw_plan_dft_r2c_3d(res(3),res(2),res(1),fftw_debug_real,fftw_debug_complex,fftw_planner_flag) !reversed order + if (debugFFTW) fftw_debug_backward= fftw_plan_dft_c2r_3d(res(3),res(2),res(1),fftw_debug_complex,fftw_debug_real,fftw_planner_flag) !reversed order + if (debugGeneral) then write (6,*) 'FFTW initialized' endif - if (restartReadStep==1_pInt) then ! not restarting, no deformation at the beginning + if (restartReadInc==1_pInt) then ! not restarting, no deformation at the beginning do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) defgrad(i,j,k,1:3,1:3) = math_I3 defgradold(i,j,k,1:3,1:3) = math_I3 @@ -533,7 +547,7 @@ program DAMASK_spectral do j = 1_pInt, res(2) k_s(2) = j - 1_pInt if(j > res(2)/2_pInt + 1_pInt) k_s(2) = k_s(2) - res(2) - do i = 1, res(1)/2_pInt + 1_pInt + do i = 1, res1_red k_s(1) = i - 1_pInt xi(3,i,j,k) = 0.0_pReal ! 2D case if(res(3) > 1_pInt) xi(3,i,j,k) = real(k_s(3), pReal)/geomdimension(3) ! 3D case @@ -543,17 +557,17 @@ program DAMASK_spectral !remove the given highest frequencies for calculation of the gamma operator cutting_freq = (/0_pInt,0_pInt,0_pInt/) ! for 0,0,0, just the highest freq. is removed - do k = 1_pInt ,res(3); do j = 1_pInt ,res(2); do i = 1_pInt,res(1)/2_pInt + 1_pInt + do k = 1_pInt ,res(3); do j = 1_pInt ,res(2); do i = 1_pInt, res1_red if((k .gt. res(3)/2_pInt - cutting_freq(3)).and.(k .le. res(3)/2_pInt + 1_pInt + cutting_freq(3))) xi(3,i,j,k)= 0.0_pReal if((j .gt. res(2)/2_pInt - cutting_freq(2)).and.(j .le. res(2)/2_pInt + 1_pInt + cutting_freq(2))) xi(2,i,j,k)= 0.0_pReal - if((i .gt. res(1)/2_pInt - cutting_freq(1)).and.(i .le. res(2)/2_pInt + 1_pInt + cutting_freq(1))) xi(1,i,j,k)= 0.0_pReal + if((i .gt. res(1)/2_pInt - cutting_freq(1)).and.(i .le. res(1)/2_pInt + 1_pInt + cutting_freq(1))) xi(1,i,j,k)= 0.0_pReal enddo; enddo; enddo if(memory_efficient) then ! allocate just single fourth order tensor allocate (gamma_hat(1,1,1,3,3,3,3)); gamma_hat = 0.0_pReal else ! precalculation of gamma_hat field - allocate (gamma_hat(res(1)/2_pInt + 1_pInt ,res(2),res(3),3,3,3,3)); gamma_hat = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)/2_pInt + 1_pInt + allocate (gamma_hat(res1_red ,res(2),res(3),3,3,3,3)); gamma_hat = 0.0_pReal + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red if (any(xi(1:3,i,j,k) /= 0.0_pReal)) then do l = 1_pInt ,3_pInt; do m = 1_pInt,3_pInt xiDyad(l,m) = xi(l,i,j,k)*xi(m,i,j,k) @@ -592,12 +606,13 @@ program DAMASK_spectral write(538), 'materialpoint_sizeResults', materialpoint_sizeResults write(538), 'loadcases', N_Loadcases write(538), 'frequencies', bc(1:N_Loadcases)%outputfrequency ! one entry per loadcase - write(538), 'times', bc(1:N_Loadcases)%timeIncrement ! one entry per loadcase - write(538), 'logscales', bc(1:N_Loadcases)%logscale ! one entry per loadcase (0: linear, 1: log) - bc(1)%steps= bc(1)%steps + 1_pInt - write(538), 'increments', bc(1:N_Loadcases)%steps ! one entry per loadcase - bc(1)%steps= bc(1)%steps - 1_pInt - write(538), 'startingIncrement', restartReadStep -1_pInt ! start with writing out the previous step + write(538), 'times', bc(1:N_Loadcases)%time ! one entry per loadcase + write(538), 'logscales', bc(1:N_Loadcases)%logscale + bc(1)%incs = bc(1)%incs + 1_pInt + write(538), 'increments', bc(1:N_Loadcases)%incs ! one entry per loadcase + bc(1)%incs = bc(1)%incs - 1_pInt + write(538), 'startingIncrement', restartReadInc -1_pInt ! start with writing out the previous inc + write(538), 'eoh' ! end of header write(538), materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! initial (non-deformed) results !$OMP END CRITICAL (write2out) @@ -606,9 +621,9 @@ program DAMASK_spectral ! Initialization End !************************************************************* - if(totalStepsCounter >= restartReadStep) then ! Do calculations (otherwise just forwarding) + if(totalIncsCounter >= restartReadInc) then ! Do calculations (otherwise just forwarding) if(bc(loadcase)%restartFrequency>0_pInt) & - restartWrite = ( mod(step - 1_pInt,bc(loadcase)%restartFrequency)==0_pInt) ! at frequency of writing restart information + restartWrite = ( mod(inc - 1_pInt,bc(loadcase)%restartFrequency)==0_pInt) ! at frequency of writing restart information ! setting restart parameter for FEsolving (first call to CPFEM_general will write ToDo: true?) if (bc(loadcase)%velGradApplied) & ! calculate fDot from given L and current F fDot = math_mul33x33(bc(loadcase)%deformation, defgradAim) @@ -650,7 +665,7 @@ program DAMASK_spectral iter = 0_pInt err_div = 2.0_pReal * err_div_tol ! go into loop - c_prev = math_rotate_forward3x3x3x3(c_current*wgt,bc(loadcase)%rotation) ! calculate stiffness from former step + ! c_prev = math_rotate_forward3x3x3x3(c_current*wgt,bc(loadcase)%rotation) ToDo: ask Philip ! calculate stiffness from former inc if(size_reduced > 0_pInt) then ! calculate compliance in case stress BC is applied c_prev99 = math_Plain3333to99(c_prev) k = 0_pInt ! build reduced stiffness @@ -681,9 +696,9 @@ program DAMASK_spectral print '(a)', '#############################################################' - print '(A,I5.5,A,es12.6)', 'Increment ', totalStepsCounter, ' Time ',time + print '(A,I5.5,A,es12.6)', 'Increment ', totalIncsCounter, ' Time ',time if (restartWrite ) then - print '(A)', 'writing converged results of previous step for restart' + print '(A)', 'writing converged results of previous inc for restart' if(IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(defgrad))) then ! and writing deformation gradient field to file write (777,rec=1) defgrad close (777) @@ -700,7 +715,7 @@ program DAMASK_spectral print '(a)', '' print '(a)', '=============================================================' - print '(5(a,i5.5))', 'Loadcase ',loadcase,' Step ',step,'/',bc(loadcase)%steps,'@Iteration ',iter,'/',itmax + print '(5(a,i5.5))', 'Loadcase ',loadcase,' Increment ',inc,'/',bc(loadcase)%incs,'@Iteration ',iter,'/',itmax do n = 1_pInt,3_pInt; do m = 1_pInt,3_pInt defgrad_av_lab(m,n) = sum(defgrad(1:res(1),1:res(2),1:res(3),m,n)) * wgt enddo; enddo @@ -718,9 +733,14 @@ program DAMASK_spectral cstress,dsde, pstress, dPdF) enddo; enddo; enddo - workfft = 0.0_pReal ! needed because of the padding for FFTW + data_real = 0.0_pReal ! needed because of the padding for FFTW c_current = 0.0_pReal - ielem = 0_pInt + ielem = 0_pInt + if (debugFFTW) then + row = (mod(totalIncsCounter+iter+7_pInt,9_pInt))/3_pInt + 1_pInt ! go through the elements of the tensors, controlled by totalIncsCounter and iter, starting at 1 + column = (mod(totalIncsCounter+iter+1_pInt,3_pInt)) + 1_pInt + endif + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) ielem = ielem + 1_pInt call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1, @@ -729,41 +749,47 @@ program DAMASK_spectral temperature(i,j,k),timeinc,ielem,1_pInt,& cstress,dsde, pstress,dPdF) CPFEM_mode = 2_pInt - - workfft(i,j,k,1:3,1:3) = pstress ! build up average P-K stress + + data_real(i,j,k,1:3,1:3) = pstress + if (debugFFTW) fftw_debug_real(i,j,k) = pstress(row,column) ! choose an arbitrary component c_current = c_current + dPdF enddo; enddo; enddo do n = 1_pInt,3_pInt; do m = 1_pInt,3_pInt - pstress_av_lab(m,n) = sum(workfft(1:res(1),1:res(2),1:res(3),m,n)) * wgt + pstress_av_lab(m,n) = sum(data_real(1:res(1),1:res(2),1:res(3),m,n)) * wgt enddo; enddo print '(a)', '' print '(a)', '... calculating equilibrium with spectral method ............' - call dfftw_execute_dft_r2c(fftw_plan(1),workfft,workfft) ! FFT of pstress - - p_hat_avg = sqrt(maxval (math_eigenvalues3x3(math_mul33x33(workfft(1,1,1,1:3,1:3),& ! L_2 norm of average stress (freq 0,0,0) in fourier space, - math_transpose3x3(workfft(1,1,1,1:3,1:3)))))) ! ignore imaginary part as it is always zero for real only input + call fftw_execute_dft_r2c(fftw_stress,data_real,data_complex) ! FFT of pstress + if (debugFFTW) then + call fftw_execute_dft_r2c(fftw_debug_forward,fftw_debug_real,fftw_debug_complex) + print '(a,i1,x,i1)', 'checking FT results of compontent ', row, column + print '(a,2(es10.4,x))', 'max FT relative error ',& + maxval( real((fftw_debug_complex(1:res1_red,1:res(2),1:res(3))- data_complex(1:res1_red,1:res(2),1:res(3),row,column))/fftw_debug_real(1:res1_red,1:res(2),1:res(3)))), & + maxval(aimag((fftw_debug_complex(1:res1_red,1:res(2),1:res(3))- data_complex(1:res1_red,1:res(2),1:res(3),row,column))/fftw_debug_real(1:res1_red,1:res(2),1:res(3)))) + fftw_debug_complex = 0.0_pReal + endif + p_hat_avg = sqrt(maxval (math_eigenvalues3x3(math_mul33x33(real(data_complex(1,1,1,1:3,1:3)),& ! L_2 norm of average stress (freq 0,0,0) in fourier space, + math_transpose3x3(real(data_complex(1,1,1,1:3,1:3))))))) ! ignore imaginary part as it is always zero for real only input err_div_avg_complex = 0.0_pReal err_div_max = 0.0_pReal ! only important if debugDivergence == .true. - divergence = 0.0_pReal ! - '' - + divergence_complex = 0.0_pReal ! - '' - - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)/2_pInt+1_pInt - divergence_complex = math_mul33x3_complex(workfft(i*2_pInt-1_pInt,j,k,1:3,1:3)& ! calculate divergence out of the real and imaginary part of the stress - + workfft(i*2_pInt ,j,k,1:3,1:3)*img,xi(1:3,i,j,k)) + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red + temp3_Complex = math_mul33x3_complex(data_complex(i,j,k,1:3,1:3),xi(1:3,i,j,k)) ! calculate divergence without factor of 2*pi*img if(debugDivergence) then ! need divergence NOT squared - divergence(i*2_pInt-1_pInt,j,k,1:3) = -aimag(divergence_complex) ! real part at i*2-1, imaginary part at i*2, multiply by i - divergence(i*2_pInt ,j,k,1:3) = real(divergence_complex) ! ==> switch order and change sign of img part + divergence_complex(i,j,k,1:3) = temp3_Complex *img !ToDo negativ img? endif - divergence_complex = divergence_complex**2.0_pReal ! all criteria need divergence squared - if(i==1_pInt .or. i == res(1)/2_pInt + 1_pInt) then ! We are on one of the two slides without conjg. complex counterpart - err_div_avg_complex = err_div_avg_complex + sum(divergence_complex) ! RMS of L_2 norm of div(stress) in fourier space (Suquet small strain) + temp3_Complex = temp3_Complex**2.0_pReal ! all criteria need divergence squared + if(i==1_pInt .or. i == res1_red) then ! We are on one of the two slides without conjg. complex counterpart + err_div_avg_complex = err_div_avg_complex + sum(temp3_Complex) ! RMS of L_2 norm of div(stress) in fourier space (Suquet small strain) else ! Has somewhere a conj. complex counterpart. Therefore count it twice. - err_div_avg_complex = err_div_avg_complex +2.0*real(sum(divergence_complex)) ! Ignore img part (conjg. complex sum will end up 0). This and the different order + err_div_avg_complex = err_div_avg_complex +2.0*real(sum(temp3_Complex) ) ! Ignore img part (conjg. complex sum will end up 0). This and the different order endif ! compared to c2c transform results in slight numerical deviations. if(debugDivergence) then - err_div_max = max(err_div_max,abs(sqrt(sum(divergence_complex)))) ! maximum of L two norm of div(stress) in fourier space (Suquet large strain) + err_div_max = max(err_div_max,abs(sqrt(sum(temp3_Complex)))) ! maximum of L two norm of div(stress) in fourier space (Suquet large strain) endif enddo; enddo; enddo @@ -773,13 +799,13 @@ program DAMASK_spectral ! calculate additional divergence criteria and report ------------- if(debugDivergence) then - call dfftw_execute_dft_c2r(fftw_plan(3),divergence,divergence) - divergence = divergence *pi*2.0_pReal* wgt !pointwise factor 2*pi from differentation and weighting from FT + call fftw_execute_dft_c2r(fftw_divergence,divergence_complex,divergence_real) + divergence_real = divergence_real *pi*2.0_pReal* wgt !pointwise factor 2*pi from differentation and weighting from FT err_real_div_avg = 0.0_pReal err_real_div_max = 0.0_pReal do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - err_real_div_avg = err_real_div_avg + sum(divergence(i,j,k,1:3)**2.0_pReal) ! avg of L_2 norm of div(stress) in real space - err_real_div_max = max(err_real_div_max, sqrt(sum(divergence(i,j,k,1:3)**2.0_pReal))) ! maximum of L two norm of div(stress) in real space + err_real_div_avg = err_real_div_avg + sum(divergence_real(i,j,k,1:3)**2.0_pReal) ! avg of L_2 norm of div(stress) in real space + err_real_div_max = max(err_real_div_max, sqrt(sum(divergence_real(i,j,k,1:3)**2.0_pReal))) ! maximum of L two norm of div(stress) in real space enddo; enddo; enddo p_real_avg = sqrt(maxval (math_eigenvalues3x3(math_mul33x33(pstress_av_lab,& ! L_2 norm of average stress in real space, math_transpose3x3(pstress_av_lab))))) @@ -795,7 +821,7 @@ program DAMASK_spectral endif ! -------------------------- if(memory_efficient) then ! memory saving version, on-the-fly calculation of gamma_hat - do k = 1_pInt, res(3); do j = 1_pInt, res(2) ;do i = 1_pInt, res(1)/2_pInt+1_pInt + do k = 1_pInt, res(3); do j = 1_pInt, res(2) ;do i = 1_pInt, res1_red if (any(xi(1:3,i,j,k) /= 0.0_pReal)) then do l = 1_pInt,3_pInt; do m = 1_pInt,3_pInt xiDyad(l,m) = xi(l,i,j,k)*xi(m,i,j,k) @@ -810,29 +836,33 @@ program DAMASK_spectral (xiDyad(m,p) +xiDyad(p,m)) enddo; enddo; enddo; enddo do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt - temp33_Complex(m,n) = sum(gamma_hat(1,1,1,m,n,1:3,1:3) *(workfft(i*2_pInt-1_pInt,j,k,1:3,1:3)& - +workfft(i*2_pInt ,j,k,1:3,1:3)*img)) + temp33_Complex(m,n) = sum(gamma_hat(1,1,1, m,n, 1:3,1:3) * data_complex(i,j,k,1:3,1:3)) enddo; enddo - workfft(i*2_pInt-1_pInt,j,k,1:3,1:3) = real (temp33_Complex) - workfft(i*2_pInt ,j,k,1:3,1:3) = aimag(temp33_Complex) + data_complex(i,j,k,1:3,1:3) = temp33_Complex + if (debugFFTW) fftw_debug_complex(i,j,k) = temp33_Complex(row,column) enddo; enddo; enddo else ! use precalculated gamma-operator - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)/2_pInt+1_pInt + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt - temp33_Complex(m,n) = sum(gamma_hat(i,j,k, m,n,1:3,1:3) *(workfft(i*2_pInt-1_pInt,j,k,1:3,1:3)& - + workfft(i*2_pInt ,j,k,1:3,1:3)*img)) + temp33_Complex(m,n) = sum(gamma_hat(i,j,k, m,n, 1:3,1:3) * data_complex(i,j,k,1:3,1:3)) enddo; enddo - workfft(i*2_pInt-1_pInt,j,k,1:3,1:3) = real (temp33_Complex) - workfft(i*2_pInt ,j,k,1:3,1:3) = aimag(temp33_Complex) + data_complex(i,j,k,1:3,1:3) = temp33_Complex + if (debugFFTW) fftw_debug_complex(i,j,k) = temp33_Complex(row,column) enddo; enddo; enddo endif - workfft(1,1,1,1:3,1:3) = defgrad_av_lab - math_I3 ! assign zero frequency (real part) with average displacement gradient - workfft(2,1,1,1:3,1:3) = 0.0_pReal ! zero frequency (imaginary part) + data_complex(1,1,1,1:3,1:3) = defgrad_av_lab - math_I3 ! assign zero frequency (real part) with average displacement gradient + if (debugFFTW) fftw_debug_complex(1,1,1) = data_complex(1,1,1,row,column) - call dfftw_execute_dft_c2r(fftw_plan(2),workfft,workfft) ! back transform of fluct deformation gradient + call fftw_execute_dft_c2r(fftw_fluctuation,data_complex,data_real) ! back transform of fluct deformation gradient + if (debugFFTW) then + print '(a,i1,x,i1)', 'checking iFT results of compontent ', row, column + call fftw_execute_dft_c2r(fftw_debug_backward,fftw_debug_complex,fftw_debug_real) + print '(a,es10.4)', 'max iFT relative error ',& + maxval((fftw_debug_real(1:res(1),1:res(2),1:res(3))- data_real(1:res(1),1:res(2),1:res(3),row,column))/fftw_debug_real(1:res(1),1:res(2),1:res(3))) + endif - defgrad = defgrad + workfft(1:res(1),1:res(2),1:res(3),1:3,1:3)*wgt ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization + defgrad = defgrad + data_real(1:res(1),1:res(2),1:res(3),1:3,1:3)*wgt ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt defgrad_av_lab(m,n) = sum(defgrad(1:res(1),1:res(2),1:res(3),m,n))*wgt ! ToDo: check whether this needs recalculation or is equivalent to former defgrad_av @@ -885,12 +915,12 @@ program DAMASK_spectral print '(a)', '' print '(a)', '=============================================================' if(err_div > err_div_tol .or. err_stress > err_stress_tol) then - print '(A,I5.5,A)', 'increment ', totalStepsCounter, ' NOT converged' + print '(A,I5.5,A)', 'increment ', totalIncsCounter, ' NOT converged' notConvergedCounter = notConvergedCounter + 1_pInt else - print '(A,I5.5,A)', 'increment ', totalStepsCounter, ' converged' + print '(A,I5.5,A)', 'increment ', totalIncsCounter, ' converged' endif - if (mod(totalStepsCounter -1_pInt,bc(loadcase)%outputfrequency) == 0_pInt) then ! at output frequency + if (mod(totalIncsCounter -1_pInt,bc(loadcase)%outputfrequency) == 0_pInt) then ! at output frequency print '(a)', '' print '(a)', '... writing results to file .................................' write(538), materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! write result to file @@ -938,19 +968,23 @@ program DAMASK_spectral ! # end test of regridding endif ! end calculation/forwarding - enddo ! end looping over steps in current loadcase + enddo ! end looping over incs in current loadcase deallocate(c_reduced) deallocate(s_reduced) enddo ! end looping over loadcases !$OMP CRITICAL (write2out) print '(a)', '' print '(a)', '#############################################################' - print '(a,i5.5,a,i5.5,a)', 'of ', totalStepsCounter - restartReadStep + 1_pInt, ' calculated steps, ',& - notConvergedCounter, ' steps did not converge!' + print '(a,i5.5,a,i5.5,a)', 'of ', totalIncsCounter - restartReadInc + 1_pInt, ' calculated increments, ',& + notConvergedCounter, ' increments did not converge!' !$OMP END CRITICAL (write2out) close(538) - call dfftw_destroy_plan(fftw_plan(1)); call dfftw_destroy_plan(fftw_plan(2)) - if (debugDivergence) call dfftw_destroy_plan(fftw_plan(3)) + call fftw_destroy_plan(fftw_stress); call fftw_destroy_plan(fftw_fluctuation) + if (debugDivergence) call fftw_destroy_plan(fftw_divergence) + if (debugFFTW) then + call fftw_destroy_plan(fftw_debug_forward) + call fftw_destroy_plan(fftw_debug_backward) + endif end program DAMASK_spectral !******************************************************************** diff --git a/code/FEsolving.f90 b/code/FEsolving.f90 index ea0fb484b..7d1d744e5 100644 --- a/code/FEsolving.f90 +++ b/code/FEsolving.f90 @@ -25,7 +25,7 @@ use prec, only: pInt,pReal implicit none - integer(pInt) :: cycleCounter = 0_pInt, theInc = -1_pInt, restartReadStep = 1_pInt + integer(pInt) :: cycleCounter = 0_pInt, theInc = -1_pInt, restartReadInc = 1_pInt real(pReal) :: theTime = 0.0_pReal, theDelta = 0.0_pReal logical :: lastIncConverged = .false.,outdatedByNewInc = .false.,outdatedFFN1 = .false.,terminallyIll = .false. logical :: symmetricSolver = .false. @@ -74,9 +74,9 @@ if(start /= 0_pInt) then ! found something length = verify(commandLine(start:len(commandLine)),'0123456789',.false.) ! where is first non number after argument? - read(commandLine(start:start+length),'(I12)') restartReadStep ! read argument - restartRead = max(0_pInt,restartReadStep) > 0_pInt - if(restartReadStep < 0_pInt) call IO_warning(warning_ID=34_pInt) + read(commandLine(start:start+length),'(I12)') restartReadInc ! read argument + restartRead = max(0_pInt,restartReadInc) > 0_pInt + if(restartReadInc < 0_pInt) call IO_warning(warning_ID=34_pInt) endif else rewind(fileunit) diff --git a/code/IO.f90 b/code/IO.f90 index ee6e371c9..022e222c8 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -1172,6 +1172,8 @@ endfunction msg = 'opening material configuration' case (101) msg = 'opening input file' + case (102) + msg = 'precistion not suitable for FFTW' case (103) msg = 'odd resolution given' case (104) @@ -1412,8 +1414,6 @@ endfunction character(len=1024) msg select case (warning_ID) - case (33_pInt) - msg = 'cannot guess along trajectory for first step of first loadcase' case (34) msg = 'invalid restart increment given' case (47_pInt) diff --git a/code/core_modules.f90 b/code/core_modules.f90 index 9389115ce..74cdf6c0d 100644 --- a/code/core_modules.f90 +++ b/code/core_modules.f90 @@ -32,7 +32,7 @@ END MODULE prec MODULE debug use prec, only: pInt implicit none - integer(pInt), parameter :: debug_verbosity = 0_pInt + integer(pInt), parameter :: debug_verbosity = 1_pInt END MODULE debug MODULE numerics diff --git a/code/damask.core.pyf b/code/damask.core.pyf index f6affd82a..82b324018 100644 --- a/code/damask.core.pyf +++ b/code/damask.core.pyf @@ -70,36 +70,29 @@ python module core ! in real*8, dimension(res[0], res[1], res[2], 3,3),intent(in), depend(res[0],res[1],res[2]) :: defgrad ! output variables real*8, dimension(res[0], res[1], res[2], 3), intent(out), depend(res[0],res[1],res[2]) :: coords - ! variables with dimension depending on input - complex*16, dimension(res[0]/2+1,res[1],res[2],3), depend(res[0],res[1],res[2]) :: coords_fft - complex*16, dimension(res[0],res[1],res[2],3,3), depend(res[0],res[1],res[2]) :: defgrad_fft end subroutine deformed_fft - subroutine curl_fft(res,geomdim,vec_tens,field,curl_field) ! in :math:math.f90 + subroutine curl_fft(res,geomdim,vec_tens,field,curl) ! in :math:math.f90 ! input variables integer, dimension(3), intent(in) :: res real*8, dimension(3), intent(in) :: geomdim integer, intent(in) :: vec_tens real*8, dimension(res[0], res[1], res[2], 3,vec_tens), intent(in), depend(res[0],res[1],res[2],vec_tens) :: field ! output variables - real*8, dimension(res[0], res[1], res[2], 3,vec_tens), intent(out), depend(res[0],res[1],res[2],vec_tens) :: curl_field + real*8, dimension(res[0], res[1], res[2], 3,vec_tens), intent(out), depend(res[0],res[1],res[2],vec_tens) :: curl ! variables with dimension depending on input - complex*16, dimension(res[0], res[1],res[2],3,vec_tens), depend(res[0],res[1],res[2],vec_tens) :: field_fft - complex*16, dimension(res[0]/2+1,res[1],res[2],3,vec_tens), depend(res[0],res[1],res[2],vec_tens) :: curl_field_fft real*8, dimension(res[0]/2+1,res[1],res[2],3), depend(res[0],res[1],res[2]) :: xi end subroutine curl_fft - subroutine divergence_fft(res,geomdim,vec_tens,field,divergence_field) ! in :math:math.f90 + subroutine divergence_fft(res,geomdim,vec_tens,field,divergence) ! in :math:math.f90 ! input variables integer, dimension(3), intent(in) :: res real*8, dimension(3), intent(in) :: geomdim integer, intent(in) :: vec_tens real*8, dimension(res[0], res[1], res[2], 3,vec_tens), intent(in), depend(res[0],res[1],res[2],vec_tens) :: field ! output variables - real*8, dimension(res[0], res[1], res[2], vec_tens), intent(out), depend(res[0],res[1],res[2],vec_tens) :: divergence_field - ! variables with dimension depending on input - complex*16, dimension(res[0], res[1],res[2],vec_tens,3),i depend(res[0],res[1],res[2],vec_tens) :: field_fft - complex*16, dimension(res[0]/2+1,res[1],res[2],vec_tens), depend(res[0],res[1],res[2],vec_tens) :: divergence_field_fft + real*8, dimension(res[0], res[1], res[2], vec_tens), intent(out), depend(res[0],res[1],res[2],vec_tens) :: divergence + ! variables with dimension depending on input real*8, dimension(res[0]/2+1,res[1],res[2],3), depend(res[0],res[1],res[2],3) :: xi end subroutine divergence_fft diff --git a/code/fftw3.f03 b/code/fftw3.f03 new file mode 100644 index 000000000..933a53d57 --- /dev/null +++ b/code/fftw3.f03 @@ -0,0 +1,1819 @@ +! Generated automatically. DO NOT EDIT! + + integer, parameter :: C_FFTW_R2R_KIND = C_INT32_T + + integer(C_INT), parameter :: FFTW_R2HC = 0 + integer(C_INT), parameter :: FFTW_HC2R = 1 + integer(C_INT), parameter :: FFTW_DHT = 2 + integer(C_INT), parameter :: FFTW_REDFT00 = 3 + integer(C_INT), parameter :: FFTW_REDFT01 = 4 + integer(C_INT), parameter :: FFTW_REDFT10 = 5 + integer(C_INT), parameter :: FFTW_REDFT11 = 6 + integer(C_INT), parameter :: FFTW_RODFT00 = 7 + integer(C_INT), parameter :: FFTW_RODFT01 = 8 + integer(C_INT), parameter :: FFTW_RODFT10 = 9 + integer(C_INT), parameter :: FFTW_RODFT11 = 10 + integer(C_INT), parameter :: FFTW_FORWARD = -1 + integer(C_INT), parameter :: FFTW_BACKWARD = +1 + integer(C_INT), parameter :: FFTW_MEASURE = 0 + integer(C_INT), parameter :: FFTW_DESTROY_INPUT = 1 + integer(C_INT), parameter :: FFTW_UNALIGNED = 2 + integer(C_INT), parameter :: FFTW_CONSERVE_MEMORY = 4 + integer(C_INT), parameter :: FFTW_EXHAUSTIVE = 8 + integer(C_INT), parameter :: FFTW_PRESERVE_INPUT = 16 + integer(C_INT), parameter :: FFTW_PATIENT = 32 + integer(C_INT), parameter :: FFTW_ESTIMATE = 64 + integer(C_INT), parameter :: FFTW_ESTIMATE_PATIENT = 128 + integer(C_INT), parameter :: FFTW_BELIEVE_PCOST = 256 + integer(C_INT), parameter :: FFTW_NO_DFT_R2HC = 512 + integer(C_INT), parameter :: FFTW_NO_NONTHREADED = 1024 + integer(C_INT), parameter :: FFTW_NO_BUFFERING = 2048 + integer(C_INT), parameter :: FFTW_NO_INDIRECT_OP = 4096 + integer(C_INT), parameter :: FFTW_ALLOW_LARGE_GENERIC = 8192 + integer(C_INT), parameter :: FFTW_NO_RANK_SPLITS = 16384 + integer(C_INT), parameter :: FFTW_NO_VRANK_SPLITS = 32768 + integer(C_INT), parameter :: FFTW_NO_VRECURSE = 65536 + integer(C_INT), parameter :: FFTW_NO_SIMD = 131072 + integer(C_INT), parameter :: FFTW_NO_SLOW = 262144 + integer(C_INT), parameter :: FFTW_NO_FIXED_RADIX_LARGE_N = 524288 + integer(C_INT), parameter :: FFTW_ALLOW_PRUNING = 1048576 + integer(C_INT), parameter :: FFTW_WISDOM_ONLY = 2097152 + + type, bind(C) :: fftw_iodim + integer(C_INT) n, is, os + end type fftw_iodim + type, bind(C) :: fftw_iodim64 + integer(C_INTPTR_T) n, is, os + end type fftw_iodim64 + + interface + type(C_PTR) function fftw_plan_dft(rank,n,in,out,sign,flags) bind(C, name='fftw_plan_dft') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftw_plan_dft + + type(C_PTR) function fftw_plan_dft_1d(n,in,out,sign,flags) bind(C, name='fftw_plan_dft_1d') + import + integer(C_INT), value :: n + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftw_plan_dft_1d + + type(C_PTR) function fftw_plan_dft_2d(n0,n1,in,out,sign,flags) bind(C, name='fftw_plan_dft_2d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftw_plan_dft_2d + + type(C_PTR) function fftw_plan_dft_3d(n0,n1,n2,in,out,sign,flags) bind(C, name='fftw_plan_dft_3d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + integer(C_INT), value :: n2 + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftw_plan_dft_3d + + type(C_PTR) function fftw_plan_many_dft(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,sign,flags) & + bind(C, name='fftw_plan_many_dft') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + integer(C_INT), value :: howmany + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + integer(C_INT), dimension(*), intent(in) :: inembed + integer(C_INT), value :: istride + integer(C_INT), value :: idist + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), dimension(*), intent(in) :: onembed + integer(C_INT), value :: ostride + integer(C_INT), value :: odist + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftw_plan_many_dft + + type(C_PTR) function fftw_plan_guru_dft(rank,dims,howmany_rank,howmany_dims,in,out,sign,flags) & + bind(C, name='fftw_plan_guru_dft') + import + integer(C_INT), value :: rank + type(fftw_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftw_iodim), dimension(*), intent(in) :: howmany_dims + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftw_plan_guru_dft + + type(C_PTR) function fftw_plan_guru_split_dft(rank,dims,howmany_rank,howmany_dims,ri,ii,ro,io,flags) & + bind(C, name='fftw_plan_guru_split_dft') + import + integer(C_INT), value :: rank + type(fftw_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftw_iodim), dimension(*), intent(in) :: howmany_dims + real(C_DOUBLE), dimension(*), intent(out) :: ri + real(C_DOUBLE), dimension(*), intent(out) :: ii + real(C_DOUBLE), dimension(*), intent(out) :: ro + real(C_DOUBLE), dimension(*), intent(out) :: io + integer(C_INT), value :: flags + end function fftw_plan_guru_split_dft + + type(C_PTR) function fftw_plan_guru64_dft(rank,dims,howmany_rank,howmany_dims,in,out,sign,flags) & + bind(C, name='fftw_plan_guru64_dft') + import + integer(C_INT), value :: rank + type(fftw_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftw_iodim64), dimension(*), intent(in) :: howmany_dims + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftw_plan_guru64_dft + + type(C_PTR) function fftw_plan_guru64_split_dft(rank,dims,howmany_rank,howmany_dims,ri,ii,ro,io,flags) & + bind(C, name='fftw_plan_guru64_split_dft') + import + integer(C_INT), value :: rank + type(fftw_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftw_iodim64), dimension(*), intent(in) :: howmany_dims + real(C_DOUBLE), dimension(*), intent(out) :: ri + real(C_DOUBLE), dimension(*), intent(out) :: ii + real(C_DOUBLE), dimension(*), intent(out) :: ro + real(C_DOUBLE), dimension(*), intent(out) :: io + integer(C_INT), value :: flags + end function fftw_plan_guru64_split_dft + + subroutine fftw_execute_dft(p,in,out) bind(C, name='fftw_execute_dft') + import + type(C_PTR), value :: p + complex(C_DOUBLE_COMPLEX), dimension(*), intent(inout) :: in + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + end subroutine fftw_execute_dft + + subroutine fftw_execute_split_dft(p,ri,ii,ro,io) bind(C, name='fftw_execute_split_dft') + import + type(C_PTR), value :: p + real(C_DOUBLE), dimension(*), intent(inout) :: ri + real(C_DOUBLE), dimension(*), intent(inout) :: ii + real(C_DOUBLE), dimension(*), intent(out) :: ro + real(C_DOUBLE), dimension(*), intent(out) :: io + end subroutine fftw_execute_split_dft + + type(C_PTR) function fftw_plan_many_dft_r2c(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,flags) & + bind(C, name='fftw_plan_many_dft_r2c') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + integer(C_INT), value :: howmany + real(C_DOUBLE), dimension(*), intent(out) :: in + integer(C_INT), dimension(*), intent(in) :: inembed + integer(C_INT), value :: istride + integer(C_INT), value :: idist + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), dimension(*), intent(in) :: onembed + integer(C_INT), value :: ostride + integer(C_INT), value :: odist + integer(C_INT), value :: flags + end function fftw_plan_many_dft_r2c + + type(C_PTR) function fftw_plan_dft_r2c(rank,n,in,out,flags) bind(C, name='fftw_plan_dft_r2c') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + real(C_DOUBLE), dimension(*), intent(out) :: in + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftw_plan_dft_r2c + + type(C_PTR) function fftw_plan_dft_r2c_1d(n,in,out,flags) bind(C, name='fftw_plan_dft_r2c_1d') + import + integer(C_INT), value :: n + real(C_DOUBLE), dimension(*), intent(out) :: in + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftw_plan_dft_r2c_1d + + type(C_PTR) function fftw_plan_dft_r2c_2d(n0,n1,in,out,flags) bind(C, name='fftw_plan_dft_r2c_2d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + real(C_DOUBLE), dimension(*), intent(out) :: in + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftw_plan_dft_r2c_2d + + type(C_PTR) function fftw_plan_dft_r2c_3d(n0,n1,n2,in,out,flags) bind(C, name='fftw_plan_dft_r2c_3d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + integer(C_INT), value :: n2 + real(C_DOUBLE), dimension(*), intent(out) :: in + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftw_plan_dft_r2c_3d + + type(C_PTR) function fftw_plan_many_dft_c2r(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,flags) & + bind(C, name='fftw_plan_many_dft_c2r') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + integer(C_INT), value :: howmany + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + integer(C_INT), dimension(*), intent(in) :: inembed + integer(C_INT), value :: istride + integer(C_INT), value :: idist + real(C_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), dimension(*), intent(in) :: onembed + integer(C_INT), value :: ostride + integer(C_INT), value :: odist + integer(C_INT), value :: flags + end function fftw_plan_many_dft_c2r + + type(C_PTR) function fftw_plan_dft_c2r(rank,n,in,out,flags) bind(C, name='fftw_plan_dft_c2r') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + real(C_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftw_plan_dft_c2r + + type(C_PTR) function fftw_plan_dft_c2r_1d(n,in,out,flags) bind(C, name='fftw_plan_dft_c2r_1d') + import + integer(C_INT), value :: n + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + real(C_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftw_plan_dft_c2r_1d + + type(C_PTR) function fftw_plan_dft_c2r_2d(n0,n1,in,out,flags) bind(C, name='fftw_plan_dft_c2r_2d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + real(C_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftw_plan_dft_c2r_2d + + type(C_PTR) function fftw_plan_dft_c2r_3d(n0,n1,n2,in,out,flags) bind(C, name='fftw_plan_dft_c2r_3d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + integer(C_INT), value :: n2 + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + real(C_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftw_plan_dft_c2r_3d + + type(C_PTR) function fftw_plan_guru_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,out,flags) & + bind(C, name='fftw_plan_guru_dft_r2c') + import + integer(C_INT), value :: rank + type(fftw_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftw_iodim), dimension(*), intent(in) :: howmany_dims + real(C_DOUBLE), dimension(*), intent(out) :: in + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftw_plan_guru_dft_r2c + + type(C_PTR) function fftw_plan_guru_dft_c2r(rank,dims,howmany_rank,howmany_dims,in,out,flags) & + bind(C, name='fftw_plan_guru_dft_c2r') + import + integer(C_INT), value :: rank + type(fftw_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftw_iodim), dimension(*), intent(in) :: howmany_dims + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + real(C_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftw_plan_guru_dft_c2r + + type(C_PTR) function fftw_plan_guru_split_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,ro,io,flags) & + bind(C, name='fftw_plan_guru_split_dft_r2c') + import + integer(C_INT), value :: rank + type(fftw_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftw_iodim), dimension(*), intent(in) :: howmany_dims + real(C_DOUBLE), dimension(*), intent(out) :: in + real(C_DOUBLE), dimension(*), intent(out) :: ro + real(C_DOUBLE), dimension(*), intent(out) :: io + integer(C_INT), value :: flags + end function fftw_plan_guru_split_dft_r2c + + type(C_PTR) function fftw_plan_guru_split_dft_c2r(rank,dims,howmany_rank,howmany_dims,ri,ii,out,flags) & + bind(C, name='fftw_plan_guru_split_dft_c2r') + import + integer(C_INT), value :: rank + type(fftw_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftw_iodim), dimension(*), intent(in) :: howmany_dims + real(C_DOUBLE), dimension(*), intent(out) :: ri + real(C_DOUBLE), dimension(*), intent(out) :: ii + real(C_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftw_plan_guru_split_dft_c2r + + type(C_PTR) function fftw_plan_guru64_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,out,flags) & + bind(C, name='fftw_plan_guru64_dft_r2c') + import + integer(C_INT), value :: rank + type(fftw_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftw_iodim64), dimension(*), intent(in) :: howmany_dims + real(C_DOUBLE), dimension(*), intent(out) :: in + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftw_plan_guru64_dft_r2c + + type(C_PTR) function fftw_plan_guru64_dft_c2r(rank,dims,howmany_rank,howmany_dims,in,out,flags) & + bind(C, name='fftw_plan_guru64_dft_c2r') + import + integer(C_INT), value :: rank + type(fftw_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftw_iodim64), dimension(*), intent(in) :: howmany_dims + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + real(C_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftw_plan_guru64_dft_c2r + + type(C_PTR) function fftw_plan_guru64_split_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,ro,io,flags) & + bind(C, name='fftw_plan_guru64_split_dft_r2c') + import + integer(C_INT), value :: rank + type(fftw_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftw_iodim64), dimension(*), intent(in) :: howmany_dims + real(C_DOUBLE), dimension(*), intent(out) :: in + real(C_DOUBLE), dimension(*), intent(out) :: ro + real(C_DOUBLE), dimension(*), intent(out) :: io + integer(C_INT), value :: flags + end function fftw_plan_guru64_split_dft_r2c + + type(C_PTR) function fftw_plan_guru64_split_dft_c2r(rank,dims,howmany_rank,howmany_dims,ri,ii,out,flags) & + bind(C, name='fftw_plan_guru64_split_dft_c2r') + import + integer(C_INT), value :: rank + type(fftw_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftw_iodim64), dimension(*), intent(in) :: howmany_dims + real(C_DOUBLE), dimension(*), intent(out) :: ri + real(C_DOUBLE), dimension(*), intent(out) :: ii + real(C_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftw_plan_guru64_split_dft_c2r + + subroutine fftw_execute_dft_r2c(p,in,out) bind(C, name='fftw_execute_dft_r2c') + import + type(C_PTR), value :: p + real(C_DOUBLE), dimension(*), intent(inout) :: in + complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + end subroutine fftw_execute_dft_r2c + + subroutine fftw_execute_dft_c2r(p,in,out) bind(C, name='fftw_execute_dft_c2r') + import + type(C_PTR), value :: p + complex(C_DOUBLE_COMPLEX), dimension(*), intent(inout) :: in + real(C_DOUBLE), dimension(*), intent(out) :: out + end subroutine fftw_execute_dft_c2r + + subroutine fftw_execute_split_dft_r2c(p,in,ro,io) bind(C, name='fftw_execute_split_dft_r2c') + import + type(C_PTR), value :: p + real(C_DOUBLE), dimension(*), intent(inout) :: in + real(C_DOUBLE), dimension(*), intent(out) :: ro + real(C_DOUBLE), dimension(*), intent(out) :: io + end subroutine fftw_execute_split_dft_r2c + + subroutine fftw_execute_split_dft_c2r(p,ri,ii,out) bind(C, name='fftw_execute_split_dft_c2r') + import + type(C_PTR), value :: p + real(C_DOUBLE), dimension(*), intent(inout) :: ri + real(C_DOUBLE), dimension(*), intent(inout) :: ii + real(C_DOUBLE), dimension(*), intent(out) :: out + end subroutine fftw_execute_split_dft_c2r + + type(C_PTR) function fftw_plan_many_r2r(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,kind,flags) & + bind(C, name='fftw_plan_many_r2r') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + integer(C_INT), value :: howmany + real(C_DOUBLE), dimension(*), intent(out) :: in + integer(C_INT), dimension(*), intent(in) :: inembed + integer(C_INT), value :: istride + integer(C_INT), value :: idist + real(C_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), dimension(*), intent(in) :: onembed + integer(C_INT), value :: ostride + integer(C_INT), value :: odist + integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind + integer(C_INT), value :: flags + end function fftw_plan_many_r2r + + type(C_PTR) function fftw_plan_r2r(rank,n,in,out,kind,flags) bind(C, name='fftw_plan_r2r') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + real(C_DOUBLE), dimension(*), intent(out) :: in + real(C_DOUBLE), dimension(*), intent(out) :: out + integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind + integer(C_INT), value :: flags + end function fftw_plan_r2r + + type(C_PTR) function fftw_plan_r2r_1d(n,in,out,kind,flags) bind(C, name='fftw_plan_r2r_1d') + import + integer(C_INT), value :: n + real(C_DOUBLE), dimension(*), intent(out) :: in + real(C_DOUBLE), dimension(*), intent(out) :: out + integer(C_FFTW_R2R_KIND), value :: kind + integer(C_INT), value :: flags + end function fftw_plan_r2r_1d + + type(C_PTR) function fftw_plan_r2r_2d(n0,n1,in,out,kind0,kind1,flags) bind(C, name='fftw_plan_r2r_2d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + real(C_DOUBLE), dimension(*), intent(out) :: in + real(C_DOUBLE), dimension(*), intent(out) :: out + integer(C_FFTW_R2R_KIND), value :: kind0 + integer(C_FFTW_R2R_KIND), value :: kind1 + integer(C_INT), value :: flags + end function fftw_plan_r2r_2d + + type(C_PTR) function fftw_plan_r2r_3d(n0,n1,n2,in,out,kind0,kind1,kind2,flags) bind(C, name='fftw_plan_r2r_3d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + integer(C_INT), value :: n2 + real(C_DOUBLE), dimension(*), intent(out) :: in + real(C_DOUBLE), dimension(*), intent(out) :: out + integer(C_FFTW_R2R_KIND), value :: kind0 + integer(C_FFTW_R2R_KIND), value :: kind1 + integer(C_FFTW_R2R_KIND), value :: kind2 + integer(C_INT), value :: flags + end function fftw_plan_r2r_3d + + type(C_PTR) function fftw_plan_guru_r2r(rank,dims,howmany_rank,howmany_dims,in,out,kind,flags) & + bind(C, name='fftw_plan_guru_r2r') + import + integer(C_INT), value :: rank + type(fftw_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftw_iodim), dimension(*), intent(in) :: howmany_dims + real(C_DOUBLE), dimension(*), intent(out) :: in + real(C_DOUBLE), dimension(*), intent(out) :: out + integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind + integer(C_INT), value :: flags + end function fftw_plan_guru_r2r + + type(C_PTR) function fftw_plan_guru64_r2r(rank,dims,howmany_rank,howmany_dims,in,out,kind,flags) & + bind(C, name='fftw_plan_guru64_r2r') + import + integer(C_INT), value :: rank + type(fftw_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftw_iodim64), dimension(*), intent(in) :: howmany_dims + real(C_DOUBLE), dimension(*), intent(out) :: in + real(C_DOUBLE), dimension(*), intent(out) :: out + integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind + integer(C_INT), value :: flags + end function fftw_plan_guru64_r2r + + subroutine fftw_execute_r2r(p,in,out) bind(C, name='fftw_execute_r2r') + import + type(C_PTR), value :: p + real(C_DOUBLE), dimension(*), intent(inout) :: in + real(C_DOUBLE), dimension(*), intent(out) :: out + end subroutine fftw_execute_r2r + + subroutine fftw_destroy_plan(p) bind(C, name='fftw_destroy_plan') + import + type(C_PTR), value :: p + end subroutine fftw_destroy_plan + + subroutine fftw_forget_wisdom() bind(C, name='fftw_forget_wisdom') + import + end subroutine fftw_forget_wisdom + + subroutine fftw_cleanup() bind(C, name='fftw_cleanup') + import + end subroutine fftw_cleanup + + subroutine fftw_set_timelimit(t) bind(C, name='fftw_set_timelimit') + import + real(C_DOUBLE), value :: t + end subroutine fftw_set_timelimit + + subroutine fftw_plan_with_nthreads(nthreads) bind(C, name='fftw_plan_with_nthreads') + import + integer(C_INT), value :: nthreads + end subroutine fftw_plan_with_nthreads + + integer(C_INT) function fftw_init_threads() bind(C, name='fftw_init_threads') + import + end function fftw_init_threads + + subroutine fftw_cleanup_threads() bind(C, name='fftw_cleanup_threads') + import + end subroutine fftw_cleanup_threads + + integer(C_INT) function fftw_export_wisdom_to_filename(filename) bind(C, name='fftw_export_wisdom_to_filename') + import + character(C_CHAR), dimension(*), intent(in) :: filename + end function fftw_export_wisdom_to_filename + + subroutine fftw_export_wisdom_to_file(output_file) bind(C, name='fftw_export_wisdom_to_file') + import + type(C_PTR), value :: output_file + end subroutine fftw_export_wisdom_to_file + + type(C_PTR) function fftw_export_wisdom_to_string() bind(C, name='fftw_export_wisdom_to_string') + import + end function fftw_export_wisdom_to_string + + subroutine fftw_export_wisdom(write_char,data) bind(C, name='fftw_export_wisdom') + import + type(C_FUNPTR), value :: write_char + type(C_PTR), value :: data + end subroutine fftw_export_wisdom + + integer(C_INT) function fftw_import_system_wisdom() bind(C, name='fftw_import_system_wisdom') + import + end function fftw_import_system_wisdom + + integer(C_INT) function fftw_import_wisdom_from_filename(filename) bind(C, name='fftw_import_wisdom_from_filename') + import + character(C_CHAR), dimension(*), intent(in) :: filename + end function fftw_import_wisdom_from_filename + + integer(C_INT) function fftw_import_wisdom_from_file(input_file) bind(C, name='fftw_import_wisdom_from_file') + import + type(C_PTR), value :: input_file + end function fftw_import_wisdom_from_file + + integer(C_INT) function fftw_import_wisdom_from_string(input_string) bind(C, name='fftw_import_wisdom_from_string') + import + character(C_CHAR), dimension(*), intent(in) :: input_string + end function fftw_import_wisdom_from_string + + integer(C_INT) function fftw_import_wisdom(read_char,data) bind(C, name='fftw_import_wisdom') + import + type(C_FUNPTR), value :: read_char + type(C_PTR), value :: data + end function fftw_import_wisdom + + subroutine fftw_fprint_plan(p,output_file) bind(C, name='fftw_fprint_plan') + import + type(C_PTR), value :: p + type(C_PTR), value :: output_file + end subroutine fftw_fprint_plan + + subroutine fftw_print_plan(p) bind(C, name='fftw_print_plan') + import + type(C_PTR), value :: p + end subroutine fftw_print_plan + + type(C_PTR) function fftw_malloc(n) bind(C, name='fftw_malloc') + import + integer(C_SIZE_T), value :: n + end function fftw_malloc + + type(C_PTR) function fftw_alloc_real(n) bind(C, name='fftw_alloc_real') + import + integer(C_SIZE_T), value :: n + end function fftw_alloc_real + + type(C_PTR) function fftw_alloc_complex(n) bind(C, name='fftw_alloc_complex') + import + integer(C_SIZE_T), value :: n + end function fftw_alloc_complex + + subroutine fftw_free(p) bind(C, name='fftw_free') + import + type(C_PTR), value :: p + end subroutine fftw_free + + subroutine fftw_flops(p,add,mul,fmas) bind(C, name='fftw_flops') + import + type(C_PTR), value :: p + real(C_DOUBLE), intent(out) :: add + real(C_DOUBLE), intent(out) :: mul + real(C_DOUBLE), intent(out) :: fmas + end subroutine fftw_flops + + real(C_DOUBLE) function fftw_estimate_cost(p) bind(C, name='fftw_estimate_cost') + import + type(C_PTR), value :: p + end function fftw_estimate_cost + + real(C_DOUBLE) function fftw_cost(p) bind(C, name='fftw_cost') + import + type(C_PTR), value :: p + end function fftw_cost + + end interface + + type, bind(C) :: fftwf_iodim + integer(C_INT) n, is, os + end type fftwf_iodim + type, bind(C) :: fftwf_iodim64 + integer(C_INTPTR_T) n, is, os + end type fftwf_iodim64 + + interface + type(C_PTR) function fftwf_plan_dft(rank,n,in,out,sign,flags) bind(C, name='fftwf_plan_dft') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftwf_plan_dft + + type(C_PTR) function fftwf_plan_dft_1d(n,in,out,sign,flags) bind(C, name='fftwf_plan_dft_1d') + import + integer(C_INT), value :: n + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftwf_plan_dft_1d + + type(C_PTR) function fftwf_plan_dft_2d(n0,n1,in,out,sign,flags) bind(C, name='fftwf_plan_dft_2d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftwf_plan_dft_2d + + type(C_PTR) function fftwf_plan_dft_3d(n0,n1,n2,in,out,sign,flags) bind(C, name='fftwf_plan_dft_3d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + integer(C_INT), value :: n2 + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftwf_plan_dft_3d + + type(C_PTR) function fftwf_plan_many_dft(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,sign,flags) & + bind(C, name='fftwf_plan_many_dft') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + integer(C_INT), value :: howmany + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in + integer(C_INT), dimension(*), intent(in) :: inembed + integer(C_INT), value :: istride + integer(C_INT), value :: idist + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), dimension(*), intent(in) :: onembed + integer(C_INT), value :: ostride + integer(C_INT), value :: odist + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftwf_plan_many_dft + + type(C_PTR) function fftwf_plan_guru_dft(rank,dims,howmany_rank,howmany_dims,in,out,sign,flags) & + bind(C, name='fftwf_plan_guru_dft') + import + integer(C_INT), value :: rank + type(fftwf_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwf_iodim), dimension(*), intent(in) :: howmany_dims + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftwf_plan_guru_dft + + type(C_PTR) function fftwf_plan_guru_split_dft(rank,dims,howmany_rank,howmany_dims,ri,ii,ro,io,flags) & + bind(C, name='fftwf_plan_guru_split_dft') + import + integer(C_INT), value :: rank + type(fftwf_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwf_iodim), dimension(*), intent(in) :: howmany_dims + real(C_FLOAT), dimension(*), intent(out) :: ri + real(C_FLOAT), dimension(*), intent(out) :: ii + real(C_FLOAT), dimension(*), intent(out) :: ro + real(C_FLOAT), dimension(*), intent(out) :: io + integer(C_INT), value :: flags + end function fftwf_plan_guru_split_dft + + type(C_PTR) function fftwf_plan_guru64_dft(rank,dims,howmany_rank,howmany_dims,in,out,sign,flags) & + bind(C, name='fftwf_plan_guru64_dft') + import + integer(C_INT), value :: rank + type(fftwf_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwf_iodim64), dimension(*), intent(in) :: howmany_dims + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftwf_plan_guru64_dft + + type(C_PTR) function fftwf_plan_guru64_split_dft(rank,dims,howmany_rank,howmany_dims,ri,ii,ro,io,flags) & + bind(C, name='fftwf_plan_guru64_split_dft') + import + integer(C_INT), value :: rank + type(fftwf_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwf_iodim64), dimension(*), intent(in) :: howmany_dims + real(C_FLOAT), dimension(*), intent(out) :: ri + real(C_FLOAT), dimension(*), intent(out) :: ii + real(C_FLOAT), dimension(*), intent(out) :: ro + real(C_FLOAT), dimension(*), intent(out) :: io + integer(C_INT), value :: flags + end function fftwf_plan_guru64_split_dft + + subroutine fftwf_execute_dft(p,in,out) bind(C, name='fftwf_execute_dft') + import + type(C_PTR), value :: p + complex(C_FLOAT_COMPLEX), dimension(*), intent(inout) :: in + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out + end subroutine fftwf_execute_dft + + subroutine fftwf_execute_split_dft(p,ri,ii,ro,io) bind(C, name='fftwf_execute_split_dft') + import + type(C_PTR), value :: p + real(C_FLOAT), dimension(*), intent(inout) :: ri + real(C_FLOAT), dimension(*), intent(inout) :: ii + real(C_FLOAT), dimension(*), intent(out) :: ro + real(C_FLOAT), dimension(*), intent(out) :: io + end subroutine fftwf_execute_split_dft + + type(C_PTR) function fftwf_plan_many_dft_r2c(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,flags) & + bind(C, name='fftwf_plan_many_dft_r2c') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + integer(C_INT), value :: howmany + real(C_FLOAT), dimension(*), intent(out) :: in + integer(C_INT), dimension(*), intent(in) :: inembed + integer(C_INT), value :: istride + integer(C_INT), value :: idist + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), dimension(*), intent(in) :: onembed + integer(C_INT), value :: ostride + integer(C_INT), value :: odist + integer(C_INT), value :: flags + end function fftwf_plan_many_dft_r2c + + type(C_PTR) function fftwf_plan_dft_r2c(rank,n,in,out,flags) bind(C, name='fftwf_plan_dft_r2c') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + real(C_FLOAT), dimension(*), intent(out) :: in + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwf_plan_dft_r2c + + type(C_PTR) function fftwf_plan_dft_r2c_1d(n,in,out,flags) bind(C, name='fftwf_plan_dft_r2c_1d') + import + integer(C_INT), value :: n + real(C_FLOAT), dimension(*), intent(out) :: in + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwf_plan_dft_r2c_1d + + type(C_PTR) function fftwf_plan_dft_r2c_2d(n0,n1,in,out,flags) bind(C, name='fftwf_plan_dft_r2c_2d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + real(C_FLOAT), dimension(*), intent(out) :: in + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwf_plan_dft_r2c_2d + + type(C_PTR) function fftwf_plan_dft_r2c_3d(n0,n1,n2,in,out,flags) bind(C, name='fftwf_plan_dft_r2c_3d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + integer(C_INT), value :: n2 + real(C_FLOAT), dimension(*), intent(out) :: in + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwf_plan_dft_r2c_3d + + type(C_PTR) function fftwf_plan_many_dft_c2r(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,flags) & + bind(C, name='fftwf_plan_many_dft_c2r') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + integer(C_INT), value :: howmany + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in + integer(C_INT), dimension(*), intent(in) :: inembed + integer(C_INT), value :: istride + integer(C_INT), value :: idist + real(C_FLOAT), dimension(*), intent(out) :: out + integer(C_INT), dimension(*), intent(in) :: onembed + integer(C_INT), value :: ostride + integer(C_INT), value :: odist + integer(C_INT), value :: flags + end function fftwf_plan_many_dft_c2r + + type(C_PTR) function fftwf_plan_dft_c2r(rank,n,in,out,flags) bind(C, name='fftwf_plan_dft_c2r') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in + real(C_FLOAT), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwf_plan_dft_c2r + + type(C_PTR) function fftwf_plan_dft_c2r_1d(n,in,out,flags) bind(C, name='fftwf_plan_dft_c2r_1d') + import + integer(C_INT), value :: n + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in + real(C_FLOAT), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwf_plan_dft_c2r_1d + + type(C_PTR) function fftwf_plan_dft_c2r_2d(n0,n1,in,out,flags) bind(C, name='fftwf_plan_dft_c2r_2d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in + real(C_FLOAT), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwf_plan_dft_c2r_2d + + type(C_PTR) function fftwf_plan_dft_c2r_3d(n0,n1,n2,in,out,flags) bind(C, name='fftwf_plan_dft_c2r_3d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + integer(C_INT), value :: n2 + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in + real(C_FLOAT), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwf_plan_dft_c2r_3d + + type(C_PTR) function fftwf_plan_guru_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,out,flags) & + bind(C, name='fftwf_plan_guru_dft_r2c') + import + integer(C_INT), value :: rank + type(fftwf_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwf_iodim), dimension(*), intent(in) :: howmany_dims + real(C_FLOAT), dimension(*), intent(out) :: in + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwf_plan_guru_dft_r2c + + type(C_PTR) function fftwf_plan_guru_dft_c2r(rank,dims,howmany_rank,howmany_dims,in,out,flags) & + bind(C, name='fftwf_plan_guru_dft_c2r') + import + integer(C_INT), value :: rank + type(fftwf_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwf_iodim), dimension(*), intent(in) :: howmany_dims + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in + real(C_FLOAT), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwf_plan_guru_dft_c2r + + type(C_PTR) function fftwf_plan_guru_split_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,ro,io,flags) & + bind(C, name='fftwf_plan_guru_split_dft_r2c') + import + integer(C_INT), value :: rank + type(fftwf_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwf_iodim), dimension(*), intent(in) :: howmany_dims + real(C_FLOAT), dimension(*), intent(out) :: in + real(C_FLOAT), dimension(*), intent(out) :: ro + real(C_FLOAT), dimension(*), intent(out) :: io + integer(C_INT), value :: flags + end function fftwf_plan_guru_split_dft_r2c + + type(C_PTR) function fftwf_plan_guru_split_dft_c2r(rank,dims,howmany_rank,howmany_dims,ri,ii,out,flags) & + bind(C, name='fftwf_plan_guru_split_dft_c2r') + import + integer(C_INT), value :: rank + type(fftwf_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwf_iodim), dimension(*), intent(in) :: howmany_dims + real(C_FLOAT), dimension(*), intent(out) :: ri + real(C_FLOAT), dimension(*), intent(out) :: ii + real(C_FLOAT), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwf_plan_guru_split_dft_c2r + + type(C_PTR) function fftwf_plan_guru64_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,out,flags) & + bind(C, name='fftwf_plan_guru64_dft_r2c') + import + integer(C_INT), value :: rank + type(fftwf_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwf_iodim64), dimension(*), intent(in) :: howmany_dims + real(C_FLOAT), dimension(*), intent(out) :: in + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwf_plan_guru64_dft_r2c + + type(C_PTR) function fftwf_plan_guru64_dft_c2r(rank,dims,howmany_rank,howmany_dims,in,out,flags) & + bind(C, name='fftwf_plan_guru64_dft_c2r') + import + integer(C_INT), value :: rank + type(fftwf_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwf_iodim64), dimension(*), intent(in) :: howmany_dims + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: in + real(C_FLOAT), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwf_plan_guru64_dft_c2r + + type(C_PTR) function fftwf_plan_guru64_split_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,ro,io,flags) & + bind(C, name='fftwf_plan_guru64_split_dft_r2c') + import + integer(C_INT), value :: rank + type(fftwf_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwf_iodim64), dimension(*), intent(in) :: howmany_dims + real(C_FLOAT), dimension(*), intent(out) :: in + real(C_FLOAT), dimension(*), intent(out) :: ro + real(C_FLOAT), dimension(*), intent(out) :: io + integer(C_INT), value :: flags + end function fftwf_plan_guru64_split_dft_r2c + + type(C_PTR) function fftwf_plan_guru64_split_dft_c2r(rank,dims,howmany_rank,howmany_dims,ri,ii,out,flags) & + bind(C, name='fftwf_plan_guru64_split_dft_c2r') + import + integer(C_INT), value :: rank + type(fftwf_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwf_iodim64), dimension(*), intent(in) :: howmany_dims + real(C_FLOAT), dimension(*), intent(out) :: ri + real(C_FLOAT), dimension(*), intent(out) :: ii + real(C_FLOAT), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwf_plan_guru64_split_dft_c2r + + subroutine fftwf_execute_dft_r2c(p,in,out) bind(C, name='fftwf_execute_dft_r2c') + import + type(C_PTR), value :: p + real(C_FLOAT), dimension(*), intent(inout) :: in + complex(C_FLOAT_COMPLEX), dimension(*), intent(out) :: out + end subroutine fftwf_execute_dft_r2c + + subroutine fftwf_execute_dft_c2r(p,in,out) bind(C, name='fftwf_execute_dft_c2r') + import + type(C_PTR), value :: p + complex(C_FLOAT_COMPLEX), dimension(*), intent(inout) :: in + real(C_FLOAT), dimension(*), intent(out) :: out + end subroutine fftwf_execute_dft_c2r + + subroutine fftwf_execute_split_dft_r2c(p,in,ro,io) bind(C, name='fftwf_execute_split_dft_r2c') + import + type(C_PTR), value :: p + real(C_FLOAT), dimension(*), intent(inout) :: in + real(C_FLOAT), dimension(*), intent(out) :: ro + real(C_FLOAT), dimension(*), intent(out) :: io + end subroutine fftwf_execute_split_dft_r2c + + subroutine fftwf_execute_split_dft_c2r(p,ri,ii,out) bind(C, name='fftwf_execute_split_dft_c2r') + import + type(C_PTR), value :: p + real(C_FLOAT), dimension(*), intent(inout) :: ri + real(C_FLOAT), dimension(*), intent(inout) :: ii + real(C_FLOAT), dimension(*), intent(out) :: out + end subroutine fftwf_execute_split_dft_c2r + + type(C_PTR) function fftwf_plan_many_r2r(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,kind,flags) & + bind(C, name='fftwf_plan_many_r2r') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + integer(C_INT), value :: howmany + real(C_FLOAT), dimension(*), intent(out) :: in + integer(C_INT), dimension(*), intent(in) :: inembed + integer(C_INT), value :: istride + integer(C_INT), value :: idist + real(C_FLOAT), dimension(*), intent(out) :: out + integer(C_INT), dimension(*), intent(in) :: onembed + integer(C_INT), value :: ostride + integer(C_INT), value :: odist + integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind + integer(C_INT), value :: flags + end function fftwf_plan_many_r2r + + type(C_PTR) function fftwf_plan_r2r(rank,n,in,out,kind,flags) bind(C, name='fftwf_plan_r2r') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + real(C_FLOAT), dimension(*), intent(out) :: in + real(C_FLOAT), dimension(*), intent(out) :: out + integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind + integer(C_INT), value :: flags + end function fftwf_plan_r2r + + type(C_PTR) function fftwf_plan_r2r_1d(n,in,out,kind,flags) bind(C, name='fftwf_plan_r2r_1d') + import + integer(C_INT), value :: n + real(C_FLOAT), dimension(*), intent(out) :: in + real(C_FLOAT), dimension(*), intent(out) :: out + integer(C_FFTW_R2R_KIND), value :: kind + integer(C_INT), value :: flags + end function fftwf_plan_r2r_1d + + type(C_PTR) function fftwf_plan_r2r_2d(n0,n1,in,out,kind0,kind1,flags) bind(C, name='fftwf_plan_r2r_2d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + real(C_FLOAT), dimension(*), intent(out) :: in + real(C_FLOAT), dimension(*), intent(out) :: out + integer(C_FFTW_R2R_KIND), value :: kind0 + integer(C_FFTW_R2R_KIND), value :: kind1 + integer(C_INT), value :: flags + end function fftwf_plan_r2r_2d + + type(C_PTR) function fftwf_plan_r2r_3d(n0,n1,n2,in,out,kind0,kind1,kind2,flags) bind(C, name='fftwf_plan_r2r_3d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + integer(C_INT), value :: n2 + real(C_FLOAT), dimension(*), intent(out) :: in + real(C_FLOAT), dimension(*), intent(out) :: out + integer(C_FFTW_R2R_KIND), value :: kind0 + integer(C_FFTW_R2R_KIND), value :: kind1 + integer(C_FFTW_R2R_KIND), value :: kind2 + integer(C_INT), value :: flags + end function fftwf_plan_r2r_3d + + type(C_PTR) function fftwf_plan_guru_r2r(rank,dims,howmany_rank,howmany_dims,in,out,kind,flags) & + bind(C, name='fftwf_plan_guru_r2r') + import + integer(C_INT), value :: rank + type(fftwf_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwf_iodim), dimension(*), intent(in) :: howmany_dims + real(C_FLOAT), dimension(*), intent(out) :: in + real(C_FLOAT), dimension(*), intent(out) :: out + integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind + integer(C_INT), value :: flags + end function fftwf_plan_guru_r2r + + type(C_PTR) function fftwf_plan_guru64_r2r(rank,dims,howmany_rank,howmany_dims,in,out,kind,flags) & + bind(C, name='fftwf_plan_guru64_r2r') + import + integer(C_INT), value :: rank + type(fftwf_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwf_iodim64), dimension(*), intent(in) :: howmany_dims + real(C_FLOAT), dimension(*), intent(out) :: in + real(C_FLOAT), dimension(*), intent(out) :: out + integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind + integer(C_INT), value :: flags + end function fftwf_plan_guru64_r2r + + subroutine fftwf_execute_r2r(p,in,out) bind(C, name='fftwf_execute_r2r') + import + type(C_PTR), value :: p + real(C_FLOAT), dimension(*), intent(inout) :: in + real(C_FLOAT), dimension(*), intent(out) :: out + end subroutine fftwf_execute_r2r + + subroutine fftwf_destroy_plan(p) bind(C, name='fftwf_destroy_plan') + import + type(C_PTR), value :: p + end subroutine fftwf_destroy_plan + + subroutine fftwf_forget_wisdom() bind(C, name='fftwf_forget_wisdom') + import + end subroutine fftwf_forget_wisdom + + subroutine fftwf_cleanup() bind(C, name='fftwf_cleanup') + import + end subroutine fftwf_cleanup + + subroutine fftwf_set_timelimit(t) bind(C, name='fftwf_set_timelimit') + import + real(C_DOUBLE), value :: t + end subroutine fftwf_set_timelimit + + subroutine fftwf_plan_with_nthreads(nthreads) bind(C, name='fftwf_plan_with_nthreads') + import + integer(C_INT), value :: nthreads + end subroutine fftwf_plan_with_nthreads + + integer(C_INT) function fftwf_init_threads() bind(C, name='fftwf_init_threads') + import + end function fftwf_init_threads + + subroutine fftwf_cleanup_threads() bind(C, name='fftwf_cleanup_threads') + import + end subroutine fftwf_cleanup_threads + + integer(C_INT) function fftwf_export_wisdom_to_filename(filename) bind(C, name='fftwf_export_wisdom_to_filename') + import + character(C_CHAR), dimension(*), intent(in) :: filename + end function fftwf_export_wisdom_to_filename + + subroutine fftwf_export_wisdom_to_file(output_file) bind(C, name='fftwf_export_wisdom_to_file') + import + type(C_PTR), value :: output_file + end subroutine fftwf_export_wisdom_to_file + + type(C_PTR) function fftwf_export_wisdom_to_string() bind(C, name='fftwf_export_wisdom_to_string') + import + end function fftwf_export_wisdom_to_string + + subroutine fftwf_export_wisdom(write_char,data) bind(C, name='fftwf_export_wisdom') + import + type(C_FUNPTR), value :: write_char + type(C_PTR), value :: data + end subroutine fftwf_export_wisdom + + integer(C_INT) function fftwf_import_system_wisdom() bind(C, name='fftwf_import_system_wisdom') + import + end function fftwf_import_system_wisdom + + integer(C_INT) function fftwf_import_wisdom_from_filename(filename) bind(C, name='fftwf_import_wisdom_from_filename') + import + character(C_CHAR), dimension(*), intent(in) :: filename + end function fftwf_import_wisdom_from_filename + + integer(C_INT) function fftwf_import_wisdom_from_file(input_file) bind(C, name='fftwf_import_wisdom_from_file') + import + type(C_PTR), value :: input_file + end function fftwf_import_wisdom_from_file + + integer(C_INT) function fftwf_import_wisdom_from_string(input_string) bind(C, name='fftwf_import_wisdom_from_string') + import + character(C_CHAR), dimension(*), intent(in) :: input_string + end function fftwf_import_wisdom_from_string + + integer(C_INT) function fftwf_import_wisdom(read_char,data) bind(C, name='fftwf_import_wisdom') + import + type(C_FUNPTR), value :: read_char + type(C_PTR), value :: data + end function fftwf_import_wisdom + + subroutine fftwf_fprint_plan(p,output_file) bind(C, name='fftwf_fprint_plan') + import + type(C_PTR), value :: p + type(C_PTR), value :: output_file + end subroutine fftwf_fprint_plan + + subroutine fftwf_print_plan(p) bind(C, name='fftwf_print_plan') + import + type(C_PTR), value :: p + end subroutine fftwf_print_plan + + type(C_PTR) function fftwf_malloc(n) bind(C, name='fftwf_malloc') + import + integer(C_SIZE_T), value :: n + end function fftwf_malloc + + type(C_PTR) function fftwf_alloc_real(n) bind(C, name='fftwf_alloc_real') + import + integer(C_SIZE_T), value :: n + end function fftwf_alloc_real + + type(C_PTR) function fftwf_alloc_complex(n) bind(C, name='fftwf_alloc_complex') + import + integer(C_SIZE_T), value :: n + end function fftwf_alloc_complex + + subroutine fftwf_free(p) bind(C, name='fftwf_free') + import + type(C_PTR), value :: p + end subroutine fftwf_free + + subroutine fftwf_flops(p,add,mul,fmas) bind(C, name='fftwf_flops') + import + type(C_PTR), value :: p + real(C_DOUBLE), intent(out) :: add + real(C_DOUBLE), intent(out) :: mul + real(C_DOUBLE), intent(out) :: fmas + end subroutine fftwf_flops + + real(C_DOUBLE) function fftwf_estimate_cost(p) bind(C, name='fftwf_estimate_cost') + import + type(C_PTR), value :: p + end function fftwf_estimate_cost + + real(C_DOUBLE) function fftwf_cost(p) bind(C, name='fftwf_cost') + import + type(C_PTR), value :: p + end function fftwf_cost + + end interface + + type, bind(C) :: fftwl_iodim + integer(C_INT) n, is, os + end type fftwl_iodim + type, bind(C) :: fftwl_iodim64 + integer(C_INTPTR_T) n, is, os + end type fftwl_iodim64 + + interface + type(C_PTR) function fftwl_plan_dft(rank,n,in,out,sign,flags) bind(C, name='fftwl_plan_dft') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftwl_plan_dft + + type(C_PTR) function fftwl_plan_dft_1d(n,in,out,sign,flags) bind(C, name='fftwl_plan_dft_1d') + import + integer(C_INT), value :: n + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftwl_plan_dft_1d + + type(C_PTR) function fftwl_plan_dft_2d(n0,n1,in,out,sign,flags) bind(C, name='fftwl_plan_dft_2d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftwl_plan_dft_2d + + type(C_PTR) function fftwl_plan_dft_3d(n0,n1,n2,in,out,sign,flags) bind(C, name='fftwl_plan_dft_3d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + integer(C_INT), value :: n2 + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftwl_plan_dft_3d + + type(C_PTR) function fftwl_plan_many_dft(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,sign,flags) & + bind(C, name='fftwl_plan_many_dft') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + integer(C_INT), value :: howmany + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + integer(C_INT), dimension(*), intent(in) :: inembed + integer(C_INT), value :: istride + integer(C_INT), value :: idist + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), dimension(*), intent(in) :: onembed + integer(C_INT), value :: ostride + integer(C_INT), value :: odist + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftwl_plan_many_dft + + type(C_PTR) function fftwl_plan_guru_dft(rank,dims,howmany_rank,howmany_dims,in,out,sign,flags) & + bind(C, name='fftwl_plan_guru_dft') + import + integer(C_INT), value :: rank + type(fftwl_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwl_iodim), dimension(*), intent(in) :: howmany_dims + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftwl_plan_guru_dft + + type(C_PTR) function fftwl_plan_guru_split_dft(rank,dims,howmany_rank,howmany_dims,ri,ii,ro,io,flags) & + bind(C, name='fftwl_plan_guru_split_dft') + import + integer(C_INT), value :: rank + type(fftwl_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwl_iodim), dimension(*), intent(in) :: howmany_dims + real(C_LONG_DOUBLE), dimension(*), intent(out) :: ri + real(C_LONG_DOUBLE), dimension(*), intent(out) :: ii + real(C_LONG_DOUBLE), dimension(*), intent(out) :: ro + real(C_LONG_DOUBLE), dimension(*), intent(out) :: io + integer(C_INT), value :: flags + end function fftwl_plan_guru_split_dft + + type(C_PTR) function fftwl_plan_guru64_dft(rank,dims,howmany_rank,howmany_dims,in,out,sign,flags) & + bind(C, name='fftwl_plan_guru64_dft') + import + integer(C_INT), value :: rank + type(fftwl_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwl_iodim64), dimension(*), intent(in) :: howmany_dims + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: sign + integer(C_INT), value :: flags + end function fftwl_plan_guru64_dft + + type(C_PTR) function fftwl_plan_guru64_split_dft(rank,dims,howmany_rank,howmany_dims,ri,ii,ro,io,flags) & + bind(C, name='fftwl_plan_guru64_split_dft') + import + integer(C_INT), value :: rank + type(fftwl_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwl_iodim64), dimension(*), intent(in) :: howmany_dims + real(C_LONG_DOUBLE), dimension(*), intent(out) :: ri + real(C_LONG_DOUBLE), dimension(*), intent(out) :: ii + real(C_LONG_DOUBLE), dimension(*), intent(out) :: ro + real(C_LONG_DOUBLE), dimension(*), intent(out) :: io + integer(C_INT), value :: flags + end function fftwl_plan_guru64_split_dft + + subroutine fftwl_execute_dft(p,in,out) bind(C, name='fftwl_execute_dft') + import + type(C_PTR), value :: p + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(inout) :: in + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + end subroutine fftwl_execute_dft + + subroutine fftwl_execute_split_dft(p,ri,ii,ro,io) bind(C, name='fftwl_execute_split_dft') + import + type(C_PTR), value :: p + real(C_LONG_DOUBLE), dimension(*), intent(inout) :: ri + real(C_LONG_DOUBLE), dimension(*), intent(inout) :: ii + real(C_LONG_DOUBLE), dimension(*), intent(out) :: ro + real(C_LONG_DOUBLE), dimension(*), intent(out) :: io + end subroutine fftwl_execute_split_dft + + type(C_PTR) function fftwl_plan_many_dft_r2c(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,flags) & + bind(C, name='fftwl_plan_many_dft_r2c') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + integer(C_INT), value :: howmany + real(C_LONG_DOUBLE), dimension(*), intent(out) :: in + integer(C_INT), dimension(*), intent(in) :: inembed + integer(C_INT), value :: istride + integer(C_INT), value :: idist + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), dimension(*), intent(in) :: onembed + integer(C_INT), value :: ostride + integer(C_INT), value :: odist + integer(C_INT), value :: flags + end function fftwl_plan_many_dft_r2c + + type(C_PTR) function fftwl_plan_dft_r2c(rank,n,in,out,flags) bind(C, name='fftwl_plan_dft_r2c') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + real(C_LONG_DOUBLE), dimension(*), intent(out) :: in + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwl_plan_dft_r2c + + type(C_PTR) function fftwl_plan_dft_r2c_1d(n,in,out,flags) bind(C, name='fftwl_plan_dft_r2c_1d') + import + integer(C_INT), value :: n + real(C_LONG_DOUBLE), dimension(*), intent(out) :: in + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwl_plan_dft_r2c_1d + + type(C_PTR) function fftwl_plan_dft_r2c_2d(n0,n1,in,out,flags) bind(C, name='fftwl_plan_dft_r2c_2d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + real(C_LONG_DOUBLE), dimension(*), intent(out) :: in + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwl_plan_dft_r2c_2d + + type(C_PTR) function fftwl_plan_dft_r2c_3d(n0,n1,n2,in,out,flags) bind(C, name='fftwl_plan_dft_r2c_3d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + integer(C_INT), value :: n2 + real(C_LONG_DOUBLE), dimension(*), intent(out) :: in + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwl_plan_dft_r2c_3d + + type(C_PTR) function fftwl_plan_many_dft_c2r(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,flags) & + bind(C, name='fftwl_plan_many_dft_c2r') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + integer(C_INT), value :: howmany + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + integer(C_INT), dimension(*), intent(in) :: inembed + integer(C_INT), value :: istride + integer(C_INT), value :: idist + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), dimension(*), intent(in) :: onembed + integer(C_INT), value :: ostride + integer(C_INT), value :: odist + integer(C_INT), value :: flags + end function fftwl_plan_many_dft_c2r + + type(C_PTR) function fftwl_plan_dft_c2r(rank,n,in,out,flags) bind(C, name='fftwl_plan_dft_c2r') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwl_plan_dft_c2r + + type(C_PTR) function fftwl_plan_dft_c2r_1d(n,in,out,flags) bind(C, name='fftwl_plan_dft_c2r_1d') + import + integer(C_INT), value :: n + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwl_plan_dft_c2r_1d + + type(C_PTR) function fftwl_plan_dft_c2r_2d(n0,n1,in,out,flags) bind(C, name='fftwl_plan_dft_c2r_2d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwl_plan_dft_c2r_2d + + type(C_PTR) function fftwl_plan_dft_c2r_3d(n0,n1,n2,in,out,flags) bind(C, name='fftwl_plan_dft_c2r_3d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + integer(C_INT), value :: n2 + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwl_plan_dft_c2r_3d + + type(C_PTR) function fftwl_plan_guru_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,out,flags) & + bind(C, name='fftwl_plan_guru_dft_r2c') + import + integer(C_INT), value :: rank + type(fftwl_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwl_iodim), dimension(*), intent(in) :: howmany_dims + real(C_LONG_DOUBLE), dimension(*), intent(out) :: in + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwl_plan_guru_dft_r2c + + type(C_PTR) function fftwl_plan_guru_dft_c2r(rank,dims,howmany_rank,howmany_dims,in,out,flags) & + bind(C, name='fftwl_plan_guru_dft_c2r') + import + integer(C_INT), value :: rank + type(fftwl_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwl_iodim), dimension(*), intent(in) :: howmany_dims + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwl_plan_guru_dft_c2r + + type(C_PTR) function fftwl_plan_guru_split_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,ro,io,flags) & + bind(C, name='fftwl_plan_guru_split_dft_r2c') + import + integer(C_INT), value :: rank + type(fftwl_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwl_iodim), dimension(*), intent(in) :: howmany_dims + real(C_LONG_DOUBLE), dimension(*), intent(out) :: in + real(C_LONG_DOUBLE), dimension(*), intent(out) :: ro + real(C_LONG_DOUBLE), dimension(*), intent(out) :: io + integer(C_INT), value :: flags + end function fftwl_plan_guru_split_dft_r2c + + type(C_PTR) function fftwl_plan_guru_split_dft_c2r(rank,dims,howmany_rank,howmany_dims,ri,ii,out,flags) & + bind(C, name='fftwl_plan_guru_split_dft_c2r') + import + integer(C_INT), value :: rank + type(fftwl_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwl_iodim), dimension(*), intent(in) :: howmany_dims + real(C_LONG_DOUBLE), dimension(*), intent(out) :: ri + real(C_LONG_DOUBLE), dimension(*), intent(out) :: ii + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwl_plan_guru_split_dft_c2r + + type(C_PTR) function fftwl_plan_guru64_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,out,flags) & + bind(C, name='fftwl_plan_guru64_dft_r2c') + import + integer(C_INT), value :: rank + type(fftwl_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwl_iodim64), dimension(*), intent(in) :: howmany_dims + real(C_LONG_DOUBLE), dimension(*), intent(out) :: in + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwl_plan_guru64_dft_r2c + + type(C_PTR) function fftwl_plan_guru64_dft_c2r(rank,dims,howmany_rank,howmany_dims,in,out,flags) & + bind(C, name='fftwl_plan_guru64_dft_c2r') + import + integer(C_INT), value :: rank + type(fftwl_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwl_iodim64), dimension(*), intent(in) :: howmany_dims + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: in + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwl_plan_guru64_dft_c2r + + type(C_PTR) function fftwl_plan_guru64_split_dft_r2c(rank,dims,howmany_rank,howmany_dims,in,ro,io,flags) & + bind(C, name='fftwl_plan_guru64_split_dft_r2c') + import + integer(C_INT), value :: rank + type(fftwl_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwl_iodim64), dimension(*), intent(in) :: howmany_dims + real(C_LONG_DOUBLE), dimension(*), intent(out) :: in + real(C_LONG_DOUBLE), dimension(*), intent(out) :: ro + real(C_LONG_DOUBLE), dimension(*), intent(out) :: io + integer(C_INT), value :: flags + end function fftwl_plan_guru64_split_dft_r2c + + type(C_PTR) function fftwl_plan_guru64_split_dft_c2r(rank,dims,howmany_rank,howmany_dims,ri,ii,out,flags) & + bind(C, name='fftwl_plan_guru64_split_dft_c2r') + import + integer(C_INT), value :: rank + type(fftwl_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwl_iodim64), dimension(*), intent(in) :: howmany_dims + real(C_LONG_DOUBLE), dimension(*), intent(out) :: ri + real(C_LONG_DOUBLE), dimension(*), intent(out) :: ii + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), value :: flags + end function fftwl_plan_guru64_split_dft_c2r + + subroutine fftwl_execute_dft_r2c(p,in,out) bind(C, name='fftwl_execute_dft_r2c') + import + type(C_PTR), value :: p + real(C_LONG_DOUBLE), dimension(*), intent(inout) :: in + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(out) :: out + end subroutine fftwl_execute_dft_r2c + + subroutine fftwl_execute_dft_c2r(p,in,out) bind(C, name='fftwl_execute_dft_c2r') + import + type(C_PTR), value :: p + complex(C_LONG_DOUBLE_COMPLEX), dimension(*), intent(inout) :: in + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + end subroutine fftwl_execute_dft_c2r + + subroutine fftwl_execute_split_dft_r2c(p,in,ro,io) bind(C, name='fftwl_execute_split_dft_r2c') + import + type(C_PTR), value :: p + real(C_LONG_DOUBLE), dimension(*), intent(inout) :: in + real(C_LONG_DOUBLE), dimension(*), intent(out) :: ro + real(C_LONG_DOUBLE), dimension(*), intent(out) :: io + end subroutine fftwl_execute_split_dft_r2c + + subroutine fftwl_execute_split_dft_c2r(p,ri,ii,out) bind(C, name='fftwl_execute_split_dft_c2r') + import + type(C_PTR), value :: p + real(C_LONG_DOUBLE), dimension(*), intent(inout) :: ri + real(C_LONG_DOUBLE), dimension(*), intent(inout) :: ii + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + end subroutine fftwl_execute_split_dft_c2r + + type(C_PTR) function fftwl_plan_many_r2r(rank,n,howmany,in,inembed,istride,idist,out,onembed,ostride,odist,kind,flags) & + bind(C, name='fftwl_plan_many_r2r') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + integer(C_INT), value :: howmany + real(C_LONG_DOUBLE), dimension(*), intent(out) :: in + integer(C_INT), dimension(*), intent(in) :: inembed + integer(C_INT), value :: istride + integer(C_INT), value :: idist + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + integer(C_INT), dimension(*), intent(in) :: onembed + integer(C_INT), value :: ostride + integer(C_INT), value :: odist + integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind + integer(C_INT), value :: flags + end function fftwl_plan_many_r2r + + type(C_PTR) function fftwl_plan_r2r(rank,n,in,out,kind,flags) bind(C, name='fftwl_plan_r2r') + import + integer(C_INT), value :: rank + integer(C_INT), dimension(*), intent(in) :: n + real(C_LONG_DOUBLE), dimension(*), intent(out) :: in + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind + integer(C_INT), value :: flags + end function fftwl_plan_r2r + + type(C_PTR) function fftwl_plan_r2r_1d(n,in,out,kind,flags) bind(C, name='fftwl_plan_r2r_1d') + import + integer(C_INT), value :: n + real(C_LONG_DOUBLE), dimension(*), intent(out) :: in + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + integer(C_FFTW_R2R_KIND), value :: kind + integer(C_INT), value :: flags + end function fftwl_plan_r2r_1d + + type(C_PTR) function fftwl_plan_r2r_2d(n0,n1,in,out,kind0,kind1,flags) bind(C, name='fftwl_plan_r2r_2d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + real(C_LONG_DOUBLE), dimension(*), intent(out) :: in + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + integer(C_FFTW_R2R_KIND), value :: kind0 + integer(C_FFTW_R2R_KIND), value :: kind1 + integer(C_INT), value :: flags + end function fftwl_plan_r2r_2d + + type(C_PTR) function fftwl_plan_r2r_3d(n0,n1,n2,in,out,kind0,kind1,kind2,flags) bind(C, name='fftwl_plan_r2r_3d') + import + integer(C_INT), value :: n0 + integer(C_INT), value :: n1 + integer(C_INT), value :: n2 + real(C_LONG_DOUBLE), dimension(*), intent(out) :: in + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + integer(C_FFTW_R2R_KIND), value :: kind0 + integer(C_FFTW_R2R_KIND), value :: kind1 + integer(C_FFTW_R2R_KIND), value :: kind2 + integer(C_INT), value :: flags + end function fftwl_plan_r2r_3d + + type(C_PTR) function fftwl_plan_guru_r2r(rank,dims,howmany_rank,howmany_dims,in,out,kind,flags) & + bind(C, name='fftwl_plan_guru_r2r') + import + integer(C_INT), value :: rank + type(fftwl_iodim), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwl_iodim), dimension(*), intent(in) :: howmany_dims + real(C_LONG_DOUBLE), dimension(*), intent(out) :: in + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind + integer(C_INT), value :: flags + end function fftwl_plan_guru_r2r + + type(C_PTR) function fftwl_plan_guru64_r2r(rank,dims,howmany_rank,howmany_dims,in,out,kind,flags) & + bind(C, name='fftwl_plan_guru64_r2r') + import + integer(C_INT), value :: rank + type(fftwl_iodim64), dimension(*), intent(in) :: dims + integer(C_INT), value :: howmany_rank + type(fftwl_iodim64), dimension(*), intent(in) :: howmany_dims + real(C_LONG_DOUBLE), dimension(*), intent(out) :: in + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + integer(C_FFTW_R2R_KIND), dimension(*), intent(in) :: kind + integer(C_INT), value :: flags + end function fftwl_plan_guru64_r2r + + subroutine fftwl_execute_r2r(p,in,out) bind(C, name='fftwl_execute_r2r') + import + type(C_PTR), value :: p + real(C_LONG_DOUBLE), dimension(*), intent(inout) :: in + real(C_LONG_DOUBLE), dimension(*), intent(out) :: out + end subroutine fftwl_execute_r2r + + subroutine fftwl_destroy_plan(p) bind(C, name='fftwl_destroy_plan') + import + type(C_PTR), value :: p + end subroutine fftwl_destroy_plan + + subroutine fftwl_forget_wisdom() bind(C, name='fftwl_forget_wisdom') + import + end subroutine fftwl_forget_wisdom + + subroutine fftwl_cleanup() bind(C, name='fftwl_cleanup') + import + end subroutine fftwl_cleanup + + subroutine fftwl_set_timelimit(t) bind(C, name='fftwl_set_timelimit') + import + real(C_DOUBLE), value :: t + end subroutine fftwl_set_timelimit + + subroutine fftwl_plan_with_nthreads(nthreads) bind(C, name='fftwl_plan_with_nthreads') + import + integer(C_INT), value :: nthreads + end subroutine fftwl_plan_with_nthreads + + integer(C_INT) function fftwl_init_threads() bind(C, name='fftwl_init_threads') + import + end function fftwl_init_threads + + subroutine fftwl_cleanup_threads() bind(C, name='fftwl_cleanup_threads') + import + end subroutine fftwl_cleanup_threads + + integer(C_INT) function fftwl_export_wisdom_to_filename(filename) bind(C, name='fftwl_export_wisdom_to_filename') + import + character(C_CHAR), dimension(*), intent(in) :: filename + end function fftwl_export_wisdom_to_filename + + subroutine fftwl_export_wisdom_to_file(output_file) bind(C, name='fftwl_export_wisdom_to_file') + import + type(C_PTR), value :: output_file + end subroutine fftwl_export_wisdom_to_file + + type(C_PTR) function fftwl_export_wisdom_to_string() bind(C, name='fftwl_export_wisdom_to_string') + import + end function fftwl_export_wisdom_to_string + + subroutine fftwl_export_wisdom(write_char,data) bind(C, name='fftwl_export_wisdom') + import + type(C_FUNPTR), value :: write_char + type(C_PTR), value :: data + end subroutine fftwl_export_wisdom + + integer(C_INT) function fftwl_import_system_wisdom() bind(C, name='fftwl_import_system_wisdom') + import + end function fftwl_import_system_wisdom + + integer(C_INT) function fftwl_import_wisdom_from_filename(filename) bind(C, name='fftwl_import_wisdom_from_filename') + import + character(C_CHAR), dimension(*), intent(in) :: filename + end function fftwl_import_wisdom_from_filename + + integer(C_INT) function fftwl_import_wisdom_from_file(input_file) bind(C, name='fftwl_import_wisdom_from_file') + import + type(C_PTR), value :: input_file + end function fftwl_import_wisdom_from_file + + integer(C_INT) function fftwl_import_wisdom_from_string(input_string) bind(C, name='fftwl_import_wisdom_from_string') + import + character(C_CHAR), dimension(*), intent(in) :: input_string + end function fftwl_import_wisdom_from_string + + integer(C_INT) function fftwl_import_wisdom(read_char,data) bind(C, name='fftwl_import_wisdom') + import + type(C_FUNPTR), value :: read_char + type(C_PTR), value :: data + end function fftwl_import_wisdom + + subroutine fftwl_fprint_plan(p,output_file) bind(C, name='fftwl_fprint_plan') + import + type(C_PTR), value :: p + type(C_PTR), value :: output_file + end subroutine fftwl_fprint_plan + + subroutine fftwl_print_plan(p) bind(C, name='fftwl_print_plan') + import + type(C_PTR), value :: p + end subroutine fftwl_print_plan + + type(C_PTR) function fftwl_malloc(n) bind(C, name='fftwl_malloc') + import + integer(C_SIZE_T), value :: n + end function fftwl_malloc + + type(C_PTR) function fftwl_alloc_real(n) bind(C, name='fftwl_alloc_real') + import + integer(C_SIZE_T), value :: n + end function fftwl_alloc_real + + type(C_PTR) function fftwl_alloc_complex(n) bind(C, name='fftwl_alloc_complex') + import + integer(C_SIZE_T), value :: n + end function fftwl_alloc_complex + + subroutine fftwl_free(p) bind(C, name='fftwl_free') + import + type(C_PTR), value :: p + end subroutine fftwl_free + + subroutine fftwl_flops(p,add,mul,fmas) bind(C, name='fftwl_flops') + import + type(C_PTR), value :: p + real(C_DOUBLE), intent(out) :: add + real(C_DOUBLE), intent(out) :: mul + real(C_DOUBLE), intent(out) :: fmas + end subroutine fftwl_flops + + real(C_DOUBLE) function fftwl_estimate_cost(p) bind(C, name='fftwl_estimate_cost') + import + type(C_PTR), value :: p + end function fftwl_estimate_cost + + real(C_DOUBLE) function fftwl_cost(p) bind(C, name='fftwl_cost') + import + type(C_PTR), value :: p + end function fftwl_cost + + end interface diff --git a/code/kdtree2.f90 b/code/kdtree2.f90 deleted file mode 100644 index 422cef447..000000000 --- a/code/kdtree2.f90 +++ /dev/null @@ -1,1885 +0,0 @@ -! -!(c) Matthew Kennel, Institute for Nonlinear Science (2004) -! -! Licensed under the Academic Free License version 1.1 found in file LICENSE -! with additional provisions found in that same file. -! -!####################################################### -! modifications: changed precision according to prec.f90 -! k.komerla, m.diehl -!####################################################### - -module kdtree2_priority_queue_module - use prec - ! - ! maintain a priority queue (PQ) of data, pairs of 'priority/payload', - ! implemented with a binary heap. This is the type, and the 'dis' field - ! is the priority. - ! - type kdtree2_result - ! a pair of distances, indexes - real(pReal) :: dis !=0.0 - integer(pInt) :: idx !=-1 Initializers cause some bugs in compilers. - end type kdtree2_result - ! - ! A heap-based priority queue lets one efficiently implement the following - ! operations, each in log(N) time, as opposed to linear time. - ! - ! 1) add a datum (push a datum onto the queue, increasing its length) - ! 2) return the priority value of the maximum priority element - ! 3) pop-off (and delete) the element with the maximum priority, decreasing - ! the size of the queue. - ! 4) replace the datum with the maximum priority with a supplied datum - ! (of either higher or lower priority), maintaining the size of the - ! queue. - ! - ! - ! In the k-d tree case, the 'priority' is the square distance of a point in - ! the data set to a reference point. The goal is to keep the smallest M - ! distances to a reference point. The tree algorithm searches terminal - ! nodes to decide whether to add points under consideration. - ! - ! A priority queue is useful here because it lets one quickly return the - ! largest distance currently existing in the list. If a new candidate - ! distance is smaller than this, then the new candidate ought to replace - ! the old candidate. In priority queue terms, this means removing the - ! highest priority element, and inserting the new one. - ! - ! Algorithms based on Cormen, Leiserson, Rivest, _Introduction - ! to Algorithms_, 1990, with further optimization by the author. - ! - ! Originally informed by a C implementation by Sriranga Veeraraghavan. - ! - ! This module is not written in the most clear way, but is implemented such - ! for speed, as it its operations will be called many times during searches - ! of large numbers of neighbors. - ! - type pq - ! - ! The priority queue consists of elements - ! priority(1:heap_size), with associated payload(:). - ! - ! There are heap_size active elements. - ! Assumes the allocation is always sufficient. Will NOT increase it - ! to match. - integer(pInt) :: heap_size = 0 - type(kdtree2_result), pointer :: elems(:) - end type pq - - public :: kdtree2_result - - public :: pq - public :: pq_create - public :: pq_delete, pq_insert - public :: pq_extract_max, pq_max, pq_replace_max, pq_maxpri - private - -contains - - - function pq_create(results_in) result(res) - ! - ! Create a priority queue from ALREADY allocated - ! array pointers for storage. NOTE! It will NOT - ! add any alements to the heap, i.e. any existing - ! data in the input arrays will NOT be used and may - ! be overwritten. - ! - ! usage: - ! real(pReal), pointer :: x(:) - ! integer(pInt), pointer :: k(:) - ! allocate(x(1000),k(1000)) - ! pq => pq_create(x,k) - ! - type(kdtree2_result), target:: results_in(:) - type(pq) :: res - ! - ! - integer(pInt) :: nalloc - - nalloc = size(results_in,1) - if (nalloc .lt. 1) then - write (*,*) 'PQ_CREATE: error, input arrays must be allocated.' - end if - res%elems => results_in - res%heap_size = 0 - return - end function pq_create - - ! - ! operations for getting parents and left + right children - ! of elements in a binary heap. - ! - -! -! These are written inline for speed. -! -! integer(pInt) function parent(i) -! integer(pInt), intent(in) :: i -! parent = (i/2) -! return -! end function parent - -! integer(pInt) function left(i) -! integer(pInt), intent(in) ::i -! left = (2*i) -! return -! end function left - -! integer(pInt) function right(i) -! integer(pInt), intent(in) :: i -! right = (2*i)+1 -! return -! end function right - -! logical function compare_priority(p1,p2) -! real(pReal), intent(in) :: p1, p2 -! -! compare_priority = (p1 .gt. p2) -! return -! end function compare_priority - - subroutine heapify(a,i_in) - ! - ! take a heap rooted at 'i' and force it to be in the - ! heap canonical form. This is performance critical - ! and has been tweaked a little to reflect this. - ! - type(pq),pointer :: a - integer(pInt), intent(in) :: i_in - ! - integer(pInt) :: i, l, r, largest - - real(pReal) :: pri_i, pri_l, pri_r, pri_largest - - - type(kdtree2_result) :: temp - - i = i_in - -bigloop: do - l = 2*i ! left(i) - r = l+1 ! right(i) - ! - ! set 'largest' to the index of either i, l, r - ! depending on whose priority is largest. - ! - ! note that l or r can be larger than the heap size - ! in which case they do not count. - - - ! does left child have higher priority? - if (l .gt. a%heap_size) then - ! we know that i is the largest as both l and r are invalid. - exit - else - pri_i = a%elems(i)%dis - pri_l = a%elems(l)%dis - if (pri_l .gt. pri_i) then - largest = l - pri_largest = pri_l - else - largest = i - pri_largest = pri_i - endif - - ! - ! between i and l we have a winner - ! now choose between that and r. - ! - if (r .le. a%heap_size) then - pri_r = a%elems(r)%dis - if (pri_r .gt. pri_largest) then - largest = r - endif - endif - endif - - if (largest .ne. i) then - ! swap data in nodes largest and i, then heapify - - temp = a%elems(i) - a%elems(i) = a%elems(largest) - a%elems(largest) = temp - ! - ! Canonical heapify() algorithm has tail-ecursive call: - ! - ! call heapify(a,largest) - ! we will simulate with cycle - ! - i = largest - cycle bigloop ! continue the loop - else - return ! break from the loop - end if - enddo bigloop - return - end subroutine heapify - - subroutine pq_max(a,e) - ! - ! return the priority and its payload of the maximum priority element - ! on the queue, which should be the first one, if it is - ! in heapified form. - ! - type(pq),pointer :: a - type(kdtree2_result),intent(out) :: e - - if (a%heap_size .gt. 0) then - e = a%elems(1) - else - write (*,*) 'PQ_MAX: ERROR, heap_size < 1' - stop - endif - return - end subroutine pq_max - - real(pReal) function pq_maxpri(a) - type(pq), pointer :: a - - if (a%heap_size .gt. 0) then - pq_maxpri = a%elems(1)%dis - else - write (*,*) 'PQ_MAX_PRI: ERROR, heapsize < 1' - stop - endif - return - end function pq_maxpri - - subroutine pq_extract_max(a,e) - ! - ! return the priority and payload of maximum priority - ! element, and remove it from the queue. - ! (equivalent to 'pop()' on a stack) - ! - type(pq),pointer :: a - type(kdtree2_result), intent(out) :: e - - if (a%heap_size .ge. 1) then - ! - ! return max as first element - ! - e = a%elems(1) - - ! - ! move last element to first - ! - a%elems(1) = a%elems(a%heap_size) - a%heap_size = a%heap_size-1 - call heapify(a,1) - return - else - write (*,*) 'PQ_EXTRACT_MAX: error, attempted to pop non-positive PQ' - stop - end if - - end subroutine pq_extract_max - - - real(pReal) function pq_insert(a,dis,idx) - ! - ! Insert a new element and return the new maximum priority, - ! which may or may not be the same as the old maximum priority. - ! - type(pq),pointer :: a - real(pReal), intent(in) :: dis - integer(pInt), intent(in) :: idx - ! type(kdtree2_result), intent(in) :: e - ! - integer(pInt) :: i, isparent - real(pReal) :: parentdis - ! - - ! if (a%heap_size .ge. a%max_elems) then - ! write (*,*) 'PQ_INSERT: error, attempt made to insert element on full PQ' - ! stop - ! else - a%heap_size = a%heap_size + 1 - i = a%heap_size - - do while (i .gt. 1) - isparent = int(i/2) - parentdis = a%elems(isparent)%dis - if (dis .gt. parentdis) then - ! move what was in i's parent into i. - a%elems(i)%dis = parentdis - a%elems(i)%idx = a%elems(isparent)%idx - i = isparent - else - exit - endif - end do - - ! insert the element at the determined position - a%elems(i)%dis = dis - a%elems(i)%idx = idx - - pq_insert = a%elems(1)%dis - return - ! end if - - end function pq_insert - - subroutine pq_adjust_heap(a,i) - type(pq),pointer :: a - integer(pInt), intent(in) :: i - ! - ! nominally arguments (a,i), but specialize for a=1 - ! - ! This routine assumes that the trees with roots 2 and 3 are already heaps, i.e. - ! the children of '1' are heaps. When the procedure is completed, the - ! tree rooted at 1 is a heap. - real(pReal) :: prichild - integer(pInt) :: parent, child, N - - type(kdtree2_result) :: e - - e = a%elems(i) - - parent = i - child = 2*i - N = a%heap_size - - do while (child .le. N) - if (child .lt. N) then - if (a%elems(child)%dis .lt. a%elems(child+1)%dis) then - child = child+1 - endif - endif - prichild = a%elems(child)%dis - if (e%dis .ge. prichild) then - exit - else - ! move child into parent. - a%elems(parent) = a%elems(child) - parent = child - child = 2*parent - end if - end do - a%elems(parent) = e - return - end subroutine pq_adjust_heap - - - real(pReal) function pq_replace_max(a,dis,idx) - ! - ! Replace the extant maximum priority element - ! in the PQ with (dis,idx). Return - ! the new maximum priority, which may be larger - ! or smaller than the old one. - ! - type(pq),pointer :: a - real(pReal), intent(in) :: dis - integer(pInt), intent(in) :: idx -! type(kdtree2_result), intent(in) :: e - ! not tested as well! - - integer(pInt) :: parent, child, N - real(pReal) :: prichild, prichildp1 - - type(kdtree2_result) :: etmp - - if (.true.) then - N=a%heap_size - if (N .ge. 1) then - parent =1 - child=2 - - loop: do while (child .le. N) - prichild = a%elems(child)%dis - - ! - ! posibly child+1 has higher priority, and if - ! so, get it, and increment child. - ! - - if (child .lt. N) then - prichildp1 = a%elems(child+1)%dis - if (prichild .lt. prichildp1) then - child = child+1 - prichild = prichildp1 - endif - endif - - if (dis .ge. prichild) then - exit loop - ! we have a proper place for our new element, - ! bigger than either children's priority. - else - ! move child into parent. - a%elems(parent) = a%elems(child) - parent = child - child = 2*parent - end if - end do loop - a%elems(parent)%dis = dis - a%elems(parent)%idx = idx - pq_replace_max = a%elems(1)%dis - else - a%elems(1)%dis = dis - a%elems(1)%idx = idx - pq_replace_max = dis - endif - else - ! - ! slower version using elementary pop and push operations. - ! - call pq_extract_max(a,etmp) - etmp%dis = dis - etmp%idx = idx - pq_replace_max = pq_insert(a,dis,idx) - endif - return - end function pq_replace_max - - subroutine pq_delete(a,i) - ! - ! delete item with index 'i' - ! - type(pq),pointer :: a - integer(pInt) :: i - - if ((i .lt. 1) .or. (i .gt. a%heap_size)) then - write (*,*) 'PQ_DELETE: error, attempt to remove out of bounds element.' - stop - endif - - ! swap the item to be deleted with the last element - ! and shorten heap by one. - a%elems(i) = a%elems(a%heap_size) - a%heap_size = a%heap_size - 1 - - call heapify(a,i) - - end subroutine pq_delete - -end module kdtree2_priority_queue_module - - -module kdtree2_module - use prec - use kdtree2_priority_queue_module - ! K-D tree routines in Fortran 90 by Matt Kennel. - ! Original program was written in Sather by Steve Omohundro and - ! Matt Kennel. Only the Euclidean metric is supported. - ! - ! - ! This module is identical to 'kd_tree', except that the order - ! of subscripts is reversed in the data file. - ! In otherwords for an embedding of N D-dimensional vectors, the - ! data file is here, in natural Fortran order data(1:D, 1:N) - ! because Fortran lays out columns first, - ! - ! whereas conventionally (C-style) it is data(1:N,1:D) - ! as in the original kd_tree module. - ! - !-------------DATA TYPE, CREATION, DELETION--------------------- - public :: pReal - public :: kdtree2, kdtree2_result, tree_node, kdtree2_create, kdtree2_destroy - !--------------------------------------------------------------- - !-------------------SEARCH ROUTINES----------------------------- - public :: kdtree2_n_nearest,kdtree2_n_nearest_around_point - ! Return fixed number of nearest neighbors around arbitrary vector, - ! or extant point in dataset, with decorrelation window. - ! - public :: kdtree2_r_nearest, kdtree2_r_nearest_around_point - ! Return points within a fixed ball of arb vector/extant point - ! - public :: kdtree2_sort_results - ! Sort, in order of increasing distance, rseults from above. - ! - public :: kdtree2_r_count, kdtree2_r_count_around_point - ! Count points within a fixed ball of arb vector/extant point - ! - public :: kdtree2_n_nearest_brute_force, kdtree2_r_nearest_brute_force - ! brute force of kdtree2_[n|r]_nearest - !---------------------------------------------------------------- - - - integer(pInt), parameter :: bucket_size = 12 - ! The maximum number of points to keep in a terminal node. - - type interval - real(pReal) :: lower,upper - end type interval - - type :: tree_node - ! an internal tree node - private - integer(pInt) :: cut_dim - ! the dimension to cut - real(pReal) :: cut_val - ! where to cut the dimension - real(pReal) :: cut_val_left, cut_val_right - ! improved cutoffs knowing the spread in child boxes. - integer(pInt) :: l, u - type (tree_node), pointer :: left, right - type(interval), pointer :: box(:) => null() - ! child pointers - ! Points included in this node are indexes[k] with k \in [l,u] - - - end type tree_node - - type :: kdtree2 - ! Global information about the tree, one per tree - integer(pInt) :: dimen=0, n=0 - ! dimensionality and total # of points - real(pReal), pointer :: the_data(:,:) => null() - ! pointer to the actual data array - ! - ! IMPORTANT NOTE: IT IS DIMENSIONED the_data(1:d,1:N) - ! which may be opposite of what may be conventional. - ! This is, because in Fortran, the memory layout is such that - ! the first dimension is in sequential order. Hence, with - ! (1:d,1:N), all components of the vector will be in consecutive - ! memory locations. The search time is dominated by the - ! evaluation of distances in the terminal nodes. Putting all - ! vector components in consecutive memory location improves - ! memory cache locality, and hence search speed, and may enable - ! vectorization on some processors and compilers. - - integer(pInt), pointer :: ind(:) => null() - ! permuted index into the data, so that indexes[l..u] of some - ! bucket represent the indexes of the actual points in that - ! bucket. - logical :: sort = .false. - ! do we always sort output results? - logical :: rearrange = .false. - real(pReal), pointer :: rearranged_data(:,:) => null() - ! if (rearrange .eqv. .true.) then rearranged_data has been - ! created so that rearranged_data(:,i) = the_data(:,ind(i)), - ! permitting search to use more cache-friendly rearranged_data, at - ! some initial computation and storage cost. - type (tree_node), pointer :: root => null() - ! root pointer of the tree - end type kdtree2 - - - type :: tree_search_record - ! - ! One of these is created for each search. - ! - private - ! - ! Many fields are copied from the tree structure, in order to - ! speed up the search. - ! - integer(pInt) :: dimen - integer(pInt) :: nn, nfound - real(pReal) :: ballsize - integer(pInt) :: centeridx=999, correltime=9999 - ! exclude points within 'correltime' of 'centeridx', iff centeridx >= 0 - integer(pInt) :: nalloc ! how much allocated for results(:)? - logical :: rearrange ! are the data rearranged or original? - ! did the # of points found overflow the storage provided? - logical :: overflow - real(pReal), pointer :: qv(:) ! query vector - type(kdtree2_result), pointer :: results(:) ! results - type(pq) :: pq - real(pReal), pointer :: data(:,:) ! temp pointer to data - integer(pInt), pointer :: ind(:) ! temp pointer to indexes - end type tree_search_record - - private - ! everything else is private. - - type(tree_search_record), save, target :: sr ! A GLOBAL VARIABLE for search - -contains - - function kdtree2_create(input_data,dim,sort,rearrange) result (mr) - ! - ! create the actual tree structure, given an input array of data. - ! - ! Note, input data is input_data(1:d,1:N), NOT the other way around. - ! THIS IS THE REVERSE OF THE PREVIOUS VERSION OF THIS MODULE. - ! The reason for it is cache friendliness, improving performance. - ! - ! Optional arguments: If 'dim' is specified, then the tree - ! will only search the first 'dim' components - ! of input_data, otherwise, dim is inferred - ! from SIZE(input_data,1). - ! - ! if sort .eqv. .true. then output results - ! will be sorted by increasing distance. - ! default=.false., as it is faster to not sort. - ! - ! if rearrange .eqv. .true. then an internal - ! copy of the data, rearranged by terminal node, - ! will be made for cache friendliness. - ! default=.true., as it speeds searches, but - ! building takes longer, and extra memory is used. - ! - ! .. Function Return Cut_value .. - type (kdtree2), pointer :: mr - integer(pInt), intent(in), optional :: dim - logical, intent(in), optional :: sort - logical, intent(in), optional :: rearrange - ! .. - ! .. Array Arguments .. - real(pReal), target :: input_data(:,:) - ! - integer(pInt) :: i - ! .. - allocate (mr) - mr%the_data => input_data - ! pointer assignment - - if (present(dim)) then - mr%dimen = dim - else - mr%dimen = size(input_data,1) - end if - mr%n = size(input_data,2) - - if (mr%dimen > mr%n) then - ! unlikely to be correct - write (*,*) 'KD_TREE_TRANS: likely user error.' - write (*,*) 'KD_TREE_TRANS: You passed in matrix with D=',mr%dimen - write (*,*) 'KD_TREE_TRANS: and N=',mr%n - write (*,*) 'KD_TREE_TRANS: note, that new format is data(1:D,1:N)' - write (*,*) 'KD_TREE_TRANS: with usually N >> D. If N =approx= D, then a k-d tree' - write (*,*) 'KD_TREE_TRANS: is not an appropriate data structure.' - stop - end if - - call build_tree(mr) - - if (present(sort)) then - mr%sort = sort - else - mr%sort = .false. - endif - - if (present(rearrange)) then - mr%rearrange = rearrange - else - mr%rearrange = .true. - endif - - if (mr%rearrange) then - allocate(mr%rearranged_data(mr%dimen,mr%n)) - do i=1,mr%n - mr%rearranged_data(:,i) = mr%the_data(:, & - mr%ind(i)) - enddo - else - nullify(mr%rearranged_data) - endif - - end function kdtree2_create - - subroutine build_tree(tp) - type (kdtree2), pointer :: tp - ! .. - integer(pInt) :: j - type(tree_node), pointer :: dummy => null() - ! .. - allocate (tp%ind(tp%n)) - forall (j=1:tp%n) - tp%ind(j) = j - end forall - tp%root => build_tree_for_range(tp,1,tp%n, dummy) - end subroutine build_tree - - recursive function build_tree_for_range(tp,l,u,parent) result (res) - ! .. Function Return Cut_value .. - type (tree_node), pointer :: res - ! .. - ! .. Structure Arguments .. - type (kdtree2), pointer :: tp - type (tree_node),pointer :: parent - ! .. - ! .. Scalar Arguments .. - integer(pInt), intent (In) :: l, u - ! .. - ! .. Local Scalars .. - integer(pInt) :: i, c, m, dimen - logical :: recompute - real(pReal) :: average - -!!$ If (.False.) Then -!!$ If ((l .Lt. 1) .Or. (l .Gt. tp%n)) Then -!!$ Stop 'illegal L value in build_tree_for_range' -!!$ End If -!!$ If ((u .Lt. 1) .Or. (u .Gt. tp%n)) Then -!!$ Stop 'illegal u value in build_tree_for_range' -!!$ End If -!!$ If (u .Lt. l) Then -!!$ Stop 'U is less than L, thats illegal.' -!!$ End If -!!$ Endif -!!$ - ! first compute min and max - dimen = tp%dimen - allocate (res) - allocate(res%box(dimen)) - - ! First, compute an APPROXIMATE bounding box of all points associated with this node. - if ( u < l ) then - ! no points in this box - nullify(res) - return - end if - - if ((u-l)<=bucket_size) then - ! - ! always compute true bounding box for terminal nodes. - ! - do i=1,dimen - call spread_in_coordinate(tp,i,l,u,res%box(i)) - end do - res%cut_dim = 0 - res%cut_val = 0.0 - res%l = l - res%u = u - res%left =>null() - res%right => null() - else - ! - ! modify approximate bounding box. This will be an - ! overestimate of the true bounding box, as we are only recomputing - ! the bounding box for the dimension that the parent split on. - ! - ! Going to a true bounding box computation would significantly - ! increase the time necessary to build the tree, and usually - ! has only a very small difference. This box is not used - ! for searching but only for deciding which coordinate to split on. - ! - do i=1,dimen - recompute=.true. - if (associated(parent)) then - if (i .ne. parent%cut_dim) then - recompute=.false. - end if - endif - if (recompute) then - call spread_in_coordinate(tp,i,l,u,res%box(i)) - else - res%box(i) = parent%box(i) - endif - end do - - - c = maxloc(res%box(1:dimen)%upper-res%box(1:dimen)%lower,1) - ! - ! c is the identity of which coordinate has the greatest spread. - ! - - if (.false.) then - ! select exact median to have fully balanced tree. - m = (l+u)/2 - call select_on_coordinate(tp%the_data,tp%ind,c,m,l,u) - else - ! - ! select point halfway between min and max, as per A. Moore, - ! who says this helps in some degenerate cases, or - ! actual arithmetic average. - ! - if (.true.) then - ! actually compute average - average = sum(tp%the_data(c,tp%ind(l:u))) / real(u-l+1,pReal) - else - average = (res%box(c)%upper + res%box(c)%lower)/2.0 - endif - - res%cut_val = average - m = select_on_coordinate_value(tp%the_data,tp%ind,c,average,l,u) - endif - - ! moves indexes around - res%cut_dim = c - res%l = l - res%u = u -! res%cut_val = tp%the_data(c,tp%ind(m)) - - res%left => build_tree_for_range(tp,l,m,res) - res%right => build_tree_for_range(tp,m+1,u,res) - - if (associated(res%right) .eqv. .false.) then - res%box = res%left%box - res%cut_val_left = res%left%box(c)%upper - res%cut_val = res%cut_val_left - elseif (associated(res%left) .eqv. .false.) then - res%box = res%right%box - res%cut_val_right = res%right%box(c)%lower - res%cut_val = res%cut_val_right - else - res%cut_val_right = res%right%box(c)%lower - res%cut_val_left = res%left%box(c)%upper - res%cut_val = (res%cut_val_left + res%cut_val_right)/2 - - - ! now remake the true bounding box for self. - ! Since we are taking unions (in effect) of a tree structure, - ! this is much faster than doing an exhaustive - ! search over all points - res%box%upper = max(res%left%box%upper,res%right%box%upper) - res%box%lower = min(res%left%box%lower,res%right%box%lower) - endif - end if - end function build_tree_for_range - - integer(pInt) function select_on_coordinate_value(v,ind,c,alpha,li,ui) & - result(res) - ! Move elts of ind around between l and u, so that all points - ! <= than alpha (in c cooordinate) are first, and then - ! all points > alpha are second. - - ! - ! Algorithm (matt kennel). - ! - ! Consider the list as having three parts: on the left, - ! the points known to be <= alpha. On the right, the points - ! known to be > alpha, and in the middle, the currently unknown - ! points. The algorithm is to scan the unknown points, starting - ! from the left, and swapping them so that they are added to - ! the left stack or the right stack, as appropriate. - ! - ! The algorithm finishes when the unknown stack is empty. - ! - ! .. Scalar Arguments .. - integer(pInt), intent (In) :: c, li, ui - real(pReal), intent(in) :: alpha - ! .. - real(pReal) :: v(1:,1:) - integer(pInt) :: ind(1:) - integer(pInt) :: tmp - ! .. - integer(pInt) :: lb, rb - ! - ! The points known to be <= alpha are in - ! [l,lb-1] - ! - ! The points known to be > alpha are in - ! [rb+1,u]. - ! - ! Therefore we add new points into lb or - ! rb as appropriate. When lb=rb - ! we are done. We return the location of the last point <= alpha. - ! - ! - lb = li; rb = ui - - do while (lb < rb) - if ( v(c,ind(lb)) <= alpha ) then - ! it is good where it is. - lb = lb+1 - else - ! swap it with rb. - tmp = ind(lb); ind(lb) = ind(rb); ind(rb) = tmp - rb = rb-1 - endif - end do - - ! now lb .eq. ub - if (v(c,ind(lb)) <= alpha) then - res = lb - else - res = lb-1 - endif - - end function select_on_coordinate_value - - subroutine select_on_coordinate(v,ind,c,k,li,ui) - ! Move elts of ind around between l and u, so that the kth - ! element - ! is >= those below, <= those above, in the coordinate c. - ! .. Scalar Arguments .. - integer(pInt), intent (In) :: c, k, li, ui - ! .. - integer(pInt) :: i, l, m, s, t, u - ! .. - real(pReal) :: v(:,:) - integer(pInt) :: ind(:) - ! .. - l = li - u = ui - do while (l=k) u = m - 1 - end do - end subroutine select_on_coordinate - - subroutine spread_in_coordinate(tp,c,l,u,interv) - ! the spread in coordinate 'c', between l and u. - ! - ! Return lower bound in 'smin', and upper in 'smax', - ! .. - ! .. Structure Arguments .. - type (kdtree2), pointer :: tp - type(interval), intent(out) :: interv - ! .. - ! .. Scalar Arguments .. - integer(pInt), intent (In) :: c, l, u - ! .. - ! .. Local Scalars .. - real(pReal) :: last, lmax, lmin, t, smin,smax - integer(pInt) :: i, ulocal - ! .. - ! .. Local Arrays .. - real(pReal), pointer :: v(:,:) - integer(pInt), pointer :: ind(:) - ! .. - v => tp%the_data(1:,1:) - ind => tp%ind(1:) - smin = v(c,ind(l)) - smax = smin - - ulocal = u - - do i = l + 2, ulocal, 2 - lmin = v(c,ind(i-1)) - lmax = v(c,ind(i)) - if (lmin>lmax) then - t = lmin - lmin = lmax - lmax = t - end if - if (smin>lmin) smin = lmin - if (smaxlast) smin = last - if (smax qv - sr%nn = nn - sr%nfound = 0 - sr%centeridx = -1 - sr%correltime = 0 - sr%overflow = .false. - - sr%results => results - - sr%nalloc = nn ! will be checked - - sr%ind => tp%ind - sr%rearrange = tp%rearrange - if (tp%rearrange) then - sr%Data => tp%rearranged_data - else - sr%Data => tp%the_data - endif - sr%dimen = tp%dimen - - call validate_query_storage(nn) - sr%pq = pq_create(results) - - call search(tp%root) - - if (tp%sort) then - call kdtree2_sort_results(nn, results) - endif -! deallocate(sr%pqp) - return - end subroutine kdtree2_n_nearest - - subroutine kdtree2_n_nearest_around_point(tp,idxin,correltime,nn,results) - ! Find the 'nn' vectors in the tree nearest to point 'idxin', - ! with correlation window 'correltime', returing results in - ! results(:), which must be pre-allocated upon entry. - type (kdtree2), pointer :: tp - integer(pInt), intent (In) :: idxin, correltime, nn - type(kdtree2_result), target :: results(:) - - allocate (sr%qv(tp%dimen)) - sr%qv = tp%the_data(:,idxin) ! copy the vector - sr%ballsize = huge(1.0) ! the largest real(pReal) number - sr%centeridx = idxin - sr%correltime = correltime - - sr%nn = nn - sr%nfound = 0 - - sr%dimen = tp%dimen - sr%nalloc = nn - - sr%results => results - - sr%ind => tp%ind - sr%rearrange = tp%rearrange - - if (sr%rearrange) then - sr%Data => tp%rearranged_data - else - sr%Data => tp%the_data - endif - - call validate_query_storage(nn) - sr%pq = pq_create(results) - - call search(tp%root) - - if (tp%sort) then - call kdtree2_sort_results(nn, results) - endif - deallocate (sr%qv) - return - end subroutine kdtree2_n_nearest_around_point - - subroutine kdtree2_r_nearest(tp,qv,r2,nfound,nalloc,results) - ! find the nearest neighbors to point 'idxin', within SQUARED - ! Euclidean distance 'r2'. Upon ENTRY, nalloc must be the - ! size of memory allocated for results(1:nalloc). Upon - ! EXIT, nfound is the number actually found within the ball. - ! - ! Note that if nfound .gt. nalloc then more neighbors were found - ! than there were storage to store. The resulting list is NOT - ! the smallest ball inside norm r^2 - ! - ! Results are NOT sorted unless tree was created with sort option. - type (kdtree2), pointer :: tp - real(pReal), target, intent (In) :: qv(:) - real(pReal), intent(in) :: r2 - integer(pInt), intent(out) :: nfound - integer(pInt), intent (In) :: nalloc - type(kdtree2_result), target :: results(:) - - ! - sr%qv => qv - sr%ballsize = r2 - sr%nn = 0 ! flag for fixed ball search - sr%nfound = 0 - sr%centeridx = -1 - sr%correltime = 0 - - sr%results => results - - call validate_query_storage(nalloc) - sr%nalloc = nalloc - sr%overflow = .false. - sr%ind => tp%ind - sr%rearrange= tp%rearrange - - if (tp%rearrange) then - sr%Data => tp%rearranged_data - else - sr%Data => tp%the_data - endif - sr%dimen = tp%dimen - - ! - !sr%dsl = Huge(sr%dsl) ! set to huge positive values - !sr%il = -1 ! set to invalid indexes - ! - - call search(tp%root) - nfound = sr%nfound - if (tp%sort) then - call kdtree2_sort_results(nfound, results) - endif - - if (sr%overflow) then - write (*,*) 'KD_TREE_TRANS: warning! return from kdtree2_r_nearest found more neighbors' - write (*,*) 'KD_TREE_TRANS: than storage was provided for. Answer is NOT smallest ball' - write (*,*) 'KD_TREE_TRANS: with that number of neighbors! I.e. it is wrong.' - endif - - return - end subroutine kdtree2_r_nearest - - subroutine kdtree2_r_nearest_around_point(tp,idxin,correltime,r2,& - nfound,nalloc,results) - ! - ! Like kdtree2_r_nearest, but around a point 'idxin' already existing - ! in the data set. - ! - ! Results are NOT sorted unless tree was created with sort option. - ! - type (kdtree2), pointer :: tp - integer(pInt), intent (In) :: idxin, correltime, nalloc - real(pReal), intent(in) :: r2 - integer(pInt), intent(out) :: nfound - type(kdtree2_result), target :: results(:) - ! .. - ! .. Intrinsic Functions .. - intrinsic HUGE - ! .. - allocate (sr%qv(tp%dimen)) - sr%qv = tp%the_data(:,idxin) ! copy the vector - sr%ballsize = r2 - sr%nn = 0 ! flag for fixed r search - sr%nfound = 0 - sr%centeridx = idxin - sr%correltime = correltime - - sr%results => results - - sr%nalloc = nalloc - sr%overflow = .false. - - call validate_query_storage(nalloc) - - ! sr%dsl = HUGE(sr%dsl) ! set to huge positive values - ! sr%il = -1 ! set to invalid indexes - - sr%ind => tp%ind - sr%rearrange = tp%rearrange - - if (tp%rearrange) then - sr%Data => tp%rearranged_data - else - sr%Data => tp%the_data - endif - sr%rearrange = tp%rearrange - sr%dimen = tp%dimen - - ! - !sr%dsl = Huge(sr%dsl) ! set to huge positive values - !sr%il = -1 ! set to invalid indexes - ! - - call search(tp%root) - nfound = sr%nfound - if (tp%sort) then - call kdtree2_sort_results(nfound,results) - endif - - if (sr%overflow) then - write (*,*) 'KD_TREE_TRANS: warning! return from kdtree2_r_nearest found more neighbors' - write (*,*) 'KD_TREE_TRANS: than storage was provided for. Answer is NOT smallest ball' - write (*,*) 'KD_TREE_TRANS: with that number of neighbors! I.e. it is wrong.' - endif - - deallocate (sr%qv) - return - end subroutine kdtree2_r_nearest_around_point - - function kdtree2_r_count(tp,qv,r2) result(nfound) - ! Count the number of neighbors within square distance 'r2'. - type (kdtree2), pointer :: tp - real(pReal), target, intent (In) :: qv(:) - real(pReal), intent(in) :: r2 - integer(pInt) :: nfound - ! .. - ! .. Intrinsic Functions .. - intrinsic HUGE - ! .. - sr%qv => qv - sr%ballsize = r2 - - sr%nn = 0 ! flag for fixed r search - sr%nfound = 0 - sr%centeridx = -1 - sr%correltime = 0 - - nullify(sr%results) ! for some reason, FTN 95 chokes on '=> null()' - - sr%nalloc = 0 ! we do not allocate any storage but that's OK - ! for counting. - sr%ind => tp%ind - sr%rearrange = tp%rearrange - if (tp%rearrange) then - sr%Data => tp%rearranged_data - else - sr%Data => tp%the_data - endif - sr%dimen = tp%dimen - - ! - !sr%dsl = Huge(sr%dsl) ! set to huge positive values - !sr%il = -1 ! set to invalid indexes - ! - sr%overflow = .false. - - call search(tp%root) - - nfound = sr%nfound - - return - end function kdtree2_r_count - - function kdtree2_r_count_around_point(tp,idxin,correltime,r2) & - result(nfound) - ! Count the number of neighbors within square distance 'r2' around - ! point 'idxin' with decorrelation time 'correltime'. - ! - type (kdtree2), pointer :: tp - integer(pInt), intent (In) :: correltime, idxin - real(pReal), intent(in) :: r2 - integer(pInt) :: nfound - ! .. - ! .. - ! .. Intrinsic Functions .. - intrinsic HUGE - ! .. - allocate (sr%qv(tp%dimen)) - sr%qv = tp%the_data(:,idxin) - sr%ballsize = r2 - - sr%nn = 0 ! flag for fixed r search - sr%nfound = 0 - sr%centeridx = idxin - sr%correltime = correltime - nullify(sr%results) - - sr%nalloc = 0 ! we do not allocate any storage but that's OK - ! for counting. - - sr%ind => tp%ind - sr%rearrange = tp%rearrange - - if (sr%rearrange) then - sr%Data => tp%rearranged_data - else - sr%Data => tp%the_data - endif - sr%dimen = tp%dimen - - ! - !sr%dsl = Huge(sr%dsl) ! set to huge positive values - !sr%il = -1 ! set to invalid indexes - ! - sr%overflow = .false. - - call search(tp%root) - - nfound = sr%nfound - - return - end function kdtree2_r_count_around_point - - - subroutine validate_query_storage(n) - ! - ! make sure we have enough storage for n - ! - integer(pInt), intent(in) :: n - - if (size(sr%results,1) .lt. n) then - write (*,*) 'KD_TREE_TRANS: you did not provide enough storage for results(1:n)' - stop - return - endif - - return - end subroutine validate_query_storage - - function square_distance(d, iv,qv) result (res) - ! distance between iv[1:n] and qv[1:n] - ! .. Function Return Value .. - ! re-implemented to improve vectorization. - real(pReal) :: res - ! .. - ! .. - ! .. Scalar Arguments .. - integer(pInt) :: d - ! .. - ! .. Array Arguments .. - real(pReal) :: iv(:),qv(:) - ! .. - ! .. - res = sum( (iv(1:d)-qv(1:d))**2 ) - end function square_distance - - recursive subroutine search(node) - ! - ! This is the innermost core routine of the kd-tree search. Along - ! with "process_terminal_node", it is the performance bottleneck. - ! - ! This version uses a logically complete secondary search of - ! "box in bounds", whether the sear - ! - type (Tree_node), pointer :: node - ! .. - type(tree_node),pointer :: ncloser, nfarther - ! - integer(pInt) :: cut_dim, i - ! .. - real(pReal) :: qval, dis - real(pReal) :: ballsize - real(pReal), pointer :: qv(:) - type(interval), pointer :: box(:) - - if ((associated(node%left) .and. associated(node%right)) .eqv. .false.) then - ! we are on a terminal node - if (sr%nn .eq. 0) then - call process_terminal_node_fixedball(node) - else - call process_terminal_node(node) - endif - else - ! we are not on a terminal node - qv => sr%qv(1:) - cut_dim = node%cut_dim - qval = qv(cut_dim) - - if (qval < node%cut_val) then - ncloser => node%left - nfarther => node%right - dis = (node%cut_val_right - qval)**2 -! extra = node%cut_val - qval - else - ncloser => node%right - nfarther => node%left - dis = (node%cut_val_left - qval)**2 -! extra = qval- node%cut_val_left - endif - - if (associated(ncloser)) call search(ncloser) - - ! we may need to search the second node. - if (associated(nfarther)) then - ballsize = sr%ballsize -! dis=extra**2 - if (dis <= ballsize) then - ! - ! we do this separately as going on the first cut dimen is often - ! a good idea. - ! note that if extra**2 < sr%ballsize, then the next - ! check will also be false. - ! - box => node%box(1:) - do i=1,sr%dimen - if (i .ne. cut_dim) then - dis = dis + dis2_from_bnd(qv(i),box(i)%lower,box(i)%upper) - if (dis > ballsize) then - return - endif - endif - end do - - ! - ! if we are still here then we need to search mroe. - ! - call search(nfarther) - endif - endif - end if - end subroutine search - - - real(pReal) function dis2_from_bnd(x,amin,amax) result (res) - real(pReal), intent(in) :: x, amin,amax - - if (x > amax) then - res = (x-amax)**2; - return - else - if (x < amin) then - res = (amin-x)**2; - return - else - res = 0.0 - return - endif - endif - return - end function dis2_from_bnd - - logical function box_in_search_range(node, sr) result(res) - ! - ! Return the distance from 'qv' to the CLOSEST corner of node's - ! bounding box - ! for all coordinates outside the box. Coordinates inside the box - ! contribute nothing to the distance. - ! - type (tree_node), pointer :: node - type (tree_search_record), pointer :: sr - - integer(pInt) :: dimen, i - real(pReal) :: dis, ballsize - real(pReal) :: l, u - - dimen = sr%dimen - ballsize = sr%ballsize - dis = 0.0 - res = .true. - do i=1,dimen - l = node%box(i)%lower - u = node%box(i)%upper - dis = dis + (dis2_from_bnd(sr%qv(i),l,u)) - if (dis > ballsize) then - res = .false. - return - endif - end do - res = .true. - return - end function box_in_search_range - - - subroutine process_terminal_node(node) - ! - ! Look for actual near neighbors in 'node', and update - ! the search results on the sr data structure. - ! - type (tree_node), pointer :: node - ! - real(pReal), pointer :: qv(:) - integer(pInt), pointer :: ind(:) - real(pReal), pointer :: data(:,:) - ! - integer(pInt) :: dimen, i, indexofi, k, centeridx, correltime - real(pReal) :: ballsize, sd, newpri - logical :: rearrange - type(pq), pointer :: pqp - ! - ! copy values from sr to local variables - ! - ! - ! Notice, making local pointers with an EXPLICIT lower bound - ! seems to generate faster code. - ! why? I don't know. - qv => sr%qv(1:) - pqp => sr%pq - dimen = sr%dimen - ballsize = sr%ballsize - rearrange = sr%rearrange - ind => sr%ind(1:) - data => sr%Data(1:,1:) - centeridx = sr%centeridx - correltime = sr%correltime - - ! doing_correl = (centeridx >= 0) ! Do we have a decorrelation window? - ! include_point = .true. ! by default include all points - ! search through terminal bucket. - - mainloop: do i = node%l, node%u - if (rearrange) then - sd = 0.0 - do k = 1,dimen - sd = sd + (data(k,i) - qv(k))**2 - if (sd>ballsize) cycle mainloop - end do - indexofi = ind(i) ! only read it if we have not broken out - else - indexofi = ind(i) - sd = 0.0 - do k = 1,dimen - sd = sd + (data(k,indexofi) - qv(k))**2 - if (sd>ballsize) cycle mainloop - end do - endif - - if (centeridx > 0) then ! doing correlation interval? - if (abs(indexofi-centeridx) < correltime) cycle mainloop - endif - - - ! - ! two choices for any point. The list so far is either undersized, - ! or it is not. - ! - ! If it is undersized, then add the point and its distance - ! unconditionally. If the point added fills up the working - ! list then set the sr%ballsize, maximum distance bound (largest distance on - ! list) to be that distance, instead of the initialized +infinity. - ! - ! If the running list is full size, then compute the - ! distance but break out immediately if it is larger - ! than sr%ballsize, "best squared distance" (of the largest element), - ! as it cannot be a good neighbor. - ! - ! Once computed, compare to best_square distance. - ! if it is smaller, then delete the previous largest - ! element and add the new one. - - if (sr%nfound .lt. sr%nn) then - ! - ! add this point unconditionally to fill list. - ! - sr%nfound = sr%nfound +1 - newpri = pq_insert(pqp,sd,indexofi) - if (sr%nfound .eq. sr%nn) ballsize = newpri - ! we have just filled the working list. - ! put the best square distance to the maximum value - ! on the list, which is extractable from the PQ. - - else - ! - ! now, if we get here, - ! we know that the current node has a squared - ! distance smaller than the largest one on the list, and - ! belongs on the list. - ! Hence we replace that with the current one. - ! - ballsize = pq_replace_max(pqp,sd,indexofi) - endif - end do mainloop - ! - ! Reset sr variables which may have changed during loop - ! - sr%ballsize = ballsize - - end subroutine process_terminal_node - - subroutine process_terminal_node_fixedball(node) - ! - ! Look for actual near neighbors in 'node', and update - ! the search results on the sr data structure, i.e. - ! save all within a fixed ball. - ! - type (tree_node), pointer :: node - ! - real(pReal), pointer :: qv(:) - integer(pInt), pointer :: ind(:) - real(pReal), pointer :: data(:,:) - ! - integer(pInt) :: nfound - integer(pInt) :: dimen, i, indexofi, k - integer(pInt) :: centeridx, correltime, nn - real(pReal) :: ballsize, sd - logical :: rearrange - - ! - ! copy values from sr to local variables - ! - qv => sr%qv(1:) - dimen = sr%dimen - ballsize = sr%ballsize - rearrange = sr%rearrange - ind => sr%ind(1:) - data => sr%Data(1:,1:) - centeridx = sr%centeridx - correltime = sr%correltime - nn = sr%nn ! number to search for - nfound = sr%nfound - - ! search through terminal bucket. - mainloop: do i = node%l, node%u - - ! - ! two choices for any point. The list so far is either undersized, - ! or it is not. - ! - ! If it is undersized, then add the point and its distance - ! unconditionally. If the point added fills up the working - ! list then set the sr%ballsize, maximum distance bound (largest distance on - ! list) to be that distance, instead of the initialized +infinity. - ! - ! If the running list is full size, then compute the - ! distance but break out immediately if it is larger - ! than sr%ballsize, "best squared distance" (of the largest element), - ! as it cannot be a good neighbor. - ! - ! Once computed, compare to best_square distance. - ! if it is smaller, then delete the previous largest - ! element and add the new one. - - ! which index to the point do we use? - - if (rearrange) then - sd = 0.0 - do k = 1,dimen - sd = sd + (data(k,i) - qv(k))**2 - if (sd>ballsize) cycle mainloop - end do - indexofi = ind(i) ! only read it if we have not broken out - else - indexofi = ind(i) - sd = 0.0 - do k = 1,dimen - sd = sd + (data(k,indexofi) - qv(k))**2 - if (sd>ballsize) cycle mainloop - end do - endif - - if (centeridx > 0) then ! doing correlation interval? - if (abs(indexofi-centeridx) 1)then - ileft=ileft-1 - value=a(ileft); ivalue=ind(ileft) - else - value=a(iright); ivalue=ind(iright) - a(iright)=a(1); ind(iright)=ind(1) - iright=iright-1 - if (iright == 1) then - a(1)=value;ind(1)=ivalue - return - endif - endif - i=ileft - j=2*ileft - do while (j <= iright) - if(j < iright) then - if(a(j) < a(j+1)) j=j+1 - endif - if(value < a(j)) then - a(i)=a(j); ind(i)=ind(j) - i=j - j=j+j - else - j=iright+1 - endif - end do - a(i)=value; ind(i)=ivalue - end do - end subroutine heapsort - - subroutine heapsort_struct(a,n) - ! - ! Sort a(1:n) in ascending order - ! - ! - integer(pInt),intent(in) :: n - type(kdtree2_result),intent(inout) :: a(:) - - ! - ! - type(kdtree2_result) :: value ! temporary value - - integer(pInt) :: i,j - integer(pInt) :: ileft,iright - - ileft=n/2+1 - iright=n - - ! do i=1,n - ! ind(i)=i - ! Generate initial idum array - ! end do - - if(n.eq.1) return - - do - if(ileft > 1)then - ileft=ileft-1 - value=a(ileft) - else - value=a(iright) - a(iright)=a(1) - iright=iright-1 - if (iright == 1) then - a(1) = value - return - endif - endif - i=ileft - j=2*ileft - do while (j <= iright) - if(j < iright) then - if(a(j)%dis < a(j+1)%dis) j=j+1 - endif - if(value%dis < a(j)%dis) then - a(i)=a(j); - i=j - j=j+j - else - j=iright+1 - endif - end do - a(i)=value - end do - end subroutine heapsort_struct - -end module kdtree2_module - diff --git a/code/makefile b/code/makefile index ea00386ac..657692596 100644 --- a/code/makefile +++ b/code/makefile @@ -16,12 +16,12 @@ # make # make install # + specify in the "pathinfo:FFTW" where FFTW was installed. -# We essentially look for two library files "lib/libfftw3_threads.a" and "lib/libfftw3.a", so you can copy those, for instance, +# We essentially look for two library files "lib/libfftw3_threads" and "lib/libfftw3", so you can copy those, for instance, # into DAMASK_ROOT/lib/fftw/lib/ and specify "./fftw/" as pathinfo:FFTW # Use --enable-float in above configure for single precision... # Uses linux threads to parallelize fftw3 # -# Instead of the AMD Core Math Library a standard "liblapack.a/dylib/etc." can be used by leaving pathinfo:ACML blank +# Instead of the AMD Core Math Library a standard "liblapack.a/dylib/etc." can be used by leaving pathinfo:ACML and pathinfo:IKML blank ######################################################################################## # OPTIONS = standard (alternative): meaning #------------------------------------------------------------- @@ -58,10 +58,10 @@ DEBUG5 =-stand std03/std95 ######################################################################################## #auto values will be set by setup_code.py -FFTWROOT := /nethome/m.diehl/DAMASK/lib/fftw -IKMLROOT := +FFTWROOT := $(DAMASK_ROOT)/lib/fftw +IKMLROOT := ACMLROOT := /opt/acml4.4.0 -LAPACKROOT := +LAPACKROOT := F90 ?= ifort COMPILERNAME ?= $(F90) @@ -114,8 +114,8 @@ LIB_DIRS += -L$(ACMLROOT)/$(F90)64$(ACML_ARCH)/lib LIBRARIES += -lacml$(ACML_ARCH) else ifdef LAPACKROOT -LIB_DIRS += -L$(LAPACKROOT)/lib -L$(LAPACKROOT)/lib64 -LIBRARIES += -llapack +LIB_DIRS += -L$(LAPACKROOT)/lib64 -L$(LAPACKROOT)/lib +LIBRARIES += -llapack endif endif endif @@ -206,11 +206,8 @@ FEsolving.o: FEsolving.f90 basics.a -basics.a: kdtree2.o - ar rc basics.a kdtree2.o math.o debug.o numerics.o IO.o DAMASK_spectral_interface.o prec.o - -kdtree2.o: kdtree2.f90 math.o - $(PREFIX) $(COMPILERNAME) $(COMPILE) kdtree2.f90 $(SUFFIX) +basics.a: math.o + ar rc basics.a math.o debug.o numerics.o IO.o DAMASK_spectral_interface.o prec.o math.o: math.f90 debug.o $(PREFIX) $(COMPILERNAME) $(COMPILE) math.f90 $(SUFFIX) diff --git a/code/math.f90 b/code/math.f90 index c6a72fb86..6a0ce4baa 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -23,7 +23,9 @@ !############################################################## + use, intrinsic :: iso_c_binding use prec, only: pReal,pInt + use IO, only: IO_error implicit none real(pReal), parameter :: pi = 3.14159265358979323846264338327950288419716939937510_pReal @@ -126,6 +128,8 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = & 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal & /),(/4,36/)) + include 'fftw3.f03' + CONTAINS !************************************************************************** @@ -2884,6 +2888,7 @@ end subroutine ! Routine to calculate the mismatch between volume of reconstructed (compatible ! cube and determinant of defgrad at the FP + use debug, only: debug_verbosity implicit none ! input variables @@ -2898,9 +2903,11 @@ end subroutine integer(pInt) i,j,k real(pReal) vol_initial - print*, 'Calculating volume mismatch' - print '(a,/,e12.5,e12.5,e12.5)', ' Dimension:', geomdim - print '(a,/,i5,i5,i5)', ' Resolution:', res + if (debug_verbosity) then + print*, 'Calculating volume mismatch' + print '(a,/,e12.5,e12.5,e12.5)', ' Dimension:', geomdim + print '(a,/,i5,i5,i5)', ' Resolution:', res + endif vol_initial = geomdim(1)*geomdim(2)*geomdim(3)/(real(res(1)*res(2)*res(3), pReal)) do k = 1_pInt,res(3) @@ -2933,6 +2940,8 @@ subroutine shape_compare(res,geomdim,defgrad,nodes,centroids,shape_mismatch) ! Routine to calculate the mismatch between the vectors from the central point to ! the corners of reconstructed (combatible) volume element and the vectors calculated by deforming ! the initial volume element with the current deformation gradient + + use debug, only: debug_verbosity implicit none ! input variables @@ -2947,9 +2956,11 @@ subroutine shape_compare(res,geomdim,defgrad,nodes,centroids,shape_mismatch) real(pReal), dimension(8,3) :: coords_initial integer(pInt) i,j,k - print*, 'Calculating shape mismatch' - print '(a,/,e12.5,e12.5,e12.5)', ' Dimension:', geomdim - print '(a,/,i5,i5,i5)', ' Resolution:', res + if (debug_verbosity) then + print*, 'Calculating shape mismatch' + print '(a,/,e12.5,e12.5,e12.5)', ' Dimension:', geomdim + print '(a,/,i5,i5,i5)', ' Resolution:', res + endif coords_initial(1,1:3) = (/-geomdim(1)/2.0_pReal/real(res(1),pReal),& -geomdim(2)/2.0_pReal/real(res(2),pReal),& @@ -3007,6 +3018,7 @@ subroutine mesh_regular_grid(res,geomdim,defgrad_av,centroids,nodes) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Routine to build mesh of (distoreted) cubes for given coordinates (= center of the cubes) ! + use debug, only: debug_verbosity implicit none ! input variables integer(pInt), intent(in), dimension(3) :: res @@ -3032,9 +3044,11 @@ subroutine mesh_regular_grid(res,geomdim,defgrad_av,centroids,nodes) 0_pInt, 1_pInt, 1_pInt & /), & (/3,8/)) - print*, 'Meshing cubes around centroids' - print '(a,/,e12.5,e12.5,e12.5)', ' Dimension:', geomdim - print '(a,/,i5,i5,i5)', ' Resolution:', res + if (debug_verbosity) then + print*, 'Meshing cubes around centroids' + print '(a,/,e12.5,e12.5,e12.5)', ' Dimension:', geomdim + print '(a,/,i5,i5,i5)', ' Resolution:', res + endif nodes = 0.0_pReal wrappedCentroids = 0.0_pReal @@ -3074,6 +3088,7 @@ subroutine deformed_linear(res,geomdim,defgrad_av,defgrad,coord_avgCorner) ! Routine to calculate coordinates in current configuration for given defgrad ! using linear interpolation (blurres out high frequency defomation) ! + use debug, only: debug_verbosity implicit none ! input variables integer(pInt), intent(in), dimension(3) :: res @@ -3121,10 +3136,12 @@ subroutine deformed_linear(res,geomdim,defgrad_av,defgrad,coord_avgCorner) /), & (/3,6/)) - print*, 'Restore geometry using linear integration' - print '(a,/,e12.5,e12.5,e12.5)', ' Dimension:', geomdim - print '(a,/,i5,i5,i5)', ' Resolution:', res - + if (debug_verbosity) then + print*, 'Restore geometry using linear integration' + print '(a,/,e12.5,e12.5,e12.5)', ' Dimension:', geomdim + print '(a,/,i5,i5,i5)', ' Resolution:', res + endif + coord_avgOrder = 0.0_pReal do s = 0_pInt, 7_pInt ! corners (from 0 to 7) @@ -3185,6 +3202,7 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords) ! using integration in Fourier space (more accurate than deformed(...)) ! use numerics, only: fftw_timelimit, fftw_planner_flag + use debug, only: debug_verbosity implicit none ! input variables integer(pInt), intent(in), dimension(3) :: res @@ -3194,53 +3212,79 @@ subroutine deformed_fft(res,geomdim,defgrad_av,scaling,defgrad,coords) real(pReal), intent(in), dimension(res(1), res(2),res(3),3,3) :: defgrad ! output variables real(pReal), intent(out), dimension(res(1), res(2),res(3),3) :: coords - ! variables with dimension depending on input - complex(pReal), dimension(res(1)/2_pInt+1_pInt,res(2),res(3),3) :: coords_fft - complex(pReal), dimension(res(1), res(2),res(3),3,3) :: defgrad_fft +! allocatable arrays for fftw c routines + type(C_PTR) :: fftw_forth, fftw_back + type(C_PTR) :: coords_fftw, defgrad_fftw + real(pReal), dimension(:,:,:,:,:), pointer :: defgrad_real + complex(pReal), dimension(:,:,:,:,:), pointer :: defgrad_complex + real(pReal), dimension(:,:,:,:), pointer :: coords_real + complex(pReal), dimension(:,:,:,:), pointer :: coords_complex ! other variables - integer(pInt) :: i, j, k + integer(pInt) :: i, j, k, res1_red integer(pInt), dimension(3) :: k_s complex(pReal), parameter :: integration_factor = cmplx(0.0_pReal,pi*2.0_pReal) real(pReal), dimension(3) :: step, offset_coords - integer*8, dimension(2) :: plan_fft - print*, 'Restore geometry using FFT-based integration' - print '(a,/,e12.5,e12.5,e12.5)', ' Dimension:', geomdim - print '(a,/,i5,i5,i5)', ' Resolution:', res - - call dfftw_set_timelimit(fftw_timelimit) - call dfftw_plan_many_dft(plan_fft(1),3,(/res(1),res(2),res(3)/),9,& - defgrad_fft,(/res(1),res(2),res(3)/),1,res(1)*res(2)*res(3),& - defgrad_fft,(/res(1),res(2),res(3)/),1,res(1)*res(2)*res(3),-1,fftw_planner_flag) ! -1 = FFTW_FORWARD - call dfftw_plan_many_dft_c2r(plan_fft(2),3,(/res(1),res(2),res(3)/),3,& - coords_fft,(/res(1)/2_pInt+1_pInt,res(2),res(3)/),1,(res(1)/2_pInt+1_pInt)*res(2)*res(3),& - coords, (/res(1), res(2),res(3)/),1, res(1)* res(2)*res(3),fftw_planner_flag) - - coords_fft = 0.0 - defgrad_fft = defgrad ! cannot do memory efficient r2c transform as input field is destroyed during plan creation + if (debug_verbosity) then + print*, 'Restore geometry using FFT-based integration' + print '(a,/,e12.5,e12.5,e12.5)', ' Dimension:', geomdim + print '(a,/,i5,i5,i5)', ' Resolution:', res + endif + res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) step = geomdim/real(res, pReal) - call dfftw_execute_dft(plan_fft(1), defgrad_fft, defgrad_fft) - + if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=102) + call fftw_set_timelimit(fftw_timelimit) + defgrad_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*9_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8) + call c_f_pointer(defgrad_fftw, defgrad_real, [res(1)+2_pInt,res(2),res(3),3,3]) + call c_f_pointer(defgrad_fftw, defgrad_complex,[res1_red ,res(2),res(3),3,3]) + coords_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*3_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8) + call c_f_pointer(coords_fftw, coords_real, [res(1)+2_pInt,res(2),res(3),3]) + call c_f_pointer(coords_fftw, coords_complex, [res1_red ,res(2),res(3),3]) + + fftw_forth = fftw_plan_many_dft_r2c(3,(/res(3),res(2) ,res(1)/),9_pInt,& ! dimensions , length in each dimension in reversed order + defgrad_real,(/res(3),res(2) ,res(1)+2_pInt/),& ! input data , physical length in each dimension in reversed order + 1, res(3)*res(2)*(res(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions + defgrad_complex,(/res(3),res(2) ,res1_red/),& + 1, res(3)*res(2)* res1_red,fftw_planner_flag) + + fftw_back = fftw_plan_many_dft_c2r(3,(/res(3),res(2) ,res(1)/),3_pInt,& + coords_complex,(/res(3),res(2) ,res1_red/),& + 1, res(3)*res(2)* res1_red,& + coords_real,(/res(3),res(2) ,res(1)+2_pInt/),& + 1, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) + + + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + defgrad_real(i,j,k,1:3,1:3) = defgrad(i,j,k,1:3,1:3) ! ensure that data is aligned properly (fftw_alloc) + enddo; enddo; enddo + + call fftw_execute_dft_r2c(fftw_forth, defgrad_real, defgrad_complex) + + coords_complex = 0.0 do k = 1_pInt, res(3) k_s(3) = k-1_pInt if(k > res(3)/2_pInt+1_pInt) k_s(3) = k_s(3)-res(3) do j = 1_pInt, res(2) k_s(2) = j-1_pInt if(j > res(2)/2_pInt+1_pInt) k_s(2) = k_s(2)-res(2) - do i = 1_pInt, res(1)/2_pInt+1_pInt + do i = 1_pInt, res1_red k_s(1) = i-1_pInt - if(i/=1_pInt) coords_fft(i,j,k,1:3) = coords_fft(i,j,k,1:3)& - + defgrad_fft(i,j,k,1:3,1)*geomdim(1)/(real(k_s(1),pReal)*integration_factor) - if(j/=1_pInt) coords_fft(i,j,k,1:3) = coords_fft(i,j,k,1:3)& - + defgrad_fft(i,j,k,1:3,2)*geomdim(2)/(real(k_s(2),pReal)*integration_factor) - if(k/=1_pInt) coords_fft(i,j,k,1:3) = coords_fft(i,j,k,1:3)& - + defgrad_fft(i,j,k,1:3,3)*geomdim(3)/(real(k_s(3),pReal)*integration_factor) + if(i/=1_pInt) coords_complex(i,j,k,1:3) = coords_complex(i,j,k,1:3)& + + defgrad_complex(i,j,k,1:3,1)*geomdim(1)/(real(k_s(1),pReal)*integration_factor) + if(j/=1_pInt) coords_complex(i,j,k,1:3) = coords_complex(i,j,k,1:3)& + + defgrad_complex(i,j,k,1:3,2)*geomdim(2)/(real(k_s(2),pReal)*integration_factor) + if(k/=1_pInt) coords_complex(i,j,k,1:3) = coords_complex(i,j,k,1:3)& + + defgrad_complex(i,j,k,1:3,3)*geomdim(3)/(real(k_s(3),pReal)*integration_factor) enddo; enddo; enddo - call dfftw_execute_dft_c2r(plan_fft(2), coords_fft, coords) - coords = coords/real(res(1)*res(2)*res(3)) + call fftw_execute_dft_c2r(fftw_back,coords_complex,coords_real) + coords_real = coords_real/real(res(1)*res(2)*res(3)) + + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + coords(i,j,k,1:3) = coords_real(i,j,k,1:3) ! ensure that data is aligned properly (fftw_alloc) + enddo; enddo; enddo offset_coords = matmul(defgrad(1,1,1,1:3,1:3),step/2.0_pReal) - scaling*coords(1,1,1,1:3) do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) @@ -3255,12 +3299,13 @@ end subroutine deformed_fft !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -subroutine curl_fft(res,geomdim,vec_tens,field,curl_field) +subroutine curl_fft(res,geomdim,vec_tens,field,curl) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! calculates curl field using differentation in Fourier space ! use vec_tens to decide if tensor (3) or vector (1) use numerics, only: fftw_timelimit, fftw_planner_flag + use debug, only: debug_verbosity implicit none ! input variables integer(pInt), intent(in), dimension(3) :: res @@ -3268,31 +3313,54 @@ subroutine curl_fft(res,geomdim,vec_tens,field,curl_field) integer(pInt), intent(in) :: vec_tens real(pReal), intent(in), dimension(res(1), res(2),res(3),3,vec_tens) :: field ! output variables - real(pReal), intent(out), dimension(res(1), res(2),res(3),3,vec_tens) :: curl_field + real(pReal), intent(out), dimension(res(1), res(2),res(3),3,vec_tens) :: curl ! variables with dimension depending on input - complex(pReal), dimension(res(1), res(2),res(3),3,vec_tens) :: field_fft - complex(pReal), dimension(res(1)/2_pInt+1_pInt,res(2),res(3),3,vec_tens) :: curl_field_fft real(pReal), dimension(res(1)/2_pInt+1_pInt,res(2),res(3),3) :: xi + ! allocatable arrays for fftw c routines + type(C_PTR) :: fftw_forth, fftw_back + type(C_PTR) :: field_fftw, curl_fftw + real(pReal), dimension(:,:,:,:,:), pointer :: field_real + complex(pReal), dimension(:,:,:,:,:), pointer :: field_complex + real(pReal), dimension(:,:,:,:,:), pointer :: curl_real + complex(pReal), dimension(:,:,:,:,:), pointer :: curl_complex ! other variables - integer(pInt) i, j, k, l - integer*8 :: plan_fft(2) + integer(pInt) i, j, k, l, res1_red complex(pReal), parameter :: img =cmplx(0.0_pReal,1.0_pReal) - print*, 'Calculating curl of vector/tensor field' - print '(a,/,e12.5,e12.5,e12.5)', ' Dimension:', geomdim - print '(a,/,i5,i5,i5)', ' Resolution:', res + if (debug_verbosity) then + print*, 'Calculating curl of vector/tensor field' + print '(a,/,e12.5,e12.5,e12.5)', ' Dimension:', geomdim + print '(a,/,i5,i5,i5)', ' Resolution:', res + endif - call dfftw_set_timelimit(fftw_timelimit) - call dfftw_plan_many_dft(plan_fft(1),3,(/res(1),res(2),res(3)/),vec_tens*3_pInt,& - field_fft,(/res(1),res(2),res(3)/),1,res(1)*res(2)*res(3),& - field_fft,(/res(1),res(2),res(3)/),1,res(1)*res(2)*res(3),-1,fftw_planner_flag) ! -1 = FFTW_FORWARD - call dfftw_plan_many_dft_c2r(plan_fft(2),3,(/res(1),res(2),res(3)/),vec_tens*3_pInt,& - curl_field_fft,(/res(1)/2_pInt+1_pInt,res(2),res(3)/),1,(res(1)/2_pInt+1_pInt)*res(2)*res(3),& - curl_field,(/res(1),res(2),res(3)/),1,res(1)*res(2)*res(3),fftw_planner_flag) + res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) + if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=102) + call fftw_set_timelimit(fftw_timelimit) + field_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*vec_tens*3_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8) + call c_f_pointer(field_fftw, field_real, [res(1)+2_pInt,res(2),res(3),3,vec_tens]) + call c_f_pointer(field_fftw, field_complex,[res1_red ,res(2),res(3),3,vec_tens]) + curl_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*vec_tens*3_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8) + call c_f_pointer(curl_fftw, curl_real, [res(1)+2_pInt,res(2),res(3),3,vec_tens]) + call c_f_pointer(curl_fftw, curl_complex, [res1_red ,res(2),res(3),3,vec_tens]) - field_fft = field ! cannot do memory efficient r2c transform as input field is destroyed during plan creation + fftw_forth = fftw_plan_many_dft_r2c(3,(/res(3),res(2) ,res(1)/),vec_tens*3_pInt,& ! dimensions , length in each dimension in reversed order + field_real,(/res(3),res(2) ,res(1)+2_pInt/),& ! input data , physical length in each dimension in reversed order + 1, res(3)*res(2)*(res(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions + field_complex,(/res(3),res(2) ,res1_red/),& + 1, res(3)*res(2)* res1_red,fftw_planner_flag) - call dfftw_execute_dft_r2c(plan_fft(1), field_fft, field_fft) + fftw_back = fftw_plan_many_dft_c2r(3,(/res(3),res(2) ,res(1)/),vec_tens*3_pInt,& + curl_complex,(/res(3),res(2) ,res1_red/),& + 1, res(3)*res(2)* res1_red,& + curl_real,(/res(3),res(2) ,res(1)+2_pInt/),& + 1, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) + + + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + field_real(i,j,k,1:3,1:vec_tens) = field(i,j,k,1:3,1:vec_tens) ! ensure that data is aligned properly (fftw_alloc) + enddo; enddo; enddo + + call fftw_execute_dft_r2c(fftw_forth, field_real, field_complex) do k = 0_pInt, res(3)-1_pInt do j = 0_pInt, res(2)-1_pInt @@ -3305,30 +3373,34 @@ subroutine curl_fft(res,geomdim,vec_tens,field,curl_field) do k = 1, res(3) do j = 1, res(2) - do i = 1, res(1)/2+1 + do i = 1, res1_red do l = 1, vec_tens - curl_field_fft(i,j,k,1,l) =( field_fft(i,j,k,3,l)*xi(i,j,k,2) - field_fft(i,j,k,2,l)*xi(i,j,k,3))& + curl_complex(i,j,k,1,l) =( field_complex(i,j,k,3,l)*xi(i,j,k,2) - field_complex(i,j,k,2,l)*xi(i,j,k,3))& *img*pi*2.0_pReal - curl_field_fft(i,j,k,2,l) =(- field_fft(i,j,k,3,l)*xi(i,j,k,1) + field_fft(i,j,k,1,l)*xi(i,j,k,3))& + curl_complex(i,j,k,2,l) =(- field_complex(i,j,k,3,l)*xi(i,j,k,1) + field_complex(i,j,k,1,l)*xi(i,j,k,3))& *img*pi*2.0_pReal - curl_field_fft(i,j,k,3,l) =( field_fft(i,j,k,2,l)*xi(i,j,k,1) - field_fft(i,j,k,1,l)*xi(i,j,k,2))& + curl_complex(i,j,k,3,l) =( field_complex(i,j,k,2,l)*xi(i,j,k,1) - field_complex(i,j,k,1,l)*xi(i,j,k,2))& *img*pi*2.0_pReal enddo enddo; enddo; enddo - call dfftw_execute_dft_c2r(plan_fft(2), curl_field_fft, curl_field) + call fftw_execute_dft_c2r(fftw_back, curl_complex, curl_real) + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + curl(i,j,k,1:3,1:vec_tens) = curl_real(i,j,k,1:3,1:vec_tens) ! ensure that data is aligned properly (fftw_alloc) + enddo; enddo; enddo end subroutine curl_fft !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -subroutine divergence_fft(res,geomdim,vec_tens,field,divergence_field) +subroutine divergence_fft(res,geomdim,vec_tens,field,divergence) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! calculates divergence field using integration in Fourier space ! use vec_tens to decide if tensor (3) or vector (1) use numerics, only: fftw_timelimit, fftw_planner_flag + use debug, only: debug_verbosity implicit none ! input variables integer(pInt), intent(in), dimension(3) :: res @@ -3336,31 +3408,53 @@ subroutine divergence_fft(res,geomdim,vec_tens,field,divergence_field) integer(pInt), intent(in) :: vec_tens real(pReal), intent(in), dimension(res(1), res(2),res(3),vec_tens,3) :: field ! output variables - real(pReal), intent(out), dimension(res(1), res(2),res(3),vec_tens) :: divergence_field + real(pReal), intent(out), dimension(res(1), res(2),res(3),vec_tens) :: divergence ! variables with dimension depending on input - complex(pReal), dimension(res(1) ,res(2),res(3),vec_tens,3) :: field_fft - complex(pReal), dimension(res(1)/2_pInt+1_pInt,res(2),res(3),vec_tens) :: divergence_field_fft real(pReal), dimension(res(1)/2_pInt+1_pInt,res(2),res(3),3) :: xi + ! allocatable arrays for fftw c routines + type(C_PTR) :: fftw_forth, fftw_back + type(C_PTR) :: field_fftw, divergence_fftw + real(pReal), dimension(:,:,:,:,:), pointer :: field_real + complex(pReal), dimension(:,:,:,:,:), pointer :: field_complex + real(pReal), dimension(:,:,:,:), pointer :: divergence_real + complex(pReal), dimension(:,:,:,:), pointer :: divergence_complex ! other variables - integer(pInt) :: i, j, k + integer(pInt) :: i, j, k, res1_red complex(pReal), parameter :: img = cmplx(0.0_pReal,1.0_pReal) - integer*8, dimension(2) :: plan_fft - call dfftw_set_timelimit(fftw_timelimit) - call dfftw_plan_many_dft(plan_fft(1),3,(/res(1),res(2),res(3)/),vec_tens*3_pInt,& - field_fft,(/res(1),res(2),res(3)/),1,res(1)*res(2)*res(3),& - field_fft,(/res(1),res(2),res(3)/),1,res(1)*res(2)*res(3),-1,fftw_planner_flag) ! -1 = FFTW_FORWARD - call dfftw_plan_many_dft_c2r(plan_fft(2),3,(/res(1),res(2),res(3)/),vec_tens,& - divergence_field_fft,(/res(1)/2_pInt+1_pInt,res(2),res(3)/),1,(res(1)/2_pInt+1_pInt)*res(2)*res(3),& - divergence_field,(/res(1),res(2),res(3)/),1,res(1)*res(2)*res(3),fftw_planner_flag) + if (debug_verbosity) then + print '(a)', 'Calculating divergence of tensor/vector field using FFT' + print '(a,/,e12.5,e12.5,e12.5)', ' Dimension:', geomdim + print '(a,/,i5,i5,i5)', ' Resolution:', res + endif - print*, 'Calculating divergence of tensor/vector field using FFT' - print '(a,/,e12.5,e12.5,e12.5)', ' Dimension:', geomdim - print '(a,/,i5,i5,i5)', ' Resolution:', res + res1_red = res(1)/2_pInt + 1_pInt ! size of complex array in first dimension (c2r, r2c) + if (pReal /= C_DOUBLE .or. pInt /= C_INT) call IO_error(error_ID=102) + call fftw_set_timelimit(fftw_timelimit) + field_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*vec_tens*3_pInt,C_SIZE_T)) !C_SIZE_T is of type integer(8) + call c_f_pointer(field_fftw, field_real, [res(1)+2_pInt,res(2),res(3),3,vec_tens]) + call c_f_pointer(field_fftw, field_complex, [res1_red ,res(2),res(3),3,vec_tens]) + divergence_fftw = fftw_alloc_complex(int(res1_red *res(2)*res(3)*vec_tens,C_SIZE_T)) !C_SIZE_T is of type integer(8) + call c_f_pointer(divergence_fftw, divergence_real, [res(1)+2_pInt,res(2),res(3),vec_tens]) + call c_f_pointer(divergence_fftw, divergence_complex,[res1_red ,res(2),res(3),vec_tens]) - field_fft = field ! cannot do memory efficient r2c transform as input field is destroyed during plan creation + fftw_forth = fftw_plan_many_dft_r2c(3,(/res(3),res(2) ,res(1)/),vec_tens*3_pInt,& ! dimensions , length in each dimension in reversed order + field_real,(/res(3),res(2) ,res(1)+2_pInt/),& ! input data , physical length in each dimension in reversed order + 1, res(3)*res(2)*(res(1)+2_pInt),& ! striding , product of physical lenght in the 3 dimensions + field_complex,(/res(3),res(2) ,res1_red/),& + 1, res(3)*res(2)* res1_red,fftw_planner_flag) - call dfftw_execute_dft_r2c(plan_fft(1), field_fft, field_fft) + fftw_back = fftw_plan_many_dft_c2r(3,(/res(3),res(2) ,res(1)/),vec_tens,& + divergence_complex,(/res(3),res(2) ,res1_red/),& + 1, res(3)*res(2)* res1_red,& + divergence_real,(/res(3),res(2) ,res(1)+2_pInt/),& + 1, res(3)*res(2)*(res(1)+2_pInt),fftw_planner_flag) + + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + field_real(i,j,k,1:3,1:vec_tens) = field(i,j,k,1:3,1:vec_tens) ! ensure that data is aligned properly (fftw_alloc) + enddo; enddo; enddo + + call fftw_execute_dft_r2c(fftw_forth, field_real, field_complex) ! Alternative calculation of discrete frequencies k_s, ordered as in FFTW (wrap around) ! do k = 0,res(3)/2 -1 @@ -3394,26 +3488,32 @@ subroutine divergence_fft(res,geomdim,vec_tens,field,divergence_field) enddo; enddo; enddo do k = 1_pInt, res(3) do j = 1_pInt, res(2) - do i = 1_pInt, res(1)/2_pInt+1_pInt - divergence_field_fft(i,j,k,1) = sum(field_fft(i,j,k,1,1:3)*xi(i,j,k,1:3)) + do i = 1_pInt, res1_red + divergence_complex(i,j,k,1) = sum(field_complex(i,j,k,1:3,1)*xi(i,j,k,1:3)) !ToDo: check formula!!! if(vec_tens == 3_pInt) then - divergence_field_fft(i,j,k,2) = sum(field_fft(i,j,k,2,1:3)*xi(i,j,k,1:3)) - divergence_field_fft(i,j,k,3) = sum(field_fft(i,j,k,3,1:3)*xi(i,j,k,1:3)) + divergence_complex(i,j,k,2) = sum(field_complex(i,j,k,1:3,2)*xi(i,j,k,1:3)) + divergence_complex(i,j,k,3) = sum(field_complex(i,j,k,1:3,3)*xi(i,j,k,1:3)) endif enddo; enddo; enddo - divergence_field_fft = divergence_field_fft*img*2.0_pReal*pi + divergence_complex = divergence_complex*img*2.0_pReal*pi - call dfftw_execute_dft_c2r(plan_fft(2), divergence_field_fft, divergence_field) + call fftw_execute_dft_c2r(fftw_back, divergence_complex, divergence_real) + + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + divergence(i,j,k,1:vec_tens) = divergence_real(i,j,k,1:vec_tens) ! ensure that data is aligned properly (fftw_alloc) + enddo; enddo; enddo ! why not weighting the divergence field? + end subroutine divergence_fft !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - subroutine divergence_fdm(res,geomdim,vec_tens,order,field,divergence_field) + subroutine divergence_fdm(res,geomdim,vec_tens,order,field,divergence) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! calculates divergence field using FDM with variable accuracy ! use vec_tes to decide if tensor (3) or vector (1) + use debug, only: debug_verbosity implicit none integer(pInt), intent(in), dimension(3) :: res integer(pInt), intent(in) :: vec_tens @@ -3421,7 +3521,7 @@ subroutine divergence_fft(res,geomdim,vec_tens,field,divergence_field) real(pReal), intent(in), dimension(3) :: geomdim real(pReal), intent(in), dimension(res(1),res(2),res(3),vec_tens,3) :: field ! output variables - real(pReal), intent(out), dimension(res(1),res(2),res(3),vec_tens) :: divergence_field + real(pReal), intent(out), dimension(res(1),res(2),res(3),vec_tens) :: divergence ! other variables integer(pInt), dimension(6,3) :: coordinates integer(pInt) i, j, k, m, l @@ -3432,11 +3532,13 @@ subroutine divergence_fft(res,geomdim,vec_tens,field,divergence_field) 4.0_pReal/5.0_pReal,-1.0_pReal/ 5.0_pReal,4.0_pReal/105.0_pReal,-1.0_pReal/280.0_pReal/),& (/4,4/)) - print*, 'Calculating divergence of tensor/vector field using FDM' - print '(a,/,e12.5,e12.5,e12.5)', ' Dimension:', geomdim - print '(a,/,i5,i5,i5)', ' Resolution:', res + if (debug_verbosity) then + print*, 'Calculating divergence of tensor/vector field using FDM' + print '(a,/,e12.5,e12.5,e12.5)', ' Dimension:', geomdim + print '(a,/,i5,i5,i5)', ' Resolution:', res + endif - divergence_field = 0.0_pReal + divergence = 0.0_pReal order = order + 1_pInt do k = 0_pInt, res(3)-1_pInt; do j = 0_pInt, res(2)-1_pInt; do i = 0_pInt, res(1)-1_pInt do m = 1_pInt, order @@ -3453,7 +3555,7 @@ subroutine divergence_fft(res,geomdim,vec_tens,field,divergence_field) coordinates(6,1:3) = mesh_location(mesh_index((/i,j,k-m/),(/res(1),res(2),res(3)/)),(/res(1),res(2),res(3)/))& + (/1_pInt,1_pInt,1_pInt/) do l = 1_pInt, vec_tens - divergence_field(i+1_pInt,j+1_pInt,k+1_pInt,l) = divergence_field(i+1_pInt,j+1_pInt,k+1_pInt,l) + FDcoefficient(m,order) * & + divergence(i+1_pInt,j+1_pInt,k+1_pInt,l) = divergence(i+1_pInt,j+1_pInt,k+1_pInt,l) + FDcoefficient(m,order) * & ((field(coordinates(1,1),coordinates(1,2),coordinates(1,3),l,1)- & field(coordinates(2,1),coordinates(2,2),coordinates(2,3),l,1))*real(res(1),pReal)/geomdim(1) +& (field(coordinates(3,1),coordinates(3,2),coordinates(3,3),l,2)- & @@ -3576,3 +3678,1893 @@ subroutine calculate_cauchy(res,defgrad,p_stress,c_stress) end subroutine calculate_cauchy END MODULE math + +!############################################################################################################################# +! BEGIN KDTREE2 +!############################################################################################################################# +!(c) Matthew Kennel, Institute for Nonlinear Science (2004) +! +! Licensed under the Academic Free License version 1.1 found in file LICENSE +! with additional provisions found in that same file. +! +!####################################################### +! modifications: changed precision according to prec.f90 +! k.komerla, m.diehl +!####################################################### + +module kdtree2_priority_queue_module + use prec + ! + ! maintain a priority queue (PQ) of data, pairs of 'priority/payload', + ! implemented with a binary heap. This is the type, and the 'dis' field + ! is the priority. + ! + type kdtree2_result + ! a pair of distances, indexes + real(pReal) :: dis !=0.0 + integer(pInt) :: idx !=-1 Initializers cause some bugs in compilers. + end type kdtree2_result + ! + ! A heap-based priority queue lets one efficiently implement the following + ! operations, each in log(N) time, as opposed to linear time. + ! + ! 1) add a datum (push a datum onto the queue, increasing its length) + ! 2) return the priority value of the maximum priority element + ! 3) pop-off (and delete) the element with the maximum priority, decreasing + ! the size of the queue. + ! 4) replace the datum with the maximum priority with a supplied datum + ! (of either higher or lower priority), maintaining the size of the + ! queue. + ! + ! + ! In the k-d tree case, the 'priority' is the square distance of a point in + ! the data set to a reference point. The goal is to keep the smallest M + ! distances to a reference point. The tree algorithm searches terminal + ! nodes to decide whether to add points under consideration. + ! + ! A priority queue is useful here because it lets one quickly return the + ! largest distance currently existing in the list. If a new candidate + ! distance is smaller than this, then the new candidate ought to replace + ! the old candidate. In priority queue terms, this means removing the + ! highest priority element, and inserting the new one. + ! + ! Algorithms based on Cormen, Leiserson, Rivest, _Introduction + ! to Algorithms_, 1990, with further optimization by the author. + ! + ! Originally informed by a C implementation by Sriranga Veeraraghavan. + ! + ! This module is not written in the most clear way, but is implemented such + ! for speed, as it its operations will be called many times during searches + ! of large numbers of neighbors. + ! + type pq + ! + ! The priority queue consists of elements + ! priority(1:heap_size), with associated payload(:). + ! + ! There are heap_size active elements. + ! Assumes the allocation is always sufficient. Will NOT increase it + ! to match. + integer(pInt) :: heap_size = 0 + type(kdtree2_result), pointer :: elems(:) + end type pq + + public :: kdtree2_result + + public :: pq + public :: pq_create + public :: pq_delete, pq_insert + public :: pq_extract_max, pq_max, pq_replace_max, pq_maxpri + private + +contains + + + function pq_create(results_in) result(res) + ! + ! Create a priority queue from ALREADY allocated + ! array pointers for storage. NOTE! It will NOT + ! add any alements to the heap, i.e. any existing + ! data in the input arrays will NOT be used and may + ! be overwritten. + ! + ! usage: + ! real(pReal), pointer :: x(:) + ! integer(pInt), pointer :: k(:) + ! allocate(x(1000),k(1000)) + ! pq => pq_create(x,k) + ! + type(kdtree2_result), target:: results_in(:) + type(pq) :: res + ! + ! + integer(pInt) :: nalloc + + nalloc = size(results_in,1) + if (nalloc .lt. 1) then + write (*,*) 'PQ_CREATE: error, input arrays must be allocated.' + end if + res%elems => results_in + res%heap_size = 0 + return + end function pq_create + + ! + ! operations for getting parents and left + right children + ! of elements in a binary heap. + ! + +! +! These are written inline for speed. +! +! integer(pInt) function parent(i) +! integer(pInt), intent(in) :: i +! parent = (i/2) +! return +! end function parent + +! integer(pInt) function left(i) +! integer(pInt), intent(in) ::i +! left = (2*i) +! return +! end function left + +! integer(pInt) function right(i) +! integer(pInt), intent(in) :: i +! right = (2*i)+1 +! return +! end function right + +! logical function compare_priority(p1,p2) +! real(pReal), intent(in) :: p1, p2 +! +! compare_priority = (p1 .gt. p2) +! return +! end function compare_priority + + subroutine heapify(a,i_in) + ! + ! take a heap rooted at 'i' and force it to be in the + ! heap canonical form. This is performance critical + ! and has been tweaked a little to reflect this. + ! + type(pq),pointer :: a + integer(pInt), intent(in) :: i_in + ! + integer(pInt) :: i, l, r, largest + + real(pReal) :: pri_i, pri_l, pri_r, pri_largest + + + type(kdtree2_result) :: temp + + i = i_in + +bigloop: do + l = 2*i ! left(i) + r = l+1 ! right(i) + ! + ! set 'largest' to the index of either i, l, r + ! depending on whose priority is largest. + ! + ! note that l or r can be larger than the heap size + ! in which case they do not count. + + + ! does left child have higher priority? + if (l .gt. a%heap_size) then + ! we know that i is the largest as both l and r are invalid. + exit + else + pri_i = a%elems(i)%dis + pri_l = a%elems(l)%dis + if (pri_l .gt. pri_i) then + largest = l + pri_largest = pri_l + else + largest = i + pri_largest = pri_i + endif + + ! + ! between i and l we have a winner + ! now choose between that and r. + ! + if (r .le. a%heap_size) then + pri_r = a%elems(r)%dis + if (pri_r .gt. pri_largest) then + largest = r + endif + endif + endif + + if (largest .ne. i) then + ! swap data in nodes largest and i, then heapify + + temp = a%elems(i) + a%elems(i) = a%elems(largest) + a%elems(largest) = temp + ! + ! Canonical heapify() algorithm has tail-ecursive call: + ! + ! call heapify(a,largest) + ! we will simulate with cycle + ! + i = largest + cycle bigloop ! continue the loop + else + return ! break from the loop + end if + enddo bigloop + return + end subroutine heapify + + subroutine pq_max(a,e) + ! + ! return the priority and its payload of the maximum priority element + ! on the queue, which should be the first one, if it is + ! in heapified form. + ! + type(pq),pointer :: a + type(kdtree2_result),intent(out) :: e + + if (a%heap_size .gt. 0) then + e = a%elems(1) + else + write (*,*) 'PQ_MAX: ERROR, heap_size < 1' + stop + endif + return + end subroutine pq_max + + real(pReal) function pq_maxpri(a) + type(pq), pointer :: a + + if (a%heap_size .gt. 0) then + pq_maxpri = a%elems(1)%dis + else + write (*,*) 'PQ_MAX_PRI: ERROR, heapsize < 1' + stop + endif + return + end function pq_maxpri + + subroutine pq_extract_max(a,e) + ! + ! return the priority and payload of maximum priority + ! element, and remove it from the queue. + ! (equivalent to 'pop()' on a stack) + ! + type(pq),pointer :: a + type(kdtree2_result), intent(out) :: e + + if (a%heap_size .ge. 1) then + ! + ! return max as first element + ! + e = a%elems(1) + + ! + ! move last element to first + ! + a%elems(1) = a%elems(a%heap_size) + a%heap_size = a%heap_size-1 + call heapify(a,1) + return + else + write (*,*) 'PQ_EXTRACT_MAX: error, attempted to pop non-positive PQ' + stop + end if + + end subroutine pq_extract_max + + + real(pReal) function pq_insert(a,dis,idx) + ! + ! Insert a new element and return the new maximum priority, + ! which may or may not be the same as the old maximum priority. + ! + type(pq),pointer :: a + real(pReal), intent(in) :: dis + integer(pInt), intent(in) :: idx + ! type(kdtree2_result), intent(in) :: e + ! + integer(pInt) :: i, isparent + real(pReal) :: parentdis + ! + + ! if (a%heap_size .ge. a%max_elems) then + ! write (*,*) 'PQ_INSERT: error, attempt made to insert element on full PQ' + ! stop + ! else + a%heap_size = a%heap_size + 1 + i = a%heap_size + + do while (i .gt. 1) + isparent = int(i/2) + parentdis = a%elems(isparent)%dis + if (dis .gt. parentdis) then + ! move what was in i's parent into i. + a%elems(i)%dis = parentdis + a%elems(i)%idx = a%elems(isparent)%idx + i = isparent + else + exit + endif + end do + + ! insert the element at the determined position + a%elems(i)%dis = dis + a%elems(i)%idx = idx + + pq_insert = a%elems(1)%dis + return + ! end if + + end function pq_insert + + subroutine pq_adjust_heap(a,i) + type(pq),pointer :: a + integer(pInt), intent(in) :: i + ! + ! nominally arguments (a,i), but specialize for a=1 + ! + ! This routine assumes that the trees with roots 2 and 3 are already heaps, i.e. + ! the children of '1' are heaps. When the procedure is completed, the + ! tree rooted at 1 is a heap. + real(pReal) :: prichild + integer(pInt) :: parent, child, N + + type(kdtree2_result) :: e + + e = a%elems(i) + + parent = i + child = 2*i + N = a%heap_size + + do while (child .le. N) + if (child .lt. N) then + if (a%elems(child)%dis .lt. a%elems(child+1)%dis) then + child = child+1 + endif + endif + prichild = a%elems(child)%dis + if (e%dis .ge. prichild) then + exit + else + ! move child into parent. + a%elems(parent) = a%elems(child) + parent = child + child = 2*parent + end if + end do + a%elems(parent) = e + return + end subroutine pq_adjust_heap + + + real(pReal) function pq_replace_max(a,dis,idx) + ! + ! Replace the extant maximum priority element + ! in the PQ with (dis,idx). Return + ! the new maximum priority, which may be larger + ! or smaller than the old one. + ! + type(pq),pointer :: a + real(pReal), intent(in) :: dis + integer(pInt), intent(in) :: idx +! type(kdtree2_result), intent(in) :: e + ! not tested as well! + + integer(pInt) :: parent, child, N + real(pReal) :: prichild, prichildp1 + + type(kdtree2_result) :: etmp + + if (.true.) then + N=a%heap_size + if (N .ge. 1) then + parent =1 + child=2 + + loop: do while (child .le. N) + prichild = a%elems(child)%dis + + ! + ! posibly child+1 has higher priority, and if + ! so, get it, and increment child. + ! + + if (child .lt. N) then + prichildp1 = a%elems(child+1)%dis + if (prichild .lt. prichildp1) then + child = child+1 + prichild = prichildp1 + endif + endif + + if (dis .ge. prichild) then + exit loop + ! we have a proper place for our new element, + ! bigger than either children's priority. + else + ! move child into parent. + a%elems(parent) = a%elems(child) + parent = child + child = 2*parent + end if + end do loop + a%elems(parent)%dis = dis + a%elems(parent)%idx = idx + pq_replace_max = a%elems(1)%dis + else + a%elems(1)%dis = dis + a%elems(1)%idx = idx + pq_replace_max = dis + endif + else + ! + ! slower version using elementary pop and push operations. + ! + call pq_extract_max(a,etmp) + etmp%dis = dis + etmp%idx = idx + pq_replace_max = pq_insert(a,dis,idx) + endif + return + end function pq_replace_max + + subroutine pq_delete(a,i) + ! + ! delete item with index 'i' + ! + type(pq),pointer :: a + integer(pInt) :: i + + if ((i .lt. 1) .or. (i .gt. a%heap_size)) then + write (*,*) 'PQ_DELETE: error, attempt to remove out of bounds element.' + stop + endif + + ! swap the item to be deleted with the last element + ! and shorten heap by one. + a%elems(i) = a%elems(a%heap_size) + a%heap_size = a%heap_size - 1 + + call heapify(a,i) + + end subroutine pq_delete + +end module kdtree2_priority_queue_module + + +module kdtree2_module + use prec + use kdtree2_priority_queue_module + ! K-D tree routines in Fortran 90 by Matt Kennel. + ! Original program was written in Sather by Steve Omohundro and + ! Matt Kennel. Only the Euclidean metric is supported. + ! + ! + ! This module is identical to 'kd_tree', except that the order + ! of subscripts is reversed in the data file. + ! In otherwords for an embedding of N D-dimensional vectors, the + ! data file is here, in natural Fortran order data(1:D, 1:N) + ! because Fortran lays out columns first, + ! + ! whereas conventionally (C-style) it is data(1:N,1:D) + ! as in the original kd_tree module. + ! + !-------------DATA TYPE, CREATION, DELETION--------------------- + public :: pReal + public :: kdtree2, kdtree2_result, tree_node, kdtree2_create, kdtree2_destroy + !--------------------------------------------------------------- + !-------------------SEARCH ROUTINES----------------------------- + public :: kdtree2_n_nearest,kdtree2_n_nearest_around_point + ! Return fixed number of nearest neighbors around arbitrary vector, + ! or extant point in dataset, with decorrelation window. + ! + public :: kdtree2_r_nearest, kdtree2_r_nearest_around_point + ! Return points within a fixed ball of arb vector/extant point + ! + public :: kdtree2_sort_results + ! Sort, in order of increasing distance, rseults from above. + ! + public :: kdtree2_r_count, kdtree2_r_count_around_point + ! Count points within a fixed ball of arb vector/extant point + ! + public :: kdtree2_n_nearest_brute_force, kdtree2_r_nearest_brute_force + ! brute force of kdtree2_[n|r]_nearest + !---------------------------------------------------------------- + + + integer(pInt), parameter :: bucket_size = 12 + ! The maximum number of points to keep in a terminal node. + + type interval + real(pReal) :: lower,upper + end type interval + + type :: tree_node + ! an internal tree node + private + integer(pInt) :: cut_dim + ! the dimension to cut + real(pReal) :: cut_val + ! where to cut the dimension + real(pReal) :: cut_val_left, cut_val_right + ! improved cutoffs knowing the spread in child boxes. + integer(pInt) :: l, u + type (tree_node), pointer :: left, right + type(interval), pointer :: box(:) => null() + ! child pointers + ! Points included in this node are indexes[k] with k \in [l,u] + + + end type tree_node + + type :: kdtree2 + ! Global information about the tree, one per tree + integer(pInt) :: dimen=0, n=0 + ! dimensionality and total # of points + real(pReal), pointer :: the_data(:,:) => null() + ! pointer to the actual data array + ! + ! IMPORTANT NOTE: IT IS DIMENSIONED the_data(1:d,1:N) + ! which may be opposite of what may be conventional. + ! This is, because in Fortran, the memory layout is such that + ! the first dimension is in sequential order. Hence, with + ! (1:d,1:N), all components of the vector will be in consecutive + ! memory locations. The search time is dominated by the + ! evaluation of distances in the terminal nodes. Putting all + ! vector components in consecutive memory location improves + ! memory cache locality, and hence search speed, and may enable + ! vectorization on some processors and compilers. + + integer(pInt), pointer :: ind(:) => null() + ! permuted index into the data, so that indexes[l..u] of some + ! bucket represent the indexes of the actual points in that + ! bucket. + logical :: sort = .false. + ! do we always sort output results? + logical :: rearrange = .false. + real(pReal), pointer :: rearranged_data(:,:) => null() + ! if (rearrange .eqv. .true.) then rearranged_data has been + ! created so that rearranged_data(:,i) = the_data(:,ind(i)), + ! permitting search to use more cache-friendly rearranged_data, at + ! some initial computation and storage cost. + type (tree_node), pointer :: root => null() + ! root pointer of the tree + end type kdtree2 + + + type :: tree_search_record + ! + ! One of these is created for each search. + ! + private + ! + ! Many fields are copied from the tree structure, in order to + ! speed up the search. + ! + integer(pInt) :: dimen + integer(pInt) :: nn, nfound + real(pReal) :: ballsize + integer(pInt) :: centeridx=999, correltime=9999 + ! exclude points within 'correltime' of 'centeridx', iff centeridx >= 0 + integer(pInt) :: nalloc ! how much allocated for results(:)? + logical :: rearrange ! are the data rearranged or original? + ! did the # of points found overflow the storage provided? + logical :: overflow + real(pReal), pointer :: qv(:) ! query vector + type(kdtree2_result), pointer :: results(:) ! results + type(pq) :: pq + real(pReal), pointer :: data(:,:) ! temp pointer to data + integer(pInt), pointer :: ind(:) ! temp pointer to indexes + end type tree_search_record + + private + ! everything else is private. + + type(tree_search_record), save, target :: sr ! A GLOBAL VARIABLE for search + +contains + + function kdtree2_create(input_data,dim,sort,rearrange) result (mr) + ! + ! create the actual tree structure, given an input array of data. + ! + ! Note, input data is input_data(1:d,1:N), NOT the other way around. + ! THIS IS THE REVERSE OF THE PREVIOUS VERSION OF THIS MODULE. + ! The reason for it is cache friendliness, improving performance. + ! + ! Optional arguments: If 'dim' is specified, then the tree + ! will only search the first 'dim' components + ! of input_data, otherwise, dim is inferred + ! from SIZE(input_data,1). + ! + ! if sort .eqv. .true. then output results + ! will be sorted by increasing distance. + ! default=.false., as it is faster to not sort. + ! + ! if rearrange .eqv. .true. then an internal + ! copy of the data, rearranged by terminal node, + ! will be made for cache friendliness. + ! default=.true., as it speeds searches, but + ! building takes longer, and extra memory is used. + ! + ! .. Function Return Cut_value .. + type (kdtree2), pointer :: mr + integer(pInt), intent(in), optional :: dim + logical, intent(in), optional :: sort + logical, intent(in), optional :: rearrange + ! .. + ! .. Array Arguments .. + real(pReal), target :: input_data(:,:) + ! + integer(pInt) :: i + ! .. + allocate (mr) + mr%the_data => input_data + ! pointer assignment + + if (present(dim)) then + mr%dimen = dim + else + mr%dimen = size(input_data,1) + end if + mr%n = size(input_data,2) + + if (mr%dimen > mr%n) then + ! unlikely to be correct + write (*,*) 'KD_TREE_TRANS: likely user error.' + write (*,*) 'KD_TREE_TRANS: You passed in matrix with D=',mr%dimen + write (*,*) 'KD_TREE_TRANS: and N=',mr%n + write (*,*) 'KD_TREE_TRANS: note, that new format is data(1:D,1:N)' + write (*,*) 'KD_TREE_TRANS: with usually N >> D. If N =approx= D, then a k-d tree' + write (*,*) 'KD_TREE_TRANS: is not an appropriate data structure.' + stop + end if + + call build_tree(mr) + + if (present(sort)) then + mr%sort = sort + else + mr%sort = .false. + endif + + if (present(rearrange)) then + mr%rearrange = rearrange + else + mr%rearrange = .true. + endif + + if (mr%rearrange) then + allocate(mr%rearranged_data(mr%dimen,mr%n)) + do i=1,mr%n + mr%rearranged_data(:,i) = mr%the_data(:, & + mr%ind(i)) + enddo + else + nullify(mr%rearranged_data) + endif + + end function kdtree2_create + + subroutine build_tree(tp) + type (kdtree2), pointer :: tp + ! .. + integer(pInt) :: j + type(tree_node), pointer :: dummy => null() + ! .. + allocate (tp%ind(tp%n)) + forall (j=1:tp%n) + tp%ind(j) = j + end forall + tp%root => build_tree_for_range(tp,1,tp%n, dummy) + end subroutine build_tree + + recursive function build_tree_for_range(tp,l,u,parent) result (res) + ! .. Function Return Cut_value .. + type (tree_node), pointer :: res + ! .. + ! .. Structure Arguments .. + type (kdtree2), pointer :: tp + type (tree_node),pointer :: parent + ! .. + ! .. Scalar Arguments .. + integer(pInt), intent (In) :: l, u + ! .. + ! .. Local Scalars .. + integer(pInt) :: i, c, m, dimen + logical :: recompute + real(pReal) :: average + +!!$ If (.False.) Then +!!$ If ((l .Lt. 1) .Or. (l .Gt. tp%n)) Then +!!$ Stop 'illegal L value in build_tree_for_range' +!!$ End If +!!$ If ((u .Lt. 1) .Or. (u .Gt. tp%n)) Then +!!$ Stop 'illegal u value in build_tree_for_range' +!!$ End If +!!$ If (u .Lt. l) Then +!!$ Stop 'U is less than L, thats illegal.' +!!$ End If +!!$ Endif +!!$ + ! first compute min and max + dimen = tp%dimen + allocate (res) + allocate(res%box(dimen)) + + ! First, compute an APPROXIMATE bounding box of all points associated with this node. + if ( u < l ) then + ! no points in this box + nullify(res) + return + end if + + if ((u-l)<=bucket_size) then + ! + ! always compute true bounding box for terminal nodes. + ! + do i=1,dimen + call spread_in_coordinate(tp,i,l,u,res%box(i)) + end do + res%cut_dim = 0 + res%cut_val = 0.0 + res%l = l + res%u = u + res%left =>null() + res%right => null() + else + ! + ! modify approximate bounding box. This will be an + ! overestimate of the true bounding box, as we are only recomputing + ! the bounding box for the dimension that the parent split on. + ! + ! Going to a true bounding box computation would significantly + ! increase the time necessary to build the tree, and usually + ! has only a very small difference. This box is not used + ! for searching but only for deciding which coordinate to split on. + ! + do i=1,dimen + recompute=.true. + if (associated(parent)) then + if (i .ne. parent%cut_dim) then + recompute=.false. + end if + endif + if (recompute) then + call spread_in_coordinate(tp,i,l,u,res%box(i)) + else + res%box(i) = parent%box(i) + endif + end do + + + c = maxloc(res%box(1:dimen)%upper-res%box(1:dimen)%lower,1) + ! + ! c is the identity of which coordinate has the greatest spread. + ! + + if (.false.) then + ! select exact median to have fully balanced tree. + m = (l+u)/2 + call select_on_coordinate(tp%the_data,tp%ind,c,m,l,u) + else + ! + ! select point halfway between min and max, as per A. Moore, + ! who says this helps in some degenerate cases, or + ! actual arithmetic average. + ! + if (.true.) then + ! actually compute average + average = sum(tp%the_data(c,tp%ind(l:u))) / real(u-l+1,pReal) + else + average = (res%box(c)%upper + res%box(c)%lower)/2.0 + endif + + res%cut_val = average + m = select_on_coordinate_value(tp%the_data,tp%ind,c,average,l,u) + endif + + ! moves indexes around + res%cut_dim = c + res%l = l + res%u = u +! res%cut_val = tp%the_data(c,tp%ind(m)) + + res%left => build_tree_for_range(tp,l,m,res) + res%right => build_tree_for_range(tp,m+1,u,res) + + if (associated(res%right) .eqv. .false.) then + res%box = res%left%box + res%cut_val_left = res%left%box(c)%upper + res%cut_val = res%cut_val_left + elseif (associated(res%left) .eqv. .false.) then + res%box = res%right%box + res%cut_val_right = res%right%box(c)%lower + res%cut_val = res%cut_val_right + else + res%cut_val_right = res%right%box(c)%lower + res%cut_val_left = res%left%box(c)%upper + res%cut_val = (res%cut_val_left + res%cut_val_right)/2 + + + ! now remake the true bounding box for self. + ! Since we are taking unions (in effect) of a tree structure, + ! this is much faster than doing an exhaustive + ! search over all points + res%box%upper = max(res%left%box%upper,res%right%box%upper) + res%box%lower = min(res%left%box%lower,res%right%box%lower) + endif + end if + end function build_tree_for_range + + integer(pInt) function select_on_coordinate_value(v,ind,c,alpha,li,ui) & + result(res) + ! Move elts of ind around between l and u, so that all points + ! <= than alpha (in c cooordinate) are first, and then + ! all points > alpha are second. + + ! + ! Algorithm (matt kennel). + ! + ! Consider the list as having three parts: on the left, + ! the points known to be <= alpha. On the right, the points + ! known to be > alpha, and in the middle, the currently unknown + ! points. The algorithm is to scan the unknown points, starting + ! from the left, and swapping them so that they are added to + ! the left stack or the right stack, as appropriate. + ! + ! The algorithm finishes when the unknown stack is empty. + ! + ! .. Scalar Arguments .. + integer(pInt), intent (In) :: c, li, ui + real(pReal), intent(in) :: alpha + ! .. + real(pReal) :: v(1:,1:) + integer(pInt) :: ind(1:) + integer(pInt) :: tmp + ! .. + integer(pInt) :: lb, rb + ! + ! The points known to be <= alpha are in + ! [l,lb-1] + ! + ! The points known to be > alpha are in + ! [rb+1,u]. + ! + ! Therefore we add new points into lb or + ! rb as appropriate. When lb=rb + ! we are done. We return the location of the last point <= alpha. + ! + ! + lb = li; rb = ui + + do while (lb < rb) + if ( v(c,ind(lb)) <= alpha ) then + ! it is good where it is. + lb = lb+1 + else + ! swap it with rb. + tmp = ind(lb); ind(lb) = ind(rb); ind(rb) = tmp + rb = rb-1 + endif + end do + + ! now lb .eq. ub + if (v(c,ind(lb)) <= alpha) then + res = lb + else + res = lb-1 + endif + + end function select_on_coordinate_value + + subroutine select_on_coordinate(v,ind,c,k,li,ui) + ! Move elts of ind around between l and u, so that the kth + ! element + ! is >= those below, <= those above, in the coordinate c. + ! .. Scalar Arguments .. + integer(pInt), intent (In) :: c, k, li, ui + ! .. + integer(pInt) :: i, l, m, s, t, u + ! .. + real(pReal) :: v(:,:) + integer(pInt) :: ind(:) + ! .. + l = li + u = ui + do while (l=k) u = m - 1 + end do + end subroutine select_on_coordinate + + subroutine spread_in_coordinate(tp,c,l,u,interv) + ! the spread in coordinate 'c', between l and u. + ! + ! Return lower bound in 'smin', and upper in 'smax', + ! .. + ! .. Structure Arguments .. + type (kdtree2), pointer :: tp + type(interval), intent(out) :: interv + ! .. + ! .. Scalar Arguments .. + integer(pInt), intent (In) :: c, l, u + ! .. + ! .. Local Scalars .. + real(pReal) :: last, lmax, lmin, t, smin,smax + integer(pInt) :: i, ulocal + ! .. + ! .. Local Arrays .. + real(pReal), pointer :: v(:,:) + integer(pInt), pointer :: ind(:) + ! .. + v => tp%the_data(1:,1:) + ind => tp%ind(1:) + smin = v(c,ind(l)) + smax = smin + + ulocal = u + + do i = l + 2, ulocal, 2 + lmin = v(c,ind(i-1)) + lmax = v(c,ind(i)) + if (lmin>lmax) then + t = lmin + lmin = lmax + lmax = t + end if + if (smin>lmin) smin = lmin + if (smaxlast) smin = last + if (smax qv + sr%nn = nn + sr%nfound = 0 + sr%centeridx = -1 + sr%correltime = 0 + sr%overflow = .false. + + sr%results => results + + sr%nalloc = nn ! will be checked + + sr%ind => tp%ind + sr%rearrange = tp%rearrange + if (tp%rearrange) then + sr%Data => tp%rearranged_data + else + sr%Data => tp%the_data + endif + sr%dimen = tp%dimen + + call validate_query_storage(nn) + sr%pq = pq_create(results) + + call search(tp%root) + + if (tp%sort) then + call kdtree2_sort_results(nn, results) + endif +! deallocate(sr%pqp) + return + end subroutine kdtree2_n_nearest + + subroutine kdtree2_n_nearest_around_point(tp,idxin,correltime,nn,results) + ! Find the 'nn' vectors in the tree nearest to point 'idxin', + ! with correlation window 'correltime', returing results in + ! results(:), which must be pre-allocated upon entry. + type (kdtree2), pointer :: tp + integer(pInt), intent (In) :: idxin, correltime, nn + type(kdtree2_result), target :: results(:) + + allocate (sr%qv(tp%dimen)) + sr%qv = tp%the_data(:,idxin) ! copy the vector + sr%ballsize = huge(1.0) ! the largest real(pReal) number + sr%centeridx = idxin + sr%correltime = correltime + + sr%nn = nn + sr%nfound = 0 + + sr%dimen = tp%dimen + sr%nalloc = nn + + sr%results => results + + sr%ind => tp%ind + sr%rearrange = tp%rearrange + + if (sr%rearrange) then + sr%Data => tp%rearranged_data + else + sr%Data => tp%the_data + endif + + call validate_query_storage(nn) + sr%pq = pq_create(results) + + call search(tp%root) + + if (tp%sort) then + call kdtree2_sort_results(nn, results) + endif + deallocate (sr%qv) + return + end subroutine kdtree2_n_nearest_around_point + + subroutine kdtree2_r_nearest(tp,qv,r2,nfound,nalloc,results) + ! find the nearest neighbors to point 'idxin', within SQUARED + ! Euclidean distance 'r2'. Upon ENTRY, nalloc must be the + ! size of memory allocated for results(1:nalloc). Upon + ! EXIT, nfound is the number actually found within the ball. + ! + ! Note that if nfound .gt. nalloc then more neighbors were found + ! than there were storage to store. The resulting list is NOT + ! the smallest ball inside norm r^2 + ! + ! Results are NOT sorted unless tree was created with sort option. + type (kdtree2), pointer :: tp + real(pReal), target, intent (In) :: qv(:) + real(pReal), intent(in) :: r2 + integer(pInt), intent(out) :: nfound + integer(pInt), intent (In) :: nalloc + type(kdtree2_result), target :: results(:) + + ! + sr%qv => qv + sr%ballsize = r2 + sr%nn = 0 ! flag for fixed ball search + sr%nfound = 0 + sr%centeridx = -1 + sr%correltime = 0 + + sr%results => results + + call validate_query_storage(nalloc) + sr%nalloc = nalloc + sr%overflow = .false. + sr%ind => tp%ind + sr%rearrange= tp%rearrange + + if (tp%rearrange) then + sr%Data => tp%rearranged_data + else + sr%Data => tp%the_data + endif + sr%dimen = tp%dimen + + ! + !sr%dsl = Huge(sr%dsl) ! set to huge positive values + !sr%il = -1 ! set to invalid indexes + ! + + call search(tp%root) + nfound = sr%nfound + if (tp%sort) then + call kdtree2_sort_results(nfound, results) + endif + + if (sr%overflow) then + write (*,*) 'KD_TREE_TRANS: warning! return from kdtree2_r_nearest found more neighbors' + write (*,*) 'KD_TREE_TRANS: than storage was provided for. Answer is NOT smallest ball' + write (*,*) 'KD_TREE_TRANS: with that number of neighbors! I.e. it is wrong.' + endif + + return + end subroutine kdtree2_r_nearest + + subroutine kdtree2_r_nearest_around_point(tp,idxin,correltime,r2,& + nfound,nalloc,results) + ! + ! Like kdtree2_r_nearest, but around a point 'idxin' already existing + ! in the data set. + ! + ! Results are NOT sorted unless tree was created with sort option. + ! + type (kdtree2), pointer :: tp + integer(pInt), intent (In) :: idxin, correltime, nalloc + real(pReal), intent(in) :: r2 + integer(pInt), intent(out) :: nfound + type(kdtree2_result), target :: results(:) + ! .. + ! .. Intrinsic Functions .. + intrinsic HUGE + ! .. + allocate (sr%qv(tp%dimen)) + sr%qv = tp%the_data(:,idxin) ! copy the vector + sr%ballsize = r2 + sr%nn = 0 ! flag for fixed r search + sr%nfound = 0 + sr%centeridx = idxin + sr%correltime = correltime + + sr%results => results + + sr%nalloc = nalloc + sr%overflow = .false. + + call validate_query_storage(nalloc) + + ! sr%dsl = HUGE(sr%dsl) ! set to huge positive values + ! sr%il = -1 ! set to invalid indexes + + sr%ind => tp%ind + sr%rearrange = tp%rearrange + + if (tp%rearrange) then + sr%Data => tp%rearranged_data + else + sr%Data => tp%the_data + endif + sr%rearrange = tp%rearrange + sr%dimen = tp%dimen + + ! + !sr%dsl = Huge(sr%dsl) ! set to huge positive values + !sr%il = -1 ! set to invalid indexes + ! + + call search(tp%root) + nfound = sr%nfound + if (tp%sort) then + call kdtree2_sort_results(nfound,results) + endif + + if (sr%overflow) then + write (*,*) 'KD_TREE_TRANS: warning! return from kdtree2_r_nearest found more neighbors' + write (*,*) 'KD_TREE_TRANS: than storage was provided for. Answer is NOT smallest ball' + write (*,*) 'KD_TREE_TRANS: with that number of neighbors! I.e. it is wrong.' + endif + + deallocate (sr%qv) + return + end subroutine kdtree2_r_nearest_around_point + + function kdtree2_r_count(tp,qv,r2) result(nfound) + ! Count the number of neighbors within square distance 'r2'. + type (kdtree2), pointer :: tp + real(pReal), target, intent (In) :: qv(:) + real(pReal), intent(in) :: r2 + integer(pInt) :: nfound + ! .. + ! .. Intrinsic Functions .. + intrinsic HUGE + ! .. + sr%qv => qv + sr%ballsize = r2 + + sr%nn = 0 ! flag for fixed r search + sr%nfound = 0 + sr%centeridx = -1 + sr%correltime = 0 + + nullify(sr%results) ! for some reason, FTN 95 chokes on '=> null()' + + sr%nalloc = 0 ! we do not allocate any storage but that's OK + ! for counting. + sr%ind => tp%ind + sr%rearrange = tp%rearrange + if (tp%rearrange) then + sr%Data => tp%rearranged_data + else + sr%Data => tp%the_data + endif + sr%dimen = tp%dimen + + ! + !sr%dsl = Huge(sr%dsl) ! set to huge positive values + !sr%il = -1 ! set to invalid indexes + ! + sr%overflow = .false. + + call search(tp%root) + + nfound = sr%nfound + + return + end function kdtree2_r_count + + function kdtree2_r_count_around_point(tp,idxin,correltime,r2) & + result(nfound) + ! Count the number of neighbors within square distance 'r2' around + ! point 'idxin' with decorrelation time 'correltime'. + ! + type (kdtree2), pointer :: tp + integer(pInt), intent (In) :: correltime, idxin + real(pReal), intent(in) :: r2 + integer(pInt) :: nfound + ! .. + ! .. + ! .. Intrinsic Functions .. + intrinsic HUGE + ! .. + allocate (sr%qv(tp%dimen)) + sr%qv = tp%the_data(:,idxin) + sr%ballsize = r2 + + sr%nn = 0 ! flag for fixed r search + sr%nfound = 0 + sr%centeridx = idxin + sr%correltime = correltime + nullify(sr%results) + + sr%nalloc = 0 ! we do not allocate any storage but that's OK + ! for counting. + + sr%ind => tp%ind + sr%rearrange = tp%rearrange + + if (sr%rearrange) then + sr%Data => tp%rearranged_data + else + sr%Data => tp%the_data + endif + sr%dimen = tp%dimen + + ! + !sr%dsl = Huge(sr%dsl) ! set to huge positive values + !sr%il = -1 ! set to invalid indexes + ! + sr%overflow = .false. + + call search(tp%root) + + nfound = sr%nfound + + return + end function kdtree2_r_count_around_point + + + subroutine validate_query_storage(n) + ! + ! make sure we have enough storage for n + ! + integer(pInt), intent(in) :: n + + if (size(sr%results,1) .lt. n) then + write (*,*) 'KD_TREE_TRANS: you did not provide enough storage for results(1:n)' + stop + return + endif + + return + end subroutine validate_query_storage + + function square_distance(d, iv,qv) result (res) + ! distance between iv[1:n] and qv[1:n] + ! .. Function Return Value .. + ! re-implemented to improve vectorization. + real(pReal) :: res + ! .. + ! .. + ! .. Scalar Arguments .. + integer(pInt) :: d + ! .. + ! .. Array Arguments .. + real(pReal) :: iv(:),qv(:) + ! .. + ! .. + res = sum( (iv(1:d)-qv(1:d))**2 ) + end function square_distance + + recursive subroutine search(node) + ! + ! This is the innermost core routine of the kd-tree search. Along + ! with "process_terminal_node", it is the performance bottleneck. + ! + ! This version uses a logically complete secondary search of + ! "box in bounds", whether the sear + ! + type (Tree_node), pointer :: node + ! .. + type(tree_node),pointer :: ncloser, nfarther + ! + integer(pInt) :: cut_dim, i + ! .. + real(pReal) :: qval, dis + real(pReal) :: ballsize + real(pReal), pointer :: qv(:) + type(interval), pointer :: box(:) + + if ((associated(node%left) .and. associated(node%right)) .eqv. .false.) then + ! we are on a terminal node + if (sr%nn .eq. 0) then + call process_terminal_node_fixedball(node) + else + call process_terminal_node(node) + endif + else + ! we are not on a terminal node + qv => sr%qv(1:) + cut_dim = node%cut_dim + qval = qv(cut_dim) + + if (qval < node%cut_val) then + ncloser => node%left + nfarther => node%right + dis = (node%cut_val_right - qval)**2 +! extra = node%cut_val - qval + else + ncloser => node%right + nfarther => node%left + dis = (node%cut_val_left - qval)**2 +! extra = qval- node%cut_val_left + endif + + if (associated(ncloser)) call search(ncloser) + + ! we may need to search the second node. + if (associated(nfarther)) then + ballsize = sr%ballsize +! dis=extra**2 + if (dis <= ballsize) then + ! + ! we do this separately as going on the first cut dimen is often + ! a good idea. + ! note that if extra**2 < sr%ballsize, then the next + ! check will also be false. + ! + box => node%box(1:) + do i=1,sr%dimen + if (i .ne. cut_dim) then + dis = dis + dis2_from_bnd(qv(i),box(i)%lower,box(i)%upper) + if (dis > ballsize) then + return + endif + endif + end do + + ! + ! if we are still here then we need to search mroe. + ! + call search(nfarther) + endif + endif + end if + end subroutine search + + + real(pReal) function dis2_from_bnd(x,amin,amax) result (res) + real(pReal), intent(in) :: x, amin,amax + + if (x > amax) then + res = (x-amax)**2; + return + else + if (x < amin) then + res = (amin-x)**2; + return + else + res = 0.0 + return + endif + endif + return + end function dis2_from_bnd + + logical function box_in_search_range(node, sr) result(res) + ! + ! Return the distance from 'qv' to the CLOSEST corner of node's + ! bounding box + ! for all coordinates outside the box. Coordinates inside the box + ! contribute nothing to the distance. + ! + type (tree_node), pointer :: node + type (tree_search_record), pointer :: sr + + integer(pInt) :: dimen, i + real(pReal) :: dis, ballsize + real(pReal) :: l, u + + dimen = sr%dimen + ballsize = sr%ballsize + dis = 0.0 + res = .true. + do i=1,dimen + l = node%box(i)%lower + u = node%box(i)%upper + dis = dis + (dis2_from_bnd(sr%qv(i),l,u)) + if (dis > ballsize) then + res = .false. + return + endif + end do + res = .true. + return + end function box_in_search_range + + + subroutine process_terminal_node(node) + ! + ! Look for actual near neighbors in 'node', and update + ! the search results on the sr data structure. + ! + type (tree_node), pointer :: node + ! + real(pReal), pointer :: qv(:) + integer(pInt), pointer :: ind(:) + real(pReal), pointer :: data(:,:) + ! + integer(pInt) :: dimen, i, indexofi, k, centeridx, correltime + real(pReal) :: ballsize, sd, newpri + logical :: rearrange + type(pq), pointer :: pqp + ! + ! copy values from sr to local variables + ! + ! + ! Notice, making local pointers with an EXPLICIT lower bound + ! seems to generate faster code. + ! why? I don't know. + qv => sr%qv(1:) + pqp => sr%pq + dimen = sr%dimen + ballsize = sr%ballsize + rearrange = sr%rearrange + ind => sr%ind(1:) + data => sr%Data(1:,1:) + centeridx = sr%centeridx + correltime = sr%correltime + + ! doing_correl = (centeridx >= 0) ! Do we have a decorrelation window? + ! include_point = .true. ! by default include all points + ! search through terminal bucket. + + mainloop: do i = node%l, node%u + if (rearrange) then + sd = 0.0 + do k = 1,dimen + sd = sd + (data(k,i) - qv(k))**2 + if (sd>ballsize) cycle mainloop + end do + indexofi = ind(i) ! only read it if we have not broken out + else + indexofi = ind(i) + sd = 0.0 + do k = 1,dimen + sd = sd + (data(k,indexofi) - qv(k))**2 + if (sd>ballsize) cycle mainloop + end do + endif + + if (centeridx > 0) then ! doing correlation interval? + if (abs(indexofi-centeridx) < correltime) cycle mainloop + endif + + + ! + ! two choices for any point. The list so far is either undersized, + ! or it is not. + ! + ! If it is undersized, then add the point and its distance + ! unconditionally. If the point added fills up the working + ! list then set the sr%ballsize, maximum distance bound (largest distance on + ! list) to be that distance, instead of the initialized +infinity. + ! + ! If the running list is full size, then compute the + ! distance but break out immediately if it is larger + ! than sr%ballsize, "best squared distance" (of the largest element), + ! as it cannot be a good neighbor. + ! + ! Once computed, compare to best_square distance. + ! if it is smaller, then delete the previous largest + ! element and add the new one. + + if (sr%nfound .lt. sr%nn) then + ! + ! add this point unconditionally to fill list. + ! + sr%nfound = sr%nfound +1 + newpri = pq_insert(pqp,sd,indexofi) + if (sr%nfound .eq. sr%nn) ballsize = newpri + ! we have just filled the working list. + ! put the best square distance to the maximum value + ! on the list, which is extractable from the PQ. + + else + ! + ! now, if we get here, + ! we know that the current node has a squared + ! distance smaller than the largest one on the list, and + ! belongs on the list. + ! Hence we replace that with the current one. + ! + ballsize = pq_replace_max(pqp,sd,indexofi) + endif + end do mainloop + ! + ! Reset sr variables which may have changed during loop + ! + sr%ballsize = ballsize + + end subroutine process_terminal_node + + subroutine process_terminal_node_fixedball(node) + ! + ! Look for actual near neighbors in 'node', and update + ! the search results on the sr data structure, i.e. + ! save all within a fixed ball. + ! + type (tree_node), pointer :: node + ! + real(pReal), pointer :: qv(:) + integer(pInt), pointer :: ind(:) + real(pReal), pointer :: data(:,:) + ! + integer(pInt) :: nfound + integer(pInt) :: dimen, i, indexofi, k + integer(pInt) :: centeridx, correltime, nn + real(pReal) :: ballsize, sd + logical :: rearrange + + ! + ! copy values from sr to local variables + ! + qv => sr%qv(1:) + dimen = sr%dimen + ballsize = sr%ballsize + rearrange = sr%rearrange + ind => sr%ind(1:) + data => sr%Data(1:,1:) + centeridx = sr%centeridx + correltime = sr%correltime + nn = sr%nn ! number to search for + nfound = sr%nfound + + ! search through terminal bucket. + mainloop: do i = node%l, node%u + + ! + ! two choices for any point. The list so far is either undersized, + ! or it is not. + ! + ! If it is undersized, then add the point and its distance + ! unconditionally. If the point added fills up the working + ! list then set the sr%ballsize, maximum distance bound (largest distance on + ! list) to be that distance, instead of the initialized +infinity. + ! + ! If the running list is full size, then compute the + ! distance but break out immediately if it is larger + ! than sr%ballsize, "best squared distance" (of the largest element), + ! as it cannot be a good neighbor. + ! + ! Once computed, compare to best_square distance. + ! if it is smaller, then delete the previous largest + ! element and add the new one. + + ! which index to the point do we use? + + if (rearrange) then + sd = 0.0 + do k = 1,dimen + sd = sd + (data(k,i) - qv(k))**2 + if (sd>ballsize) cycle mainloop + end do + indexofi = ind(i) ! only read it if we have not broken out + else + indexofi = ind(i) + sd = 0.0 + do k = 1,dimen + sd = sd + (data(k,indexofi) - qv(k))**2 + if (sd>ballsize) cycle mainloop + end do + endif + + if (centeridx > 0) then ! doing correlation interval? + if (abs(indexofi-centeridx) 1)then + ileft=ileft-1 + value=a(ileft); ivalue=ind(ileft) + else + value=a(iright); ivalue=ind(iright) + a(iright)=a(1); ind(iright)=ind(1) + iright=iright-1 + if (iright == 1) then + a(1)=value;ind(1)=ivalue + return + endif + endif + i=ileft + j=2*ileft + do while (j <= iright) + if(j < iright) then + if(a(j) < a(j+1)) j=j+1 + endif + if(value < a(j)) then + a(i)=a(j); ind(i)=ind(j) + i=j + j=j+j + else + j=iright+1 + endif + end do + a(i)=value; ind(i)=ivalue + end do + end subroutine heapsort + + subroutine heapsort_struct(a,n) + ! + ! Sort a(1:n) in ascending order + ! + ! + integer(pInt),intent(in) :: n + type(kdtree2_result),intent(inout) :: a(:) + + ! + ! + type(kdtree2_result) :: value ! temporary value + + integer(pInt) :: i,j + integer(pInt) :: ileft,iright + + ileft=n/2+1 + iright=n + + ! do i=1,n + ! ind(i)=i + ! Generate initial idum array + ! end do + + if(n.eq.1) return + + do + if(ileft > 1)then + ileft=ileft-1 + value=a(ileft) + else + value=a(iright) + a(iright)=a(1) + iright=iright-1 + if (iright == 1) then + a(1) = value + return + endif + endif + i=ileft + j=2*ileft + do while (j <= iright) + if(j < iright) then + if(a(j)%dis < a(j+1)%dis) j=j+1 + endif + if(value%dis < a(j)%dis) then + a(i)=a(j); + i=j + j=j+j + else + j=iright+1 + endif + end do + a(i)=value + end do + end subroutine heapsort_struct + +end module kdtree2_module +!############################################################################################################################# +! END KDTREE2 +!############################################################################################################################# diff --git a/code/numerics.f90 b/code/numerics.f90 index fadec4049..6e87dca61 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -71,7 +71,7 @@ real(pReal) :: relevantStrain, & ! strain fftw_timelimit, & ! sets the timelimit of plan creation for FFTW, see manual on www.fftw.org rotation_tol ! tolerance of rotation specified in loadcase character(len=64) :: fftw_planner_string ! reads the planing-rigor flag, see manual on www.fftw.org -integer*8 :: fftw_planner_flag ! conversion of fftw_planner_string to integer, basically what is usually done in the include file of fftw +integer(pInt) :: fftw_planner_flag ! conversion of fftw_planner_string to integer, basically what is usually done in the include file of fftw logical :: memory_efficient,& ! for fast execution (pre calculation of gamma_hat) divergence_correction ! correct divergence calculation in fourier space integer(pInt) :: itmax , & ! maximum number of iterations