fixing merge conflict

.gitignore

@@ -23,3 +23,4 @@ lib/*
 *.sh
 sbatch*
 *.obj
+make.inc

docs/fortran-c.rst

@@ -355,7 +355,7 @@ and the associated pressure
 Here $x_{j}$ are the source locations,
 $\sigma_{j}$ are the Stokeslet densities,
 $\nu_{j}$ are the stresslet orientation vectors, $\mu_{j}$ 
-are the stresslet densities, and rhw xollwxrion of $x$ 
+are the stresslet densities, and the locations $x$ 
 at which the velocity and its gradient are evaluated are referred to
 as the evaluation points.
 

make.inc.icc

@@ -1,7 +1,7 @@
 CC=icc
 CXX=icpc
 FC=ifort
-FFLAGS= -fPIC -O3 -march=native -funroll-loops -mkl 
+FFLAGS= -fPIC -O3 -march=native -funroll-loops -mkl -w 
 LIBS=-lm 
 CLIBS = -lm -ldl -lifcore
 

make.inc.macos.gnu

@@ -1,6 +1,6 @@
 # makefile overrides
 # OS:       macOS
-# Compiler: gfortran X.X
+# Compiler: gfortran X.X/Clang
 # OpenMP:   enabled
 #
 
@@ -22,8 +22,6 @@ OMPLIBS = -lgomp
 # MATLAB interface:
 FDIR=$$(dirname `gfortran --print-file-name libgfortran.dylib`)
 MFLAGS +=-L${FDIR}
-MEX = $(shell ls -d /Applications/MATLAB_R20**.app)/bin/mex
-#LIBS = -lm -lstdc++.6
-#MEXLIBS= -lm -lstdc++.6 -lgfortran -ldl
+MEX = $(shell ls -d /Applications/MATLAB_R* | sort | tail -1)/bin/mex
 
 

make.inc.macos.intel

@@ -7,7 +7,8 @@
 CC=icc
 CXX=icpc
 FC=ifort
-FFLAGS= -fPIC -O3 -march=native -funroll-loops -qmkl
+
+FFLAGS= -fPIC -O3 -march=native -funroll-loops -mkl -std=legacy -w 
 LIBS=
 #CLIBS = -lm -ldl -lifcore
 CLIBS = -lm -ldl
@@ -28,4 +29,9 @@ endif
 OMPFLAGS = -qopenmp
 OMPLIBS = -qopenmp
 
+# MATLAB interface:
+FDIR=$$(dirname `gfortran --print-file-name libgfortran.dylib`)
+MFLAGS +=-L${FDIR}
+MEX = $(shell ls -d /Applications/MATLAB_R* | sort | tail -1)/bin/mex
+
 

make.inc.macos.nag

@@ -0,0 +1,30 @@
+# make.inc for NAG Fortran compiler
+# Online documentation: https://www.nag.com/nagware/np/r71_doc/manual/compiler_2_4.html#OPTIONS
+
+FC=nagfor
+
+# The path of libraries by NAG compiler
+LIB_NAG = /usr/local/lib/NAG_Fortran
+
+# Brief descriptions of specified options below:
+# -PIC: produce position-independent code
+# -O2: optimization at a normal level
+# -Ounroll=2: the depth of loo-unrolling
+# -f90_sign: use the Fortran 77/90 version of the SIGN intrinsic instead of the Fortran 95 one 
+# -dcfuns: enable recognition of non-standard double precision complex intrinsic functions.  
+# -dusty: allows the compilation and execution of legacy software.
+# -w=x77: suppresses extension warnings for obsolete but common extensions to Fortran 77.
+# -w=unreffed: suppresses warning messages about variables set but never referenced.
+# -w=unused: suppresses warning messages about unused entities. 
+# -ieee=full: set the mode of IEEE arithmetic operation according to full mode.
+
+# Main compile command for NAG Fortran compiler
+FFLAGS = -PIC -O2 -Ounroll=1 -f90_sign -dcfuns -dusty -w=obs -w=x77 -w=unreffed -w=unused -ieee=full 
+
+# Flags overwritten in makefile
+OMPFLAGS = -openmp
+# OMPLIBS = -lf71omp64 -L$(LIB_NAG)
+OMPLIBS = -lf71omp64 -lf71rts -L$(LIB_NAG)
+LIBS = -lf71rts -L$(LIB_NAG)
+CLIBS = -lm -ldl -L$(LIB_NAG)
+FFLAGS_DYN = -PIC

make.inc.windows.mingw

@@ -4,7 +4,7 @@
 # OpenMP:   default enabled unless specified
 #
 
-FFLAGS= -fPIC -O3 -funroll-loops -std=legacy 
+FFLAGS= -fPIC -O3 -funroll-loops -std=legacy -w 
 
 DYNAMICLIB = $(LIBNAME).dll
 LIMPLIB = $(LIBNAME)_dll.lib

makefile

@@ -14,7 +14,8 @@ FC=gfortran
 
 
 # set compiler flags for c and fortran
-FFLAGS= -fPIC -O3 -march=native -funroll-loops -std=legacy
+FFLAGS= -fPIC -O3 -march=native -funroll-loops -std=legacy -w
+FFLAGS_DYN= -shared -fPIC
 CFLAGS= -fPIC -O3 -march=native -funroll-loops -std=c99
 CXXFLAGS= -std=c++11 -DSCTL_PROFILE=-1 -fPIC -O3 -march=native -funroll-loops
 
@@ -86,6 +87,7 @@ endif
 # vectorized kernel directory
 SRCDIR = ./vec-kernels/src
 INCDIR = ./vec-kernels/include
+FINCDIR = ./src/Helmholtz
 LIBDIR = lib-static
 
 # objects to compile
@@ -192,10 +194,10 @@ usage:
 	$(CXX) -c $(CXXFLAGS) $< -o $@
 %.o: %.c %.h
 	$(CC) -c $(CFLAGS) $< -o $@
-%.o: %.f %.h
-	$(FC) -c $(FFLAGS) $< -o $@
+%.o: %.f 
+	$(FC) -c $(FFLAGS) -I$(FINCDIR) $< -o $@
 %.o: %.f90
-	$(FC) -c $(FFLAGS) $< -o $@
+	$(FC) -c $(FFLAGS) -I$(FINCDIR) $< -o $@
 
 # build the library...
 lib: $(STATICLIB) $(DYNAMICLIB)
@@ -224,7 +226,7 @@ $(STATICLIB): $(OBJS)
 	ar rcs $(STATICLIB) $(OBJS)
 	mv $(STATICLIB) lib-static/
 $(DYNAMICLIB): $(OBJS)
-	$(FC) -shared -fPIC $(OBJS) -o $(DYNAMICLIB) $(DYLIBS)
+	$(FC) $(FFLAGS_DYN) $(OBJS) -o $(DYNAMICLIB) $(DYLIBS)
 	mv $(DYNAMICLIB) lib/
 	[ ! -f $(LIMPLIB) ] || mv $(LIMPLIB) lib/
 

src/Common/fmmcommon.f

@@ -218,8 +218,8 @@ subroutine ireorderi(ndim,n,arr,arrsort,iarr)
 c
       subroutine drescale(n,a,r)
       implicit none
-      real *8 a(n),r
       integer i,n
+      real *8 a(n),r
 
 C$OMP PARALLEL DO DEFAULT(SHARED)
       do i=1,n

src/Common/tree_routs3d.f

@@ -1125,11 +1125,11 @@ subroutine getlist4pwdirtest(dir,censrc,centrg,boxsize)
       subroutine subdividebox(pos,npts,center,boxsize,
      1           isorted,iboxfl,subcenters)
       implicit none
+      integer npts
       double precision pos(3,npts)
       double precision center(3)
       double precision subcenters(3,8)
       double precision boxsize
-      integer npts
       integer isorted(*)
       integer iboxfl(2,8)
 

src/Common/yrecursion.f

@@ -1430,7 +1430,8 @@ subroutine zylgndrbr(nmax, z, y)
 c     branch cut at (0,+i), select the lower branch 
 c     of complex square root
 c
-      if( imag(1-z*z) .gt. 0 .and. real(1-z*z) .lt. 0) u=+sqrt(1-z*z)
+c      if( imag(1-z*z) .gt. 0 .and. real(1-z*z) .lt. 0) u=+sqrt(1-z*z)
+      if( dimag(1-z*z) .gt. 0 .and. real(1-z*z) .lt. 0) u=+sqrt(1-z*z)
 ccc      call prin2('in zylgndrbr, u=*', -u, 2)
 ccc      call prin2('in zylgndrbr, 1-z^2=*', 1-z*z, 2)
 c
@@ -1498,8 +1499,10 @@ subroutine zylgndrsc(nmax, z,scale, ysc)
 c
       ztmp = 1-z*z
       u=-sqrt(ztmp)
-      if(abs(imag(z)).le.1.0d-16.and.abs(real(z)).gt.1) then
-        if(imag(u).lt.0) u = dconjg(u)
+c      if(abs(imag(z)).le.1.0d-16.and.abs(real(z)).gt.1) then
+c        if(imag(u).lt.0) u = dconjg(u)
+      if(abs(dimag(z)).le.1.0d-16.and.abs(real(z)).gt.1) then
+        if(dimag(u).lt.0) u = dconjg(u)
       endif
       ysc(0,0)=1
       do m=0, nmax

src/Helmholtz/hfmm3d.f

@@ -562,6 +562,7 @@ subroutine hfmm3dmain(nd,eps,zk,
       double complex jsort(nd,0:ntj,-ntj:ntj,nexpc)
 
 
+      integer nboxes
       integer *8 iaddr(2,nboxes), lmptot
       double precision rmlexp(lmptot)
 
@@ -575,7 +576,6 @@ subroutine hfmm3dmain(nd,eps,zk,
       integer nterms(0:nlevels)
       integer *8 ipointer(8),ltree
       integer itree(ltree)
-      integer nboxes
       double precision rscales(0:nlevels)
       double precision boxsize(0:nlevels)
       integer isrcse(2,nboxes),itargse(2,nboxes),iexpcse(2,nboxes)
@@ -752,7 +752,7 @@ subroutine hfmm3dmain(nd,eps,zk,
 
       zkiupbound = 12*pi
       zkrupbound = 16*pi
-      zi = imag(zk)
+      zi = dimag(zk)
 
       ilevcutoff = -1
 
@@ -1054,7 +1054,7 @@ subroutine hfmm3dmain(nd,eps,zk,
         allocate(iboxlexp(nd*(nterms(ilev)+1)*
      1           (2*nterms(ilev)+1),8,nthd))
         zk2 = zk*boxsize(ilev)
-        if(real(zk2).le.zkrupbound.and.imag(zk2).lt.zkiupbound.and.
+        if(real(zk2).le.zkrupbound.and.dimag(zk2).lt.zkiupbound.and.
      1        ilev.gt.ilevcutoff) then
 c             get new pw quadrature
 
@@ -1562,7 +1562,7 @@ subroutine hfmm3dmain(nd,eps,zk,
           deallocate(pgboxwexp)
 
 
-        else if((real(zk2).gt.zkrupbound.or.imag(zk2).gt.zkiupbound).
+        else if((real(zk2).gt.zkrupbound.or.dimag(zk2).gt.zkiupbound).
      1            and.ilev.gt.ilevcutoff) then
           nquad2 = nterms(ilev)*2.2
           if(ifprint.ge.1) print *, "In point and shoot regime"
@@ -1719,7 +1719,7 @@ subroutine hfmm3dmain(nd,eps,zk,
       if(ifcharge.eq.1.and.ifdipole.eq.0) then
         do ilev=1,nlevels
           zk2 = zk*boxsize(ilev)
-          if((real(zk2).gt.zkrupbound.or.imag(zk2).gt.zkiupbound).
+          if((real(zk2).gt.zkrupbound.or.dimag(zk2).gt.zkiupbound).
      1            and.ilev.gt.ilevcutoff) then
 
 C$OMP PARALLEL DO DEFAULT(SHARED)
@@ -1754,7 +1754,7 @@ subroutine hfmm3dmain(nd,eps,zk,
       if(ifcharge.eq.0.and.ifdipole.eq.1) then
         do ilev=1,nlevels
           zk2 = zk*boxsize(ilev)
-          if((real(zk2).gt.zkrupbound.or.imag(zk2).gt.zkiupbound).
+          if((real(zk2).gt.zkrupbound.or.dimag(zk2).gt.zkiupbound).
      1            and.ilev.gt.ilevcutoff) then
 
 C$OMP PARALLEL DO DEFAULT(SHARED)
@@ -1789,7 +1789,7 @@ subroutine hfmm3dmain(nd,eps,zk,
       if(ifcharge.eq.1.and.ifdipole.eq.1) then
         do ilev=1,nlevels
           zk2 = zk*boxsize(ilev)
-          if((real(zk2).gt.zkrupbound.or.imag(zk2).gt.zkiupbound).
+          if((real(zk2).gt.zkrupbound.or.dimag(zk2).gt.zkiupbound).
      1            and.ilev.gt.ilevcutoff) then
 
 C$OMP PARALLEL DO DEFAULT(SHARED)

src/Helmholtz/hfmm3d_memest.f

@@ -377,7 +377,7 @@ subroutine hfmm3d_memest(nd,eps,zk,nsource,source,ifcharge,
 
       zkiupbound = 12*pi
       zkrupbound = 16*pi
-      zi = imag(zkfmm)
+      zi = dimag(zkfmm)
 
       ilevcutoff = -1
 
@@ -391,7 +391,7 @@ subroutine hfmm3d_memest(nd,eps,zk,nsource,source,ifcharge,
 
       do ilev=2,nlevels
         zk2 = zkfmm*boxsize(ilev)
-        if(real(zk2).le.zkrupbound.and.imag(zk2).lt.zkiupbound.and.
+        if(real(zk2).le.zkrupbound.and.dimag(zk2).lt.zkiupbound.and.
      1     ilev.gt.ilevcutoff) then
 
           ier = 0

src/Helmholtz/hfmm3d_mps.f90

@@ -393,7 +393,7 @@ subroutine hfmm3dmain_mps(nd, eps, zk, &
   integer :: impolesort(nmpole)
 
   ! storage stuff for tree and multipole expansions
-  integer :: lmptemp
+  integer :: lmptemp,nboxes
   integer *8 :: iaddr(2,nboxes), lmptot
   double precision :: rmlexp(lmptot)
   double precision :: mptemp(lmptemp)
@@ -406,7 +406,6 @@ subroutine hfmm3dmain_mps(nd, eps, zk, &
   integer :: nterms(0:nlevels)
   integer *8 :: ipointer(8)
   integer :: itree(ltree)
-  integer :: nboxes
   integer :: mnbors,mnlist1, mnlist2,mnlist3,mnlist4
   integer :: isrcse(2,nmpole)
   integer, allocatable :: nlist1(:),list1(:,:)
@@ -819,7 +818,7 @@ subroutine hfmm3dmain_mps(nd, eps, zk, &
 
       ! load the necessary quadrature for plane waves
       zk2 = zk*boxsize(ilev)
-      if ( (real(zk2).le.16*pi) .and. (imag(zk2).le.12*pi) &
+      if ( (real(zk2).le.16*pi) .and. (dimag(zk2).le.12*pi) &
           .and. (ifmp .eq. 0) ) then
         ier = 0
 

src/Helmholtz/hfmm3dwrap_legacy.f

@@ -195,12 +195,12 @@ subroutine hfmm3dparttarg(ier,iprec,zk,nsource,source,
       double complex charge(nsource),dipstr(nsource)
       double precision dipvec(3,nsource)
 
+      integer ntarg
       integer ifpot,iffld,ifpottarg,iffldtarg
       double complex  pot(nsource),fld(3,nsource)
       double complex  pottarg(ntarg),fldtarg(3,ntarg)
 
       integer nd,ifpgh,ifpghtarg
-      integer ntarg
       double precision targ(3,ntarg)
       double complex, allocatable :: dipvec_in(:,:)
       double complex, allocatable :: pottmp(:),gradtmp(:,:)