additional changes to other programs

EMsoft-org · Jul 14, 2020 · 6606859 · 6606859
1 parent 30b3015
commit 6606859
Show file tree

Hide file tree

Showing 4 changed files with 108 additions and 25 deletions.
diff --git a/Source/EMOpenCLLib/Indexingmod.f90 b/Source/EMOpenCLLib/Indexingmod.f90
@@ -241,7 +241,6 @@ END SUBROUTINE ProgCallBackTypeErrorDIdriver
 real(kind=sgl),allocatable                          :: EBSDdictpatflt(:,:)
 real(kind=dbl),allocatable                          :: rdata(:,:), fdata(:,:), rrdata(:,:), ffdata(:,:), ksqarray(:,:)
 complex(kind=dbl),allocatable                       :: hpmask(:,:)
-complex(C_DOUBLE_COMPLEX),allocatable               :: inp(:,:), outp(:,:)
 real(kind=dbl)                                      :: w, Jres
 integer(kind=irg)                                   :: dims(2)
 character(11)                                       :: dstr

diff --git a/Source/EMsoftHDFLib/patternmod.f90 b/Source/EMsoftHDFLib/patternmod.f90
@@ -1921,7 +1921,7 @@ recursive subroutine PreProcessTKDPatterns(nthreads, inRAM, tkdnl, binx, biny, m
 use omp_lib
 use filters
 use timing
-
+use FFTW3mod
 use ISO_C_BINDING
 
 IMPLICIT NONE
@@ -1952,7 +1952,8 @@ recursive subroutine PreProcessTKDPatterns(nthreads, inRAM, tkdnl, binx, biny, m
 real(kind=sgl),allocatable                  :: tmpimageexpt(:), TKDPattern(:,:), imageexpt(:), exppatarray(:), TKDpat(:,:)
 real(kind=dbl),allocatable                  :: rrdata(:,:), ffdata(:,:), ksqarray(:,:)
 complex(kind=dbl),allocatable               :: hpmask(:,:)
-complex(C_DOUBLE_COMPLEX),allocatable       :: inp(:,:), outp(:,:)
+complex(C_DOUBLE_COMPLEX),pointer           :: inp(:,:), outp(:,:)
+type(c_ptr), allocatable                    :: ip, op
 type(C_PTR)                                 :: planf, HPplanf, HPplanb
 
 
@@ -2037,10 +2038,26 @@ recursive subroutine PreProcessTKDPatterns(nthreads, inRAM, tkdnl, binx, biny, m
 end if
 
 ! initialize the HiPassFilter routine (has its own FFTW plans)
-allocate(hpmask(binx,biny),inp(binx,biny),outp(binx,biny),stat=istat)
+allocate(hpmask(binx,biny),stat=istat)
 if (istat .ne. 0) stop 'could not allocate hpmask array'
+
+! use the fftw_alloc routine to create the inp and outp arrays
+! using a regular allocate can occasionally cause issues, in particular with 
+! the ifort compiler. [MDG, 7/14/20]
+ip = fftw_alloc_complex(int(binx*biny,C_SIZE_T))
+call c_f_pointer(ip, inp, [binx,biny])
+
+op = fftw_alloc_complex(int(binx*biny,C_SIZE_T))
+call c_f_pointer(op, outp, [binx,biny])
+
+inp = cmplx(0.D0,0D0)
+outp = cmplx(0.D0,0.D0)
+
 call init_HiPassFilter(w, (/ binx, biny /), hpmask, inp, outp, HPplanf, HPplanb) 
-deallocate(inp, outp)
+
+! remove these pointers again because we will need them for the parallel part
+call fftw_free(ip)
+call fftw_free(op)
 
 call Message('Starting processing of experimental patterns')
 call Time_tick(tickstart)
@@ -2052,7 +2069,7 @@ recursive subroutine PreProcessTKDPatterns(nthreads, inRAM, tkdnl, binx, biny, m
 prepexperimentalloop: do iii = iiistart,iiiend
 
 ! start the OpenMP portion
-!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(TID, jj, kk, mi, ma, istat) &
+!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(TID, jj, kk, mi, ma, istat, ip, op) &
 !$OMP& PRIVATE(imageexpt, tmpimageexpt, TKDPat, rrdata, ffdata, TKDpint, vlen, tmp, inp, outp)
 
 ! set the thread ID
@@ -2065,8 +2082,14 @@ recursive subroutine PreProcessTKDPatterns(nthreads, inRAM, tkdnl, binx, biny, m
     allocate(TKDpint(binx,biny),stat=istat)
     if (istat .ne. 0) stop 'could not allocate TKDpint array'
 
-    allocate(inp(binx,biny),outp(binx,biny),stat=istat)
-    if (istat .ne. 0) stop 'could not allocate inp, outp arrays'
+    ip = fftw_alloc_complex(int(binx*biny,C_SIZE_T))
+    call c_f_pointer(ip, inp, [binx,biny])
+
+    op = fftw_alloc_complex(int(binx*biny,C_SIZE_T))
+    call c_f_pointer(op, outp, [binx,biny])
+
+    inp = cmplx(0.D0,0D0)
+    outp = cmplx(0.D0,0.D0)
 
     tmpimageexpt = 0.0
     rrdata = 0.D0
@@ -2143,7 +2166,10 @@ recursive subroutine PreProcessTKDPatterns(nthreads, inRAM, tkdnl, binx, biny, m
       end if
     end if
 
-deallocate(tmpimageexpt, TKDPat, rrdata, ffdata, TKDpint, inp, outp)
+deallocate(tmpimageexpt, TKDPat, rrdata, ffdata, TKDpint)
+call fftw_free(ip)
+call fftw_free(op)
+call fftw_cleanup()
 !$OMP BARRIER
 !$OMP END PARALLEL
 

diff --git a/Source/EMsoftLib/filters.f90 b/Source/EMsoftLib/filters.f90
@@ -995,7 +995,7 @@ recursive subroutine init_HiPassFilter(w, dims, hpmask, inp, outp, planf, planb)
 integer(kind=irg),INTENT(IN)            :: dims(2)
 complex(kind=dbl),INTENT(INOUT)         :: hpmask(dims(1),dims(2))
 !f2py intent(in,out) ::  hpmask
-complex(C_DOUBLE_COMPLEX),INTENT(INOUT) :: inp(dims(1),dims(2)), outp(dims(1),dims(2))
+complex(C_DOUBLE_COMPLEX),pointer,INTENT(INOUT) :: inp(:,:), outp(:,:)
 !f2py intent(in,out) ::  inp
 type(C_PTR),INTENT(INOUT)               :: planf, planb
 !f2py intent(in,out) ::  planf, planb
@@ -1058,7 +1058,7 @@ recursive function applyHiPassFilter(rdata, dims, w, hpmask, inp, outp, planf, p
 real(kind=dbl),INTENT(IN)               :: w
 real(kind=dbl),INTENT(IN)               :: rdata(dims(1),dims(2))
 complex(kind=dbl),INTENT(IN)            :: hpmask(dims(1),dims(2))
-complex(C_DOUBLE_COMPLEX),INTENT(INOUT) :: inp(dims(1),dims(2)), outp(dims(1),dims(2))
+complex(C_DOUBLE_COMPLEX),pointer,INTENT(INOUT) :: inp(:,:), outp(:,:)
 !f2py intent(in,out) ::  inp
 type(C_PTR),INTENT(IN)                  :: planf, planb
 real(kind=dbl)                          :: fdata(dims(1),dims(2))

diff --git a/Source/EMsoftWrapperLib/DictionaryIndexing/EMDIwrappermod.f90 b/Source/EMsoftWrapperLib/DictionaryIndexing/EMDIwrappermod.f90
@@ -233,6 +233,8 @@ recursive subroutine EMsoftCpreprocessEBSDPatterns(ipar, fpar, spar, mask, exptI
 use EBSDDImod
 use omp_lib
 use patternmod
+use FFTW3mod
+
 
 IMPLICIT NONE
 
@@ -261,7 +263,8 @@ recursive subroutine EMsoftCpreprocessEBSDPatterns(ipar, fpar, spar, mask, exptI
 real(kind=sgl),allocatable              :: EBSDpatternintd(:,:), EBSDpat(:,:), exppatarray(:)
 real(kind=dbl),allocatable              :: ksqarray(:,:), rrdata(:,:), ffdata(:,:)
 complex(kind=dbl),allocatable           :: hpmask(:,:)
-complex(C_DOUBLE_COMPLEX),allocatable   :: inp(:,:), outp(:,:)
+complex(C_DOUBLE_COMPLEX),pointer       :: inp(:,:), outp(:,:)
+type(c_ptr),allocatable                 :: ip, op
 real(kind=sgl),allocatable              :: lstore(:,:), pstore(:,:), lp(:), cp(:)
 
 
@@ -453,10 +456,27 @@ recursive subroutine EMsoftCpreprocessEBSDPatterns(ipar, fpar, spar, mask, exptI
 deallocate(EBSDPat)
 
 ! initialize the HiPassFilter routine (has its own FFTW plans)
-allocate(hpmask(binx,biny),inp(binx,biny),outp(binx,biny),stat=istat)
+allocate(hpmask(binx,biny),stat=istat)
 ! if (istat .ne. 0) stop 'could not allocate hpmask array'
+
+! use the fftw_alloc routine to create the inp and outp arrays
+! using a regular allocate can occasionally cause issues, in particular with 
+! the ifort compiler. [MDG, 7/14/20]
+ip = fftw_alloc_complex(int(binx*biny,C_SIZE_T))
+call c_f_pointer(ip, inp, [binx,biny])
+
+op = fftw_alloc_complex(int(binx*biny,C_SIZE_T))
+call c_f_pointer(op, outp, [binx,biny])
+
+inp = cmplx(0.D0,0D0)
+outp = cmplx(0.D0,0.D0)
+
+
 call init_HiPassFilter(w, (/ binx, biny /), hpmask, inp, outp, HPplanf, HPplanb) 
-deallocate(inp, outp)
+
+! remove these pointers again because we will need them for the parallel part
+call fftw_free(ip)
+call fftw_free(op)
 
 dims3 = (/ binx, biny, ipar(26) /)
 
@@ -488,7 +508,7 @@ recursive subroutine EMsoftCpreprocessEBSDPatterns(ipar, fpar, spar, mask, exptI
 prepexperimentalloop: do iii = iiistart,iiiend
 
 ! start the OpenMP portion
-!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(TID, jj, kk, mi, ma, istat) &
+!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(TID, jj, kk, mi, ma, istat, ip, op) &
 !$OMP& PRIVATE(imageexpt, tmpimageexpt, EBSDPat, rrdata, ffdata, EBSDpint, vlen, tmp, inp, outp)
 
 ! set the thread ID
@@ -501,8 +521,14 @@ recursive subroutine EMsoftCpreprocessEBSDPatterns(ipar, fpar, spar, mask, exptI
     allocate(EBSDpint(binx,biny),stat=istat)
     ! if (istat .ne. 0) stop 'could not allocate EBSDpint array'
 
-    allocate(inp(binx,biny),outp(binx,biny),stat=istat)
-    ! if (istat .ne. 0) stop 'could not allocate inp, outp arrays'
+    ip = fftw_alloc_complex(int(binx*biny,C_SIZE_T))
+    call c_f_pointer(ip, inp, [binx,biny])
+
+    op = fftw_alloc_complex(int(binx*biny,C_SIZE_T))
+    call c_f_pointer(op, outp, [binx,biny])
+
+    inp = cmplx(0.D0,0D0)
+    outp = cmplx(0.D0,0.D0)
 
     tmpimageexpt = 0.0
     rrdata = 0.D0
@@ -605,7 +631,10 @@ recursive subroutine EMsoftCpreprocessEBSDPatterns(ipar, fpar, spar, mask, exptI
       end if
     end if
 
-deallocate(tmpimageexpt, EBSDPat, rrdata, ffdata, EBSDpint, inp, outp)
+deallocate(tmpimageexpt, EBSDPat, rrdata, ffdata, EBSDpint)
+call fftw_free(ip)
+call fftw_free(op)
+call fftw_cleanup()
 !$OMP BARRIER
 !$OMP END PARALLEL
 
@@ -674,6 +703,7 @@ recursive subroutine EMsoftCpreprocessSingleEBSDPattern(ipar, fpar, inputpattern
 use EBSDDImod
 use omp_lib
 use patternmod
+use FFTW3mod
 
 IMPLICIT NONE
 
@@ -688,8 +718,8 @@ recursive subroutine EMsoftCpreprocessSingleEBSDPattern(ipar, fpar, inputpattern
 real(kind=sgl),allocatable              :: EBSDpat(:,:)
 real(kind=dbl),allocatable              :: rrdata(:,:), ffdata(:,:)
 complex(kind=dbl),allocatable           :: hpmask(:,:)
-complex(C_DOUBLE_COMPLEX),allocatable   :: inp(:,:), outp(:,:)
-
+complex(C_DOUBLE_COMPLEX),pointer       :: inp(:,:), outp(:,:)
+type(c_ptr), allocatable                :: ip, op
 type(C_PTR)                             :: HPplanf, HPplanb
 
 HPplanf = C_NULL_PTR
@@ -715,7 +745,17 @@ recursive subroutine EMsoftCpreprocessSingleEBSDPattern(ipar, fpar, inputpattern
 biny = ipar(20)
 
 ! initialize the HiPassFilter routine (has its own FFTW plans)
-allocate(hpmask(binx,biny),inp(binx,biny),outp(binx,biny),stat=istat)
+allocate(hpmask(binx,biny),stat=istat)
+
+ip = fftw_alloc_complex(int(binx*biny,C_SIZE_T))
+call c_f_pointer(ip, inp, [binx,biny])
+
+op = fftw_alloc_complex(int(binx*biny,C_SIZE_T))
+call c_f_pointer(op, outp, [binx,biny])
+
+inp = cmplx(0.D0,0D0)
+outp = cmplx(0.D0,0.D0)
+
 call init_HiPassFilter(dble(fpar(24)), (/ binx, biny /), hpmask, inp, outp, HPplanf, HPplanb) 
 
 ! initialize thread private variables
@@ -734,10 +774,14 @@ recursive subroutine EMsoftCpreprocessSingleEBSDPattern(ipar, fpar, inputpattern
 EBSDpint = nint(((EBSDPat - mi) / (ma-mi))*255.0)
 outputpattern = float(adhisteq(ipar(28), binx, biny, EBSDpint))
 
-deallocate(inp, outp, EBSDpint, EBSDPat, rrdata, ffdata, hpmask)
+deallocate(EBSDpint, EBSDPat, rrdata, ffdata, hpmask)
 HPplanf = C_NULL_PTR
 HPplanb = C_NULL_PTR
 
+call fftw_free(ip)
+call fftw_free(op)
+call fftw_cleanup()
+
 ! that's it folks...
 end subroutine EMsoftCpreprocessSingleEBSDPattern
 
@@ -777,6 +821,7 @@ recursive subroutine EMsoftCEBSDDIpreview(ipar, fpar, spar, averagedpattern, pat
 use omp_lib
 use patternmod
 use HDF5
+use FFTW3mod
 
 IMPLICIT NONE
 
@@ -799,7 +844,8 @@ recursive subroutine EMsoftCEBSDDIpreview(ipar, fpar, spar, averagedpattern, pat
 integer(HSIZE_T)                        :: dims3(3), offset3(3)
 type(C_PTR)                             :: HPplanf, HPplanb
 complex(kind=dbl),allocatable           :: hpmask(:,:)
-complex(C_DOUBLE_COMPLEX),allocatable   :: inp(:,:), outp(:,:)
+complex(C_DOUBLE_COMPLEX),pointer       :: inp(:,:), outp(:,:)
+type(c_ptr),allocatable                 :: ip, op 
 real(kind=dbl),allocatable              :: rrdata(:,:), ffdata(:,:), ksqarray(:,:)
 
 ! parameters to deal with the input string array spar
@@ -969,8 +1015,17 @@ recursive subroutine EMsoftCEBSDDIpreview(ipar, fpar, spar, averagedpattern, pat
 end do
 
 ! next we need to set up the high-pass filter fftw plans
-allocate(hpmask(binx,biny),inp(binx,biny),outp(binx,biny))
+allocate(hpmask(binx,biny))
 allocate(rrdata(binx,biny),ffdata(binx,biny))
+ip = fftw_alloc_complex(int(binx*biny,C_SIZE_T))
+call c_f_pointer(ip, inp, [binx,biny])
+
+op = fftw_alloc_complex(int(binx*biny,C_SIZE_T))
+call c_f_pointer(op, outp, [binx,biny])
+
+inp = cmplx(0.D0,0D0)
+outp = cmplx(0.D0,0.D0)
+
 call init_HiPassFilter(dble(hpvals(1)), (/numsx, numsy /), hpmask, inp, outp, HPplanf, HPplanb) 
 
 ! the outer loop goes over the hipass filter width and is displayed horizontally in the final image
@@ -1023,9 +1078,12 @@ recursive subroutine EMsoftCEBSDDIpreview(ipar, fpar, spar, averagedpattern, pat
 end do
 
 ! clean up all auxiliary arrays
-deallocate(hpmask, inp, outp, pcopy, rrdata, ffdata, hpvals, nrvals)
+deallocate(hpmask, pcopy, rrdata, ffdata, hpvals, nrvals)
 deallocate(pattern, pint, ppp, expt)
 if (allocated(sumexpt)) deallocate(sumexpt)
+call fftw_free(ip)
+call fftw_free(op)
+call fftw_cleanup()
 
 ! this completes the computation of the patternarray output array