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