1: !
  2: !    Fortran kernels for SOR relaxations
  3: !
  4: #include <petsc/finclude/petscsysdef.h>
  5: !
  6:       subroutine FortranRelaxAIJForwardZero(n,omega,x,ai,aj,            &
  7:      &           adiag,idiag,aa,b)
  8:       implicit none
  9:       PetscScalar x(0:*),aa(0:*)
 10:       PetscScalar b(0:*),idiag(0:*)
 11:       PetscReal   omega
 12:       PetscInt    n,ai(0:*)
 13:       PetscInt    aj(0:*),adiag(0:*)

 15:       PetscInt    i,j,jstart,jend
 16:       PetscScalar sum
 17: !
 18: !     Forward Solve with zero initial guess
 19: !
 20:       x(0) = omega*b(0)*idiag(0)
 21:       do 20 i=1,n-1
 22:          jstart = ai(i)
 23:          jend   = adiag(i) - 1
 24:          sum    = b(i)
 25:          do 30 j=jstart,jend
 26:             sum  = sum -  aa(j) * x(aj(j))
 27:  30      continue
 28:          x(i) = omega*sum*idiag(i)
 29:  20   continue

 31: !     return
 32:       end
 33: !
 34: !-------------------------------------------------------------------
 35: !
 36:       subroutine FortranRelaxAIJBackwardZero(n,omega,x,ai,aj,           &
 37:      &                                  adiag,idiag,aa,b)
 38:       implicit none
 39:       PetscScalar x(0:*),aa(0:*)
 40:       PetscScalar b(0:*),idiag(0:*)
 41:       PetscReal   omega
 42:       PetscInt    aj(0:*),adiag(0:*)
 43:       PetscInt    n,ai(0:*)

 45:       PetscInt    i,j,jstart,jend
 46:       PetscScalar sum
 47: !
 48: !     Backward solve with zero initial guess
 49: !
 50:       do 40 i=n-1,0,-1
 51:          jstart  = adiag(i) + 1
 52:          jend    = ai(i+1) - 1
 53:          sum     = b(i)
 54:          do 50 j=jstart,jend
 55:             sum = sum - aa(j)* x(aj(j))
 56:  50      continue
 57:          x(i)    = omega*sum*idiag(i)
 58:  40   continue
 59:       return
 60:       end

 62: !-------------------------------------------------------------------
 63: !
 64:       subroutine FortranRelaxAIJForward(n,omega,x,ai,aj,adiag,aa,b)
 65:       implicit none
 66:       PetscScalar x(0:*),aa(0:*),b(0:*)
 67:       PetscReal   omega
 68:       PetscInt    n,ai(0:*)
 69:       PetscInt    aj(0:*),adiag(0:*)

 71:       PetscInt    i,j,jstart,jend
 72:       PetscScalar sum
 73: !
 74: !     Forward solve with non-zero initial guess
 75: !
 76:       do 40 i=0,n-1
 77:          sum    = b(i)

 79:          jstart = ai(i)
 80:          jend    = ai(i+1) - 1
 81:          do 50 j=jstart,jend
 82:             sum = sum - aa(j)* x(aj(j))
 83:  50      continue
 84:          x(i)    = (1.0 - omega)*x(i) +                                 &
 85:      &              omega*(sum + aa(adiag(i))*x(i))/ aa(adiag(i))
 86:  40   continue
 87:       return
 88:       end

 90: !-------------------------------------------------------------------
 91: !
 92:       subroutine FortranRelaxAIJBackward(n,omega,x,ai,aj,adiag,aa,b)
 93:       implicit none
 94:       PetscScalar x(0:*),aa(0:*),b(0:*)
 95:       PetscReal   omega
 96:       PetscInt    n,ai(0:*)
 97:       PetscInt    aj(0:*),adiag(0:*)

 99:       PetscInt    i,j,jstart,jend
100:       PetscScalar sum
101: !
102: !     Backward solve with non-zero initial guess
103: !
104:       do 40 i=n-1,0,-1
105:          sum    = b(i)

107:          jstart = ai(i)
108:          jend   = ai(i+1) - 1
109:          do 50 j=jstart,jend
110:             sum = sum - aa(j)* x(aj(j))
111:  50      continue
112:          x(i)    = (1.0 - omega)*x(i) +                                 &
113:      &              omega*(sum + aa(adiag(i))*x(i)) / aa(adiag(i))
114:  40   continue
115:       return
116:       end