- Fixes redundant IFAIL errors in IGT routines (see #252).

- Removes writing to stdout from Fortran files - this removes dependence on quadmath library (was causing problems in some compilation environments). Now ifail is checked in "interaction.c" - errors should appear only seldom (i.e. when `-eps 12` is used).

Author: myurkin

src/fort/propaesplibreintadda.f

@@ -13,15 +13,15 @@
 c     license: GNU GPL
 
       subroutine propaespacelibreintadda(Rij,k0a,
-     $      gridspacex,gridspacey,gridspacez,relreq,result)
+     $      gridspacex,gridspacey,gridspacez,relreq,result,IFAIL)
       implicit none
       integer i,j
 c     definition of the position of the dipole, observation, wavenumber
 c     ,wavelength, spacing lattice
       double precision result(12)
       double precision k0a,gridspacexm,gridspaceym,gridspacezm
       double precision x,y,z,gridspacex,gridspacey,gridspacez
-      double precision k0,xx0,yy0,zz0
+      double precision k0,xx0,yy0,zz0,RR
       double precision Rij(3)
 
 c     The structure of the result is the following:
@@ -43,23 +43,33 @@ subroutine propaespacelibreintadda(Rij,k0a,
       y=Rij(2)
       z=Rij(3)
       k0=k0a
-      gridspacexm=gridspacex*0.1d0
-      gridspaceym=gridspacey*0.1d0
-      gridspacezm=gridspacez*0.1d0
+      gridspacexm=gridspacex*0.01d0
+      gridspaceym=gridspacey*0.01d0
+      gridspacezm=gridspacez*0.01d0
 c     We perform the integration of the tensor
 c     definition for the integration
       MINCLS = 1000
       MAXCLS = 1000000
       KEY = 0
-      ABSREQ = 0.d0
-      
+
+c     Based on the conservative estimate from below of the largest
+c     component of the integral. For moderate and large distances the
+c     integral will actually scale as k^2/R, for small distances R^-3
+c     dependence will lead to larger values than based on its center
+      RR=x*x+y*y+z*z
+      ABSREQ = RELREQ*gridspacex*gridspacey*gridspacez/(RR**1.5d0)
+c      ABSREQ=0
+
       A(1)=-gridspacex/2.d0
       A(2)=-gridspacey/2.d0
       A(3)=-gridspacez/2.d0
       B(1)=+gridspacex/2.d0
       B(2)=+gridspacey/2.d0
       B(3)=+gridspacez/2.d0
 
+c     If some of the coordinates are small, the corresponding 
+c     non-diagonal terms of the Green's tensor are further considered
+c     null (during each integrand evaluation)
       xx0=1.d0
       yy0=1.d0
       zz0=1.d0
@@ -80,9 +90,12 @@ subroutine propaespacelibreintadda(Rij,k0a,
          result(N)=result(N)/gridspacex/gridspacey/gridspacez
       enddo
 
-      if (ifail.ne.0) then
-         write(*,*) 'IFAIL in IGT routine',IFAIL
-      endif
+c     We return IFAIL instead of testing it here to avoid output to
+c     terminal directly from Fortran, which leads to dependence on
+c     quadmath library
+c     if (ifail.ne.0) then
+c        write(*,*) 'IFAIL in IGT routine',IFAIL
+c     endif
 
       end
 c*************************************************************

src/interaction.c

@@ -68,7 +68,7 @@ static __m128d exptbl[361];
 #ifndef NO_FORTRAN
 // fort/propaesplibreintadda.f
 void propaespacelibreintadda_(const double *Rij,const double *ka,const double *gridspacex,const double *gridspacey,
-	const double *gridspacez,const double *relreq,double *result);
+	const double *gridspacez,const double *relreq,double *result,int *ifail);
 #endif
 // sinint.c
 void cisi(double x,double *ci,double *si);
@@ -974,17 +974,23 @@ static inline void InterTerm_igt(double qvec[static 3],doublecomplex result[stat
 	double rr,rn,invr3,kr,kr2; // |R|, |R/d|, |R|^-3, kR, (kR)^2
 	doublecomplex expval; // exp(ikR)/|R|^3
 	double tmp[12];
-	int comp;
+	int comp,ifail;
 
 	// the following looks complicated, but should be easy to optimize by compiler
 	if (igt_lim==UNDEF || DotProd(qvec,qvec)<=maxRectScale*maxRectScale*igt_lim*igt_lim
 		*(unitsGrid ? 1 : (gridspace*gridspace)) ) {
 		if (unitsGrid) vMultScal(gridspace,qvec,qvec);
-		/* passing complex vectors from Fortran to c is not necessarily portable (at least requires extra effort in
+		/* passing complex vectors from Fortran to C is not necessarily portable (at least requires extra effort in
 		 * the Fortran code. So we do it through double. This is not bad for performance, since double is anyway used
 		 * internally for integration in this Fortran routine.
 		 */
-		propaespacelibreintadda_(qvec,&WaveNum,&gridSpaceX,&gridSpaceY,&gridSpaceZ,&igt_eps,tmp);
+		propaespacelibreintadda_(qvec,&WaveNum,&gridSpaceX,&gridSpaceY,&gridSpaceZ,&igt_eps,tmp,&ifail);
+		if (ifail!=0) {
+			if (ifail==1) LogWarning(EC_WARN,ALL_POS,"Failed to reach relative accuracy of %g for Green's tensor "
+				"integration for distance "GFORMDEF3V,igt_eps,COMP3V(qvec));
+			else LogWarning(EC_WARN,ALL_POS,"Error in Green's tensor integration (code %d - see dcuhre.f) for distance "
+				GFORMDEF3V,ifail,COMP3V(qvec));
+		}
 		for (comp=0;comp<6;comp++) result[comp] = tmp[comp] + I*tmp[comp+6];
 	}
 	else {