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) 10*fe66ebccSMartin Diehl use, intrinsic :: ISO_C_binding 110ccf82acSMartin Diehl implicit none(type, external) 120ccf82acSMartin Diehl PetscScalar, intent(in) :: aa(0:*), b(0:*) 130ccf82acSMartin Diehl PetscInt, intent(in) :: n, ai(0:*), aj(0:*), adiag(0:*) 140ccf82acSMartin Diehl PetscScalar, intent(inout) :: x(0:*) 15c96caaccSSatish Balay 16d66e387eSMartin Diehl PetscInt i, jstart, jend 17c96caaccSSatish Balay ! 18c96caaccSSatish Balay ! Forward Solve 19c96caaccSSatish Balay ! 20c96caaccSSatish Balay x(0) = b(0) 210113e719SMartin Diehl do i = 1, n - 1 22c96caaccSSatish Balay jstart = ai(i) 23c96caaccSSatish Balay jend = adiag(i) - 1 24d66e387eSMartin Diehl x(i) = b(i) - sum(aa(jstart:jend)*x(aj(jstart:jend))) 250113e719SMartin Diehl end do 26c96caaccSSatish Balay ! 27c96caaccSSatish Balay ! Backward solve the upper triangular 28c96caaccSSatish Balay ! 290113e719SMartin Diehl do i = n - 1, 0, -1 30c96caaccSSatish Balay jstart = adiag(i) + 1 31c96caaccSSatish Balay jend = ai(i + 1) - 1 32d66e387eSMartin Diehl x(i) = x(i) - sum(aa(jstart:jend)*x(aj(jstart:jend)))*aa(adiag(i)) 330113e719SMartin Diehl end do 340113e719SMartin Diehlend subroutine FortranSolveAIJ 35