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