Data Local Iterative Methods For The Efficient Solution of Partial Differential Equations

A cooperation
between

and
.

Funded by
.

rb8.F Program Source Code

Fortran77 Source:

define(rel,$3=f($2,$1)+ce*u($2,$1)+no*u($2-1,$1)+
     &      so*u($2+1,$1)+ea*u($2,$1+1)+we*u($2,$1-1)
      u($2,$1)= $3
      ops= ops+1)
define(Crel,)
      program rb
      implicit none
      integer l,n,i,j,k,it,loops,m
      parameter (n=NMAX, m=NMELT, loops=2**28/n**2*m)
      double precision u(0:n+3,0:n), f(0:n+3,0:n),
     &  no,we,so,ea,ce,
     &  r, omega,
     &  dnull
      integer null
      integer ops
C
      read *, null
      dnull= null
      no= 0.25d0+dnull
      we= 0.25d0+dnull*dnull
      ea= 0.25d0-dnull
      so= 0.25d0+dnull+dnull*dnull
      ce= dnull
      omega= 1.5d0
      do j=0,n
      do i=0,n
         u(i,j)= 0d0
         f(i,j)= 0d0
      enddo
      enddo
      do i=0,n
         u(i,n)= 1d0
      enddo
C redBlackRelax8
      do it=0,loops/2-1
C
         do l=2,1,-1
            do j=1,n-1,2
               rel(1,j,r)
            enddo
            do i=2,l*2-1
               do j=1+mod(i+1,2),n-1,2
                  rel(i,j,r)
                  rel(i-1,j,r)
               enddo
            enddo
         enddo
C
         do i=2*2,n-1
            do j=1+mod(i+1,2),n-2-1,4
               do l=0,2*2-1,2
                  k=i-l
                  rel(k,j,r)
                  rel(k,j+2,r)
                  rel(k-1,j,r)
                  rel(k-1,j+2,r)
               enddo
            enddo
            do l=0,2*2-1,2
               k=i-l
               rel(k,n-2,r)
               rel(k-1,n-2,r)
            enddo
            i=i+1
            do j=1+mod(i+1,2),n-2-1,4
               do l=0,2*2-1,2
                  k=i-l
                  rel(k,j,r)
                  rel(k,j+2,r)
                  rel(k-1,j,r)
                  rel(k-1,j+2,r)
               enddo
            enddo
         enddo
C
         do l=1,2
            do i=(l-1)*2,1,-1
               do j=1+mod(i+1,2),n-1,2
                  rel(n-i,j,r)
                  rel(n-i-1,j,r)
               enddo
            enddo
            do j=2,n-1,2
               rel(n-1,j,r)
            enddo
         enddo
      enddo
      
      print *, u(n/2,n/2), ' SIZE', n, ' Operations=', ops
      end