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