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! 9*0ccf82acSMartin Diehlpure subroutine FortranSolveAIJ(n,x,ai,aj,adiag,aa,b) 10*0ccf82acSMartin Diehl implicit none (type, external) 11*0ccf82acSMartin Diehl PetscScalar, intent(in) :: aa(0:*),b(0:*) 12*0ccf82acSMartin Diehl PetscInt, intent(in) :: n,ai(0:*), aj(0:*),adiag(0:*) 13*0ccf82acSMartin Diehl PetscScalar, intent(inout) :: x(0:*) 14c96caaccSSatish Balay 15c96caaccSSatish Balay PetscInt i,j,jstart,jend 16c96caaccSSatish Balay PetscScalar sum 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 24c96caaccSSatish Balay sum = b(i) 250113e719SMartin Diehl do j=jstart,jend 26c96caaccSSatish Balay sum = sum - aa(j) * x(aj(j)) 270113e719SMartin Diehl end do 28c96caaccSSatish Balay x(i) = sum 290113e719SMartin Diehl end do 30c96caaccSSatish Balay ! 31c96caaccSSatish Balay ! Backward solve the upper triangular 32c96caaccSSatish Balay ! 330113e719SMartin Diehl do i=n-1,0,-1 34c96caaccSSatish Balay jstart = adiag(i) + 1 35c96caaccSSatish Balay jend = ai(i+1) - 1 36c96caaccSSatish Balay sum = x(i) 370113e719SMartin Diehl do j=jstart,jend 38c96caaccSSatish Balay sum = sum - aa(j)* x(aj(j)) 390113e719SMartin Diehl end do 40c96caaccSSatish Balay x(i) = sum * aa(adiag(i)) 410113e719SMartin Diehl end do 420113e719SMartin Diehlend subroutine FortranSolveAIJ 43