nnfc Subroutine

subroutine nnfc(N, R, C, Ic, Ia, Ja, A, Z, B, Lmax, Il, Jl, Ijl, L, D, Umax, Iu, Ju, Iju, U, Row, Tmp, Irl, Jrl, Flag)

NAME

nnfc(3f) - [M_odepack] numerical LDU-factorization of sparse nonsymmetric matrix

DESCRIPTION

numerical ldu-factorization of sparse nonsymmetric matrix and solution of system of linear equations (compressed pointer storage)

       input variables..  n, r, c, ic, ia, ja, a, b,
                          il, jl, ijl, lmax, iu, ju, iju, umax
       output variables.. z, l, d, u, flag

   parameters used internally..

 nia   - irl,  - vectors used to find the rows of  l.  at the kth step
 nia   - jrl       of the factorization,  jrl(k)  points to the head
       -           of a linked list in  jrl  of column indices j
       -           such j .lt. k and  l(k,j)  is nonzero.  zero
       -           indicates the end of the list.  irl(j)  (j.lt.k)
       -           points to the smallest i such that i .ge. k and
       -           l(i,j)  is nonzero.
       -           size of each = n.
 fia   - row   - holds intermediate values in calculation of  u and l.
       -           size = n.
 fia   - tmp   - holds new right-hand side  b*  for solution of the
       -           equation ux = b*.
       -           size = n.

internal variables..

jmin, jmax - indices of the first and last positions in a row to
  be examined.

sum - used in calculating  tmp.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: N
integer, intent(in), dimension(*) :: R
integer, intent(in), dimension(*) :: C
integer, intent(in), dimension(*) :: Ic
integer, intent(in), dimension(*) :: Ia
integer, intent(in), dimension(*) :: Ja
real(kind=dp), intent(in), dimension(*) :: A
real(kind=dp), intent(out), dimension(*) :: Z
real(kind=dp), intent(in), dimension(*) :: B
integer, intent(in) :: Lmax
integer, intent(in), dimension(*) :: Il
integer, intent(in), dimension(*) :: Jl
integer, intent(in), dimension(*) :: Ijl
real(kind=dp), intent(out), dimension(*) :: L
real(kind=dp), intent(out), dimension(*) :: D
integer, intent(in) :: Umax
integer, intent(in), dimension(*) :: Iu
integer, intent(in), dimension(*) :: Ju
integer, intent(in), dimension(*) :: Iju
real(kind=dp), intent(inout), dimension(*) :: U
real(kind=dp), intent(inout), dimension(*) :: Row
real(kind=dp), intent(inout), dimension(*) :: Tmp
integer, intent(inout), dimension(*) :: Irl
integer, intent(inout), dimension(*) :: Jrl
integer, intent(out) :: Flag

Variables

Type Visibility Attributes Name Initial
real(kind=dp), public :: dk
integer, public :: i
integer, public :: i1
integer, public :: i2
integer, public :: ijlb
integer, public :: j
integer, public :: jmax
integer, public :: jmin
integer, public :: k
real(kind=dp), public :: lki
integer, public :: mu
integer, public :: rk
real(kind=dp), public :: sum

Source Code

