1c96caaccSSatish Balay! 2c96caaccSSatish Balay! 3c96caaccSSatish Balay! Fortran kernel for sparse triangular solve in the AIJ matrix format 4c96caaccSSatish Balay! This ONLY works for factorizations in the NATURAL ORDERING, i.e. 5c96caaccSSatish Balay! with MatSolve_SeqAIJ_NaturalOrdering() 6c96caaccSSatish Balay! 7c96caaccSSatish Balay#include <petsc/finclude/petscsys.h> 8c96caaccSSatish Balay! 90ccf82acSMartin Diehlpure subroutine FortranSolveAIJ(n,x,ai,aj,adiag,aa,b) 100ccf82acSMartin Diehl implicit none (type, external) 110ccf82acSMartin Diehl PetscScalar, intent(in) :: aa(0:*),b(0:*) 120ccf82acSMartin Diehl PetscInt, intent(in) :: n,ai(0:*), aj(0:*),adiag(0:*) 130ccf82acSMartin Diehl PetscScalar, intent(inout) :: x(0:*) 14c96caaccSSatish Balay 15*d66e387eSMartin Diehl PetscInt i,jstart,jend 16c96caaccSSatish Balay ! 17c96caaccSSatish Balay ! Forward Solve 18c96caaccSSatish Balay ! 19c96caaccSSatish Balay x(0) = b(0) 200113e719SMartin Diehl do i=1,n-1 21c96caaccSSatish Balay jstart = ai(i) 22c96caaccSSatish Balay jend = adiag(i) - 1 23*d66e387eSMartin Diehl x(i) = b(i) - sum(aa(jstart:jend)*x(aj(jstart:jend))) 240113e719SMartin Diehl end do 25c96caaccSSatish Balay ! 26c96caaccSSatish Balay ! Backward solve the upper triangular 27c96caaccSSatish Balay ! 280113e719SMartin Diehl do i=n-1,0,-1 29c96caaccSSatish Balay jstart = adiag(i) + 1 30c96caaccSSatish Balay jend = ai(i+1) - 1 31*d66e387eSMartin Diehl x(i) = x(i) - sum(aa(jstart:jend)*x(aj(jstart:jend))) * aa(adiag(i)) 320113e719SMartin Diehl end do 330113e719SMartin Diehlend subroutine FortranSolveAIJ 34