1: !
  2: !
  3: !    Fortran kernel for sparse triangular solve in the AIJ matrix format
  4: ! This ONLY works for factorizations in the NATURAL ORDERING, i.e.
  5: ! with MatSolve_SeqAIJ_NaturalOrdering()
  6: !
  7: #include <petsc/finclude/petscsys.h>
  8: !
  9:       subroutine FortranSolveAIJ(n,x,ai,aj,adiag,aa,b)
 10:       implicit none
 11:       PetscScalar x(0:*),aa(0:*),b(0:*)
 12:       PetscInt n,ai(0:*)
 13:       PetscInt aj(0:*),adiag(0:*)

 15:       PetscInt i,j,jstart,jend
 16:       PetscScalar sum
 17: !
 18: !     Forward Solve
 19: !
 20:       x(0) = b(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) = sum
 29:  20   continue

 31: !
 32: !     Backward solve the upper triangular
 33: !
 34:       do 40 i=n-1,0,-1
 35:          jstart  = adiag(i) + 1
 36:          jend    = ai(i+1) - 1
 37:          sum     = x(i)
 38:          do 50 j=jstart,jend
 39:             sum = sum - aa(j)* x(aj(j))
 40:  50      continue
 41:          x(i)    = sum * aa(adiag(i))
 42:  40   continue
 43:       return
 44:       end