Actual source code: fsolve.F

      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 include/finclude/petscdef.h
      8: !
      9:       subroutine FortranSolveAIJ(n,x,ai,aj,adiag,aa,b)
     10:       implicit none
     11:       PetscScalar      x(0:*),aa(0:*),b(0:*)
     12:       PetscInt          n,ai(0:*),aj(0:*),adiag(0:*)

     14:       PetscInt          i,j,jstart,jend
     15:       PetscScalar      sum
     16: !
     17: !     Forward Solve
     18: !
     19:       x(0) = b(0)
     20:       do 20 i=1,n-1
     21:          jstart = ai(i)
     22:          jend   = adiag(i) - 1
     23:          sum    = b(i)
     24:          do 30 j=jstart,jend
     25:             sum  = sum -  aa(j) * x(aj(j))
     26:  30      continue
     27:          x(i) = sum
     28:  20   continue
     29:
     30: !
     31: !     Backward solve the upper triangular
     32: !
     33:       do 40 i=n-1,0,-1
     34:          jstart  = adiag(i) + 1
     35:          jend    = ai(i+1) - 1
     36:          sum     = x(i)
     37:          do 50 j=jstart,jend
     38:             sum = sum - aa(j)* x(aj(j))
     39:  50      continue
     40:          x(i)    = sum * aa(adiag(i))
     41:  40   continue
     42:       return
     43:       end
     44: