1*c96caaccSSatish Balay! 2*c96caaccSSatish Balay! 3*c96caaccSSatish Balay! Fortran kernel for sparse triangular solve in the AIJ matrix format 4*c96caaccSSatish Balay! This ONLY works for factorizations in the NATURAL ORDERING, i.e. 5*c96caaccSSatish Balay! with MatSolve_SeqAIJ_NaturalOrdering() 6*c96caaccSSatish Balay! 7*c96caaccSSatish Balay#include <petsc/finclude/petscsys.h> 8*c96caaccSSatish Balay! 9*c96caaccSSatish Balay subroutine FortranSolveAIJ(n,x,ai,aj,adiag,aa,b) 10*c96caaccSSatish Balay implicit none 11*c96caaccSSatish Balay PetscScalar x(0:*),aa(0:*),b(0:*) 12*c96caaccSSatish Balay PetscInt n,ai(0:*) 13*c96caaccSSatish Balay PetscInt aj(0:*),adiag(0:*) 14*c96caaccSSatish Balay 15*c96caaccSSatish Balay PetscInt i,j,jstart,jend 16*c96caaccSSatish Balay PetscScalar sum 17*c96caaccSSatish Balay! 18*c96caaccSSatish Balay! Forward Solve 19*c96caaccSSatish Balay! 20*c96caaccSSatish Balay x(0) = b(0) 21*c96caaccSSatish Balay do 20 i=1,n-1 22*c96caaccSSatish Balay jstart = ai(i) 23*c96caaccSSatish Balay jend = adiag(i) - 1 24*c96caaccSSatish Balay sum = b(i) 25*c96caaccSSatish Balay do 30 j=jstart,jend 26*c96caaccSSatish Balay sum = sum - aa(j) * x(aj(j)) 27*c96caaccSSatish Balay 30 continue 28*c96caaccSSatish Balay x(i) = sum 29*c96caaccSSatish Balay 20 continue 30*c96caaccSSatish Balay 31*c96caaccSSatish Balay! 32*c96caaccSSatish Balay! Backward solve the upper triangular 33*c96caaccSSatish Balay! 34*c96caaccSSatish Balay do 40 i=n-1,0,-1 35*c96caaccSSatish Balay jstart = adiag(i) + 1 36*c96caaccSSatish Balay jend = ai(i+1) - 1 37*c96caaccSSatish Balay sum = x(i) 38*c96caaccSSatish Balay do 50 j=jstart,jend 39*c96caaccSSatish Balay sum = sum - aa(j)* x(aj(j)) 40*c96caaccSSatish Balay 50 continue 41*c96caaccSSatish Balay x(i) = sum * aa(adiag(i)) 42*c96caaccSSatish Balay 40 continue 43*c96caaccSSatish Balay return 44*c96caaccSSatish Balay end 45*c96caaccSSatish Balay 46