Actual source code: fsolve.F90

      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: !
      9: pure subroutine FortranSolveAIJ(n, x, ai, aj, adiag, aa, b)
     10:   use, intrinsic :: ISO_C_binding
     11:   implicit none(type, external)
     12:   PetscScalar, intent(in) :: aa(0:*), b(0:*)
     13:   PetscInt, intent(in) :: n, ai(0:*), aj(0:*), adiag(0:*)
     14:   PetscScalar, intent(inout) :: x(0:*)

     16:   PetscInt i, jstart, jend
     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:     x(i) = b(i) - sum(aa(jstart:jend)*x(aj(jstart:jend)))
     25:   end do
     26:   !
     27:   ! Backward solve the upper triangular
     28:   !
     29:   do i = n - 1, 0, -1
     30:     jstart = adiag(i) + 1
     31:     jend = ai(i + 1) - 1
     32:     x(i) = x(i) - sum(aa(jstart:jend)*x(aj(jstart:jend)))*aa(adiag(i))
     33:   end do
     34: end subroutine FortranSolveAIJ