!
!    Fortran kernels for SOR relaxations
!
#include "include/finclude/petscdef.h"
!
      subroutine FortranRelaxAIJForwardZero(n,omega,x,ai,aj,            &
     &           adiag,idiag,aa,b)
      implicit none
      PetscScalar      x(0:*),aa(0:*),b(0:*),idiag(0:*)
      PetscReal        omega
      PetscInt          n,ai(0:*),aj(0:*),adiag(0:*)

      PetscInt          i,j,jstart,jend
      PetscScalar      sum
!     
!     Forward Solve with zero initial guess
!
      x(0) = omega*b(0)*idiag(0)
      do 20 i=1,n-1
         jstart = ai(i)
         jend   = adiag(i) - 1
         sum    = b(i)
         do 30 j=jstart,jend
            sum  = sum -  aa(j) * x(aj(j))
 30      continue
         x(i) = omega*sum*idiag(i)
 20   continue
      
!     return
      end
!
!-------------------------------------------------------------------
!      
      subroutine FortranRelaxAIJBackwardZero(n,omega,x,ai,aj,           &
     &                                  adiag,idiag,aa,b)
      implicit none
      PetscScalar      x(0:*),aa(0:*),b(0:*),idiag(0:*)
      PetscReal        omega
      PetscInt          n,ai(0:*),aj(0:*),adiag(0:*)

      PetscInt          i,j,jstart,jend
      PetscScalar      sum
!
!     Backward solve with zero initial guess
!
      do 40 i=n-1,0,-1
         jstart  = adiag(i) + 1
         jend    = ai(i+1) - 1
         sum     = b(i)
         do 50 j=jstart,jend
            sum = sum - aa(j)* x(aj(j)) 
 50      continue
         x(i)    = omega*sum*idiag(i)
 40   continue
      return
      end
      
!-------------------------------------------------------------------
!      
      subroutine FortranRelaxAIJForward(n,omega,x,ai,aj,adiag,aa,b)
      implicit none
      PetscScalar      x(0:*),aa(0:*),b(0:*)
      PetscReal        omega
      PetscInt          n,ai(0:*),aj(0:*),adiag(0:*)

      PetscInt          i,j,jstart,jend
      PetscScalar      sum
!
!     Forward solve with non-zero initial guess
!
      do 40 i=0,n-1
         sum    = b(i)

         jstart = ai(i)
         jend    = ai(i+1) - 1
         do 50 j=jstart,jend
            sum = sum - aa(j)* x(aj(j)) 
 50      continue
         x(i)    = (1.0 - omega)*x(i) +                                 &
     &              omega*(sum + aa(adiag(i))*x(i))/ aa(adiag(i))
 40   continue
      return
      end
      
!-------------------------------------------------------------------
!      
      subroutine FortranRelaxAIJBackward(n,omega,x,ai,aj,adiag,aa,b)
      implicit none
      PetscScalar      x(0:*),aa(0:*),b(0:*)
      PetscReal        omega
      PetscInt          n,ai(0:*),aj(0:*),adiag(0:*)

      PetscInt          i,j,jstart,jend
      PetscScalar      sum
!
!     Backward solve with non-zero initial guess
!
      do 40 i=n-1,0,-1
         sum    = b(i)

         jstart = ai(i)
         jend   = ai(i+1) - 1
         do 50 j=jstart,jend
            sum = sum - aa(j)* x(aj(j)) 
 50      continue
         x(i)    = (1.0 - omega)*x(i) +                                 &
     &              omega*(sum + aa(adiag(i))*x(i)) / aa(adiag(i))
 40   continue
      return
      end
      