subroutine nnfc(N,R,C,Ic,Ia,Ja,A,Z,B,Lmax,Il,Jl,Ijl,L,D,Umax,Iu,Ju,Iju,U,Row,Tmp,Irl,Jrl,Flag)
!
!      real  a(*), l(*), d(*), u(*), z(*), b(*), row(*)
!      real tmp(*), lki, sum, dk
integer , intent(in) :: N
integer , intent(in) , dimension(*) :: R
integer , intent(in) , dimension(*) :: C
integer , intent(in) , dimension(*) :: Ic
integer , intent(in) , dimension(*) :: Ia
integer , intent(in) , dimension(*) :: Ja
real(kind=dp) , intent(in) , dimension(*) :: A
real(kind=dp) , intent(out) , dimension(*) :: Z
real(kind=dp) , intent(in) , dimension(*) :: B
integer , intent(in) :: Lmax
integer , intent(in) , dimension(*) :: Il
integer , intent(in) , dimension(*) :: Jl
integer , intent(in) , dimension(*) :: Ijl
real(kind=dp) , intent(out) , dimension(*) :: L
real(kind=dp) , intent(out) , dimension(*) :: D
integer , intent(in) :: Umax
integer , intent(in) , dimension(*) :: Iu
integer , intent(in) , dimension(*) :: Ju
integer , intent(in) , dimension(*) :: Iju
real(kind=dp) , intent(inout) , dimension(*) :: U
real(kind=dp) , intent(inout) , dimension(*) :: Row
real(kind=dp) , intent(inout) , dimension(*) :: Tmp
integer , intent(inout) , dimension(*) :: Irl
integer , intent(inout) , dimension(*) :: Jrl
integer , intent(out) :: Flag
!
real(kind=dp) :: dk , lki , sum
integer :: i , i1 , i2 , ijlb , j , jmax , jmin , k , mu , rk
   !
   !   ******  initialize pointers and test storage  ***********************
   if ( Il(N+1)-1>Lmax ) then
   !
   !  ** error.. insufficient storage for l
      Flag = 4*N + 1
      return
   elseif ( Iu(N+1)-1>Umax ) then
   !  ** error.. insufficient storage for u
      Flag = 7*N + 1
      return
   else
      do k = 1 , N
         Irl(k) = Il(k)
         Jrl(k) = 0
      enddo
   !
   !   ******  for each row  ***********************************************
      do k = 1 , N
   !   ******  reverse jrl and zero row where kth row of l will fill in  ***
         Row(k) = 0
         i1 = 0
         if ( Jrl(k)/=0 ) then
            i = Jrl(k)
            do
               i2 = Jrl(i)
               Jrl(i) = i1
               i1 = i
               Row(i) = 0
               i = i2
               if ( i==0 ) exit
            enddo
         endif
   !   ******  set row to zero where u will fill in  ***********************
         jmin = Iju(k)
         jmax = jmin + Iu(k+1) - Iu(k) - 1
         if ( jmin<=jmax ) then
            do j = jmin , jmax
               Row(Ju(j)) = 0
            enddo
         endif
   !   ******  place kth row of a in row  **********************************
         rk = R(k)
         jmin = Ia(rk)
         jmax = Ia(rk+1) - 1
         do j = jmin , jmax
            Row(Ic(Ja(j))) = A(j)
         enddo
   !   ******  initialize sum, and link through jrl  ***********************
         sum = B(rk)
         i = i1
         if ( i/=0 ) then
            !   ******  assign the kth row of l and adjust row, sum  ****************
            do
               lki = -Row(i)
               !   ******  if l is not required, then comment out the following line  **
               L(Irl(i)) = -lki
               sum = sum + lki*Tmp(i)
               jmin = Iu(i)
               jmax = Iu(i+1) - 1
               if ( jmin<=jmax ) then
                  mu = Iju(i) - jmin
                  do j = jmin , jmax
                     Row(Ju(mu+j)) = Row(Ju(mu+j)) + lki*U(j)
                  enddo
               endif
               i = Jrl(i)
               if ( i==0 ) exit
            enddo
         endif
         !   ******  assign kth row of u and diagonal d, set tmp(k)  *************
         if ( Row(k)==0.0D0 ) then
            !  ** error.. zero pivot
            Flag = 8*N + k
         return
         else
            dk = 1.0D0/Row(k)
            D(k) = dk
            Tmp(k) = sum*dk
            if ( k==N ) cycle
            jmin = Iu(k)
            jmax = Iu(k+1) - 1
            if ( jmin<=jmax ) then
               mu = Iju(k) - jmin
               do j = jmin , jmax
                  U(j) = Row(Ju(mu+j))*dk
               enddo
            endif
            !
            !   ******  update irl and jrl, keeping jrl in decreasing order  ********
            i = i1
            if ( i/=0 )then
               do
                  Irl(i) = Irl(i) + 1
                  i1 = Jrl(i)
                  if ( Irl(i)<Il(i+1) ) then
                     ijlb = Irl(i) - Il(i) + Ijl(i)
                     j = Jl(ijlb)
                     do while ( i<=Jrl(j) )
                        j = Jrl(j)
                     enddo
                     Jrl(i) = Jrl(j)
                     Jrl(j) = i
                  endif
                  i = i1
                  if ( i==0 ) exit
               enddo
         endif
         endif
         if ( Irl(k)<Il(k+1) ) then
            j = Jl(Ijl(k))
            Jrl(k) = Jrl(j)
            Jrl(j) = k
         endif
      enddo
   !
   !   ******  solve  ux = tmp  by back substitution  **********************
      k = N
      do i = 1 , N
         sum = Tmp(k)
         jmin = Iu(k)
         jmax = Iu(k+1) - 1
         if ( jmin<=jmax ) then
            mu = Iju(k) - jmin
            do j = jmin , jmax
               sum = sum - U(j)*Tmp(Ju(mu+j))
            enddo
         endif
         Tmp(k) = sum
         Z(C(k)) = sum
         k = k - 1
      enddo
      Flag = 0
   endif
end subroutine nnfc