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:*)
14
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
30
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
45
46