C23456789012345678901234567890123456789012345678901234567890123456789012
C
C         Recurrences using the partition method
C         Marco Zagha (marcoz@cs.cmu.edu)
C         Started 2/22/93
C
C         Note that this stuff only works for exact multiples of the stride
C23456789012345678901234567890123456789012345678901234567890123456789012
      
      real function r1last(d,s,nrows,ncols)
      integer nrows,ncols
      real d(nrows,ncols),s(nrows,ncols)
      integer i,j
      real sum
      
      i = 1
      do j = 1,ncols
         d(i,j) = s(i,j)
      end do
      
      do i = 2,nrows
         do j = 1,ncols
            d(i,j) = d(i-1,j) + s(i,j)
         end do
      end do
      
      sum = 0.0
      do j = 1,ncols
         sum = sum + d(nrows,j)
      end do
      
      r1last = sum
      return
      end
      
C23456789012345678901234567890123456789012345678901234567890123456789012      
      
      subroutine r1all(d,s,nrows,ncols)
      integer nrows,ncols
      real d(nrows,ncols),s(nrows,ncols)
      integer i,j
      
      i = 1
      do j = 1,ncols
         d(i,j) = s(i,j)
      end do
      
      do i = 2,nrows
         do j = 1,ncols
            d(i,j) = d(i-1,j) + s(i,j)
         end do
      end do
      
      do j = 2,ncols
         do i = 1,nrows
            d(i,j) = d(i,j) + d(nrows,j-1)
         end do
      end do
      
      return
      end

C23456789012345678901234567890123456789012345678901234567890123456789012

      real function r1Clast(d,d2,s,b,nrows,ncols)
      integer nrows,ncols
      real d(nrows,ncols),d2(nrows,ncols)
      real s(nrows,ncols),b(nrows,ncols)
      real sum,prod
      integer i,j
      
      i = 1
      do j = 1,ncols
         d(i,j) = 0.0
         d2(i,j) = 1.0
      end do
      
      do i = 2,nrows
         do j = 1,ncols
            d(i,j) = d(i-1,j) * s(i-1,j) + b(i-1,j)
            d2(i,j) = d2(i-1,j) * s(i-1,j)
         end do
      end do
  
      sum = 0.0
      prod = 1.0
      do j = 2,ncols
         sum = (sum * s(nrows,j-1) + b(nrows,j-1)) 
     $          * d2(nrows,j) + d(nrows,j)
      end do
      r1Clast = sum

      return
      end
      
C23456789012345678901234567890123456789012345678901234567890123456789012
      
      subroutine r1Call(d,d2,s,b,nrows,ncols)
      integer nrows,ncols
      real d(nrows,ncols),d2(nrows,ncols)
      real s(nrows,ncols),b(nrows,ncols)
      real temp
      integer i,j
      
      i = 1
      do j = 1,ncols
         d(i,j) = 0.0
         d2(i,j) = 1.0
      end do
      
      do i = 2,nrows
         do j = 1,ncols
            d(i,j) = d(i-1,j) * s(i-1,j) + b(i-1,j)
            d2(i,j) = d2(i-1,j) * s(i-1,j)
         end do
      end do
      
      do j = 2,ncols
         temp = d(nrows,j-1) * s(nrows,j-1) + b(nrows,j-1)
         do i = 1,nrows
           d(i,j) = temp * d2(i,j) + d(i,j)
         end do
      end do
      
      return
      end
      
C23456789012345678901234567890123456789012345678901234567890123456789012
      
      real function r2last(d,d2,s,b,nrows,ncols)
      integer nrows,ncols
      real d(nrows,ncols),d2(nrows,ncols)
      real s(nrows,ncols),b(nrows,ncols)
      real temp, temp2, sum, sum2
      integer i,j
      
      i = 1
      do j = 1,ncols
         d(i,j)    = 0.0
         d(i+1,j)  = 1.0
         d2(i,j)   = 1.0
         d2(i+1,j) = 0.0
      end do
      
      do i = 3,nrows
         do j = 1,ncols
            d(i,j)  =  d(i-1,j) * s(i-2,j) +  d(i-2,j) * b(i-2,j)
            d2(i,j) = d2(i-1,j) * s(i-2,j) + d2(i-2,j) * b(i-2,j)
         end do
      end do

      sum = d(nrows,1)
      sum2 = d(nrows-1,1)
      do j = 2,ncols
         temp=sum*s(nrows-1,j-1)+sum2*b(nrows-1,j-1)
         temp2 = temp * s(nrows,j-1) + sum * b(nrows,j-1)
         sum2 = temp * d2(nrows-1,j) + temp2 * d(nrows-1,j)
         sum  = temp * d2(nrows,j) + temp2 * d(nrows,j)
      end do
      r2last = sum

      return
      end

C23456789012345678901234567890123456789012345678901234567890123456789012

      subroutine r2all(d,d2,s,b,nrows,ncols)
      integer nrows,ncols
      real d(nrows,ncols),d2(nrows,ncols)
      real s(nrows,ncols),b(nrows,ncols)
      real temp, temp2
      integer i,j
      
      i = 1
      do j = 1,ncols
         d(i,j)    = 0.0
         d(i+1,j)  = 1.0
         d2(i,j)   = 1.0
         d2(i+1,j) = 0.0
      end do
      
      do i = 3,nrows
         do j = 1,ncols
            d(i,j)  =  d(i-1,j) * s(i-2,j) +  d(i-2,j) * b(i-2,j)
            d2(i,j) = d2(i-1,j) * s(i-2,j) + d2(i-2,j) * b(i-2,j)
         end do
      end do
      
      do j = 2,ncols
         temp=d(nrows,j-1)*s(nrows-1,j-1)+d(nrows-1,j-1)*b(nrows-1,j-1)
         temp2 = temp * s(nrows,j-1) + d(nrows,j-1) * b(nrows,j-1)
         do i = 1,nrows
           d(i,j) = temp * d2(i,j) + temp2 * d(i,j)
         end do
      end do
      
      return
      end
 
     
