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! 9pure subroutine FortranSolveAIJ(n,x,ai,aj,adiag,aa,b) 10 implicit none (type, external) 11 PetscScalar, intent(in) :: aa(0:*),b(0:*) 12 PetscInt, intent(in) :: n,ai(0:*), aj(0:*),adiag(0:*) 13 PetscScalar, intent(inout) :: x(0:*) 14 15 PetscInt i,jstart,jend 16 ! 17 ! Forward Solve 18 ! 19 x(0) = b(0) 20 do i=1,n-1 21 jstart = ai(i) 22 jend = adiag(i) - 1 23 x(i) = b(i) - sum(aa(jstart:jend)*x(aj(jstart:jend))) 24 end do 25 ! 26 ! Backward solve the upper triangular 27 ! 28 do i=n-1,0,-1 29 jstart = adiag(i) + 1 30 jend = ai(i+1) - 1 31 x(i) = x(i) - sum(aa(jstart:jend)*x(aj(jstart:jend))) * aa(adiag(i)) 32 end do 33end subroutine FortranSolveAIJ 34