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: