program domain 
      dimension u(9,9),f(9,9) 
! $SPF ANALYSIS (PRIVATE(gran1))       
      dimension u1(0:10,0:10),f1(9,9),gran1(9) 
! $SPF ANALYSIS (PRIVATE(gran2))       
      dimension u2(0:10,0:10),f2(9,9),gran2(9) 
! $SPF ANALYSIS (PRIVATE(g))       
      dimension a(9,9),g(9) 
      net = 10 
      ntime = 2 
      n = net - 1 
      m = net / 2 
      h = 1. / net 
      r = 0.125 
      tau = r * (h * h) 
      r1 = 1 - 4 * r 
      rr = 1 + 2 * r 
      r2 = 1 - 2 * r 
! $SPF TRANSFORM(UNROLL)
      do  i = 1,n 
         x = i * h 
         s = x * (1 - x) 
         do  j = 1,n 
            y = j * h 
            q = y * (1 - y) 
            u(i,j) = 16 * s * q 
10          f(i,j) = (s + q) * 32 
         enddo   
      enddo
!$SPF TRANSFORM(FISSION(j))   
      do  j = 1,n 
         gran1(j) = u(m,j) 
         gran2(j) = u(m + 1,j) 
      enddo 
!$SPF TRANSFORM(FISSION(i))  
      do  i = 0,m 
         u1(i,0) = 0 
         u1(i,n + 1) = 0 
      enddo
!$SPF TRANSFORM(FISSION(i))   
      do  i = 1,m 
         u2(i,0) = 0 
         u2(i,n + 1) = 0 
      enddo 
!$SPF TRANSFORM(FISSION(j))  
      do  j = 1,n 
         u1(0,j) = 0 
         u1(m + 1,j) = gran2(j) 
         u2(0,j) = gran1(j) 
         u2(m,j) = 0 
      enddo 
!$SPF TRANSFORM(FISSION(i,j))   
      do  i = 1,m 
         do  j = 1,n 
            u1(i,j) = u(i,j) 
60          f1(i,j) = f(i,j) 
         enddo   
      enddo   
      do  i = 1,m - 1 
         do  j = 1,n 
            u2(i,j) = u(i + m,j) 
70          f2(i,j) = f(i + m,j) 
         enddo   
      enddo   
      do  ktau = 1,ntime 
         do  i = 1,m - 1 
            do  j = 1,n 
               w = u1(i - 1,j) + u1(i + 1,j) + u1(i,j - 1) + u1(i,j + 1) 
80             a(i,j) = r * w + r1 * u1(i,j) + tau * f1(i,j) 
            enddo   
         enddo   
         do  j = 1,n 
            w = u1(m - 1,j) + gran2(j) + u1(m,j - 1) + u1(m,j + 1) 
            a(m,j) = r * w + r1 * u1(m,j) + tau * f1(m,j) 
         enddo   
         do  i = 1,m 
            do  j = 1,n 
100            u1(i,j) = a(i,j) 
            enddo   
         enddo   
         do  j = 1,n 
             g(j) = r * gran1(j) + r2 * u2(1,j) + r * u2(2,j) +  
     &         tau * f2(1,j) 
         enddo   
         call progon(n,r,rr,r,g) 
         do  j = 1,n 
            u2(1,j) = g(j) 
         enddo   
         do  i = 2,m - 1 
            do  j = 1,n 
               g(j) = r * u2(i - 1,j) + r2 * u2(i,j) + r *  
     &           u2(i + 1,j) + tau * f2(i,j) 
            enddo   
            call progon(n,r,rr,r,g) 
            do  j = 1,n 
140            u2(i,j) = g(j) 
            enddo   
         enddo   
         do  j = 1,n 
            gran1(j) = u1(m,j) 
            gran2(j) = u2(1,j) 
         enddo   
         do  i = 1,m 
            do  j = 1,n 
160            u(i,j) = u1(i,j) 
            enddo   
         enddo   
         do  i = 1,m - 1 
            do  j = 1,n 
170            u(i + m,j) = u2(i,j) 
            enddo   
         enddo   
      enddo   
      print 190, net,ntime 
190    FORMAT(I5,I5) 
      print *, u 
200    FORMAT(F5.2) 
      end 
 
      subroutine progon (n, a, b, c, f) 
      dimension f(n) 
! $SPF ANALYSIS (PRIVATE(alfa,beta))       
 
      real, allocatable:: alfa(:),beta(:) 
      
      allocate(alfa(9),beta(9)) 
      alfa(1) = 0. 
      beta(1) = 0. 
      do  i = 1,n - 1 
         w = 1. / (b - a * alfa(i)) 
         alfa(i + 1) = c * w 
         beta(i + 1) = w * (a * beta(i) + f(i)) 
      enddo   
      f(n) = (a * beta(n) + f(n)) / (b - a * alfa(n)) 
      do  i = n - 1,1,(-(1)) 
         f(i) = alfa(i + 1) * f(i + 1) + beta(i + 1) 
      enddo   
      deallocate(alfa) 
      deallocate(beta) 
      return  
      end