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! 9c96caaccSSatish Balaysubroutine FortranSolveAIJ(n,x,ai,aj,adiag,aa,b) 10c96caaccSSatish Balay implicit none 11c96caaccSSatish Balay PetscScalar x(0:*),aa(0:*),b(0:*) 12c96caaccSSatish Balay PetscInt n,ai(0:*) 13c96caaccSSatish Balay PetscInt aj(0:*),adiag(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) 21*0113e719SMartin Diehl do i=1,n-1 22c96caaccSSatish Balay jstart = ai(i) 23c96caaccSSatish Balay jend = adiag(i) - 1 24c96caaccSSatish Balay sum = b(i) 25*0113e719SMartin Diehl do j=jstart,jend 26c96caaccSSatish Balay sum = sum - aa(j) * x(aj(j)) 27*0113e719SMartin Diehl end do 28c96caaccSSatish Balay x(i) = sum 29*0113e719SMartin Diehl end do 30c96caaccSSatish Balay ! 31c96caaccSSatish Balay ! Backward solve the upper triangular 32c96caaccSSatish Balay ! 33*0113e719SMartin 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) 37*0113e719SMartin Diehl do j=jstart,jend 38c96caaccSSatish Balay sum = sum - aa(j)* x(aj(j)) 39*0113e719SMartin Diehl end do 40c96caaccSSatish Balay x(i) = sum * aa(adiag(i)) 41*0113e719SMartin Diehl end do 42*0113e719SMartin Diehl 43*0113e719SMartin Diehlend subroutine FortranSolveAIJ 44