ivan-pi / y12m

Solution of Large and Sparse Systems of Linear Algebraic Equations
16 stars 1 forks source link

Tricky loop refactoring in y12mb #8

Open ivan-pi opened 11 months ago

ivan-pi commented 11 months ago

https://github.com/ivan-pi/y12m/blob/8a4b83a909c1f97af27fcf3ee3f1ae1f01795c15/src/legacy/y12mbe.f#L93-L114

This loop has some subtle checks inside; the order appears to matter. Only if all three conditions are false, it will signal an error:

    if(r.eq.0)go to 72 
    if(j.eq.l1)go to 72 
    if(rnr(index-1).ne.i)go to 72 
    ifail=11 
    go to 22 
 72 continue

When refactored as follows:

         if(r /= 0 .and. j /= l1 .and. rnr(index-1) == i) then
               ifail = 11
               return
         end if

it leads to an out-of-bounds error:

~/fortran/y12m/build$ ./test_y12m
At line 158 of file /home/ivan/fortran/y12m/src/y12mbe.f90
Fortran runtime error: Index '0' of dimension 1 of array 'rnr' below lower bound of 1

⚠️ I've slightly changed the labels compared to the original version on Netlib - https://www.netlib.org/y12m/:

C y12mbe.f, line ??
      do 70 i=1,n
      if(mode.eq.2)go to 60
      if(ha(i,2).eq.0)go to 60
      ha(i,11)=l4
      l4=l4+ha(i,2)
      ha(i,2)=ha(i,11)
   60 ha(i,3)=ha(i,1)+ha(i,3)-1
      l1=ha(i,1)
      l2=ha(i,3)
      do 70 j=l1,l2
      l3=snr(j)
      r=ha(l3,6)
      index=ha(l3,4)+r
      rnr(index)=i
      if(r.eq.0)go to 70
      if(j.eq.l1)go to 70
      if(rnr(index-1).ne.i)go to 70
      ifail=11
      go to 22
   70 ha(l3,6)=r+1