Code won't run on fotran mac; decimal/0/infinity issue

[Krotan]

Hi, I have a fortran code which I've successfully compiled and run on both windows PC's and PPC macs. However, the same code will compile but run incorrectly on the new intel macs. I have tried g77, gfortran, and ifort, but they all produce the same problem.

One subroutine in particular calculates and manipulates a 2D array, adjusting cells based on a previous cells value. However, it seems to be having a problem where all of the values are 0 or infinity (mostly from dividing by those zeroes) but I can't figure out why. Would you mind taking a peek at the code and looking for an obvious error? We've narrowed it down to this subroutine. The array 'a' is the problem. Thank you.

subroutine gaussj(a,n,np,b,m,mp)
implicit double precision (a-h,o-z)
parameter (nmax=100)
dimension a(np,np),b(np,mp),ipiv(nmax),indxr(nmax),indxc(nmax)
dum=0
do 5 l=1,np
do 6 m=1,np
a(l,m)=100
C print *,'l=',l,'m=',m,'a(l,m)=',a(l,m)
6 continue
5 continue
C return
do 11 j=1,n
ipiv(j)=0
11 continue
do 22 i=1,n
big=0.d0
do 13 j=1,n
if(ipiv(j).ne.1)then
do 12 k=1,n
if (ipiv(k).eq.0) then
if (dabs(a(j,k)).ge.big)then
big=dabs(a(j,k))
irow=j
icol=k
C print *,'icol=',icol,'a(j,k)=',a(j,k)
endif
else if (ipiv(k).gt.1) then
write(20,*)'ERROR(GAUSS): singular matrix-1'
C print *,-1,'n=',n,'np=',np,'m=',m,'mp=',mp,'k=',k
m=-1
return
endif
12 continue
endif
13 continue
ipiv(icol)=ipiv(icol)+1
if (irow.ne.icol) then
do 14 l=1,n
dum=a(irow,l)
a(irow,l)=a(icol,l)
a(icol,l)=dum
print *,dum
14 continue
do 15 l=1,m
dum=b(irow,l)
b(irow,l)=b(icol,l)
b(icol,l)=dum
15 continue
endif
indxr(i)=irow
indxc(i)=icol
print *,'icol=',icol,'k=',k,'a(icol,icol)=',a(icol,icol)
if (a(icol,icol).eq.0.) then
write(20,*)'ERROR(GAUSS): singular matrix-2'
print *, -2,'n=',n,'np=',np,'m=',m,'mp=',mp,'k=',k,'icol=',icol
print *,'a(icol,icol)=',a(icol,icol)
m=-1
return
endif
pivinv=1.d0/a(icol,icol)
a(icol,icol)=1.d0
do 16 l=1,n
a(icol,l)=a(icol,l)*pivinv
16 continue
do 17 l=1,m
b(icol,l)=b(icol,l)*pivinv
17 continue
do 21 ll=1,n
if(ll.ne.icol)then
dum=a(ll,icol)
a(ll,icol)=0.d0
do 18 l=1,n
a(ll,l)=a(ll,l)-a(icol,l)*dum
18 continue
do 19 l=1,m
b(ll,l)=b(ll,l)-b(icol,l)*dum
19 continue
endif
21 continue
22 continue
do 24 l=n,1,-1
if(indxr(l).ne.indxc(l))then
do 23 k=1,n
dum=a(k,indxr(l))
a(k,indxr(l))=a(k,indxc(l))
a(k,indxc(l))=dum
23 continue
endif
24 continue
return
end

 

[gullipe]

Hi Krotan

Try to change the following line (the 54th line):
if (a(icol,icol).eq.0.) then
to:
if (a(icol,icol).eq.0.d0) then

I do not know if this solves the problem, but it is worth trying.

 

[FJacq]

Try to change the following line (the 54th line):
if (a(icol,icol).eq.0.) then
to:
if (a(icol,icol).eq.0.d0) then

I do not know if this solves the problem, but it is worth trying.

No, sorry. This will not solve the problem because 0 or 0.d0 are strictly equivalent in that context. It would be better to test something like :

      if (abs(a(icol,icol)).le.1.d-30) then

Anyway, the program seems very strange : you are testing a subroutine named GAUSS, which probably tries to solve a linear system of equations (I vaguely recognize a GAUSS algorithm with a total pivot search). But you initialize all the terms of that matrix to 100 : a very very bag idea because such matrix cannot be inversed !

Now, I suspect that the initialization part you showed did not exist when the trouble was observed the first time (the matrix was certainly initialized outside the routine GAUSS). From my point of view, the original problem came from something else. Are you sure that the matrix A was exactly the same when you executed your program on PC-Windows and on MAC ?

I suggest to replace the wrong initialization phase by simply printing out that matrix at the entrance of the subroutine using a fixed format like "(10(1X,ES14.7))" for instance. Then run the code on Windows and Mac and compare the two listing.

Other question : is the argument N less than 100 ? I hope so else your local arrays ipiv(nmax),indxr(nmax),indxc(nmax) are too small.

Notice that you can replace NMAX by N : all the Fortran compilers are F95 now and automatic arrays are accepted.

François Jacq