nsfc Subroutine

subroutine nsfc(N, R, Ic, Ia, Ja, Jlmax, Il, Jl, Ijl, Jumax, Iu, Ju, Iju, Q, Ira, Jra, Irac, Irl, Jrl, Iru, Jru, Flag)

subroutine nsfc

symbolic ldu-factorization of nonsymmetric sparse matrix

  (compressed pointer storage)


   input variables.. n, r, ic, ia, ja, jlmax, jumax.
   output variables.. il, jl, ijl, iu, ju, iju, flag.

   parameters used internally..

 nia   - q     - suppose  m*  is the result of reordering  m.  if
       -           processing of the ith row of  m*  (hence the ith
       -           row of  u) is being done,  q(j)  is initially
       -           nonzero if  m*(i,j) is nonzero (j.ge.i).  since
       -           values need not be stored, each entry points to the
       -           next nonzero and  q(n+1)  points to the first.  n+1
       -           indicates the end of the list.  for example, if n=9
       -           and the 5th row of  m*  is
       -              0 x x 0 x 0 0 x 0
       -           then  q  will initially be
       -              a a a a 8 a a 10 5           (a - arbitrary).
       -           as the algorithm proceeds, other elements of  q
       -           are inserted in the list because of fillin.
       -           q  is used in an analogous manner to compute the
       -           ith column of  l.
       -           size = n+1.
 nia   - ira,  - vectors used to find the columns of  m.  at the kth
 nia   - jra,      step of the factorization,  irac(k)  points to the
 nia   - irac      head of a linked list in  jra  of row indices i
       -           such that i .ge. k and  m(i,k)  is nonzero.  zero
       -           indicates the end of the list.  ira(i)  (i.ge.k)
       -           points to the smallest j such that j .ge. k and
       -           m(i,j)  is nonzero.
       -           size of each = n.
 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.
 nia   - iru,  - vectors used in a manner analogous to  irl and jrl
 nia   - jru       to find the columns of  u.
       -           size of each = n.

internal variables..

    jlptr - points to the last position used in  jl.
    juptr - points to the last position used in  ju.
    jmin,jmax - are the indices in  a or u  of the first and last
                elements to be examined in a given row.
                for example,  jmin=ia(k), jmax=ia(k+1)-1.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: N
integer, intent(inout) :: R(*)
integer, intent(in) :: Ic(*)
integer, intent(in) :: Ia(*)
integer, intent(in) :: Ja(*)
integer, intent(in) :: Jlmax
integer, intent(inout) :: Il(*)
integer, intent(inout) :: Jl(*)
integer, intent(inout) :: Ijl(*)
integer, intent(in) :: Jumax
integer, intent(inout) :: Iu(*)
integer, intent(inout) :: Ju(*)
integer, intent(inout) :: Iju(*)
integer, intent(inout) :: Q(*)
integer, intent(inout) :: Ira(*)
integer, intent(inout) :: Jra(*)
integer, intent(inout) :: Irac(*)
integer, intent(inout) :: Irl(*)
integer, intent(inout) :: Jrl(*)
integer, intent(inout) :: Iru(*)
integer, intent(inout) :: Jru(*)
integer, intent(out) :: Flag

Variables

Type Visibility Attributes Name Initial
integer, public :: cend
integer, public :: i
integer, public :: i1
integer, public :: iak
integer, public :: irai
integer, public :: irll
integer, public :: irul
integer, public :: j
integer, public :: jaiak
integer, public :: jairai
integer, public :: jlmin
integer, public :: jlptr
integer, public :: jmax
integer, public :: jmin
integer, public :: jtmp
integer, public :: jumin
integer, public :: juptr
integer, public :: k
integer, public :: lasti
integer, public :: lastid
integer, public :: long
integer, public :: luk
integer, public :: m
integer, public :: np1
integer, public :: qm
integer, public :: rend
integer, public :: rk
integer, public :: vj

Source Code

subroutine nsfc(N,R,Ic,Ia,Ja,Jlmax,Il,Jl,Ijl,Jumax,Iu,Ju,Iju,Q,Ira,Jra,Irac,Irl,Jrl,Iru,Jru,Flag)
!
integer,intent(in)    :: N
integer,intent(inout) :: R(*)
integer,intent(in)    :: Ic(*)
integer,intent(in)    :: Ia(*)
integer,intent(in)    :: Ja(*)
integer,intent(in)    :: Jlmax
integer,intent(inout) :: Il(*)
integer,intent(inout) :: Jl(*)
integer,intent(inout) :: Ijl(*)
integer,intent(in)    :: Jumax
integer,intent(inout) :: Iu(*)
integer,intent(inout) :: Ju(*)
integer,intent(inout) :: Iju(*)
integer,intent(inout) :: Q(*)
integer,intent(inout) :: Ira(*)
integer,intent(inout) :: Jra(*)
integer,intent(inout) :: Irac(*)
integer,intent(inout) :: Irl(*)
integer,intent(inout) :: Jrl(*)
integer,intent(inout) :: Iru(*)
integer,intent(inout) :: Jru(*)
integer,intent(out)   :: Flag
!
integer :: cend, i, i1, iak, irai, irll, irul, j, jaiak, jairai, jlmin, jlptr, jmax, jmin, jtmp, jumin, juptr, k, &
         & lasti, lastid, long, luk, m, np1, qm, rend, rk, vj
   !
   !   ******  initialize pointers  ****************************************
   np1 = N + 1
   jlmin = 1
   jlptr = 0
   Il(1) = 1
   jumin = 1
   juptr = 0
   Iu(1) = 1
   do k = 1, N
      Irac(k) = 0
      Jra(k) = 0
      Jrl(k) = 0
      Jru(k) = 0
   enddo
   !   ******  initialize column pointers for a  ***************************
   do k = 1, N
      rk = R(k)
      iak = Ia(rk)
      if ( iak>=Ia(rk+1) ) then
   !
   !  ** error.. null row in a
         Flag = N + rk
         return
      else
         jaiak = Ic(Ja(iak))
         if ( jaiak>k ) then
            !  ** error.. null pivot
            Flag = 5*N + k
            return
         endif
         Jra(k) = Irac(jaiak)
         Irac(jaiak) = k
         Ira(k) = iak
      endif
   enddo
   !
   !   ******  for each column of l and row of u  **************************
   ALL: do k = 1, N
   !
   !   ******  initialize q for computing kth column of l  *****************
      Q(np1) = np1
      luk = -1
      !   ******  by filling in kth column of a  ******************************
      vj = Irac(k)
      if ( vj/=0 ) then
         qm = np1
         do
            m = qm
            qm = Q(m)
            if ( qm>=vj ) then
               if ( qm==vj )then
                  !  ** error.. duplicate entry in a
                  Flag = 2*N + rk
                  return
            endif
               luk = luk + 1
               Q(m) = vj
               Q(vj) = qm
               vj = Jra(vj)
               if ( vj==0 ) exit
               qm = np1
            endif
         enddo
      endif
      !   ******  link through jru  *******************************************
      lastid = 0
      lasti = 0
      Ijl(k) = jlptr
      i = k
      LOOP: do
         i = Jru(i)
         if ( i==0 ) then
         !   ******  lasti is the longest column merged into the kth  ************
         !   ******  see if it equals the entire kth column  *********************
            qm = Q(np1)
            if ( qm/=k ) then
               !  ** error.. null pivot
               Flag = 5*N + k
               return
         endif
            LUKZ: if ( luk/=0 ) then
               if ( lastid/=luk ) then
               !   ******  if not, see if kth column can overlap the previous one  *****
                  KTHCOL: if ( jlmin<=jlptr ) then
                     qm = Q(qm)
                     do j = jlmin, jlptr
                        if ( Jl(j)<qm ) cycle
                        if ( Jl(j)==qm ) then
                           Ijl(k) = j
                           do i = j, jlptr
                              if ( Jl(i)/=qm ) exit KTHCOL
                              qm = Q(qm)
                              if ( qm>N ) exit LUKZ
                           enddo
                           jlptr = j - 1
                        endif
                        exit
                     enddo
                  endif KTHCOL
                  !   ******  move column indices from q to jl, update vectors  ***********
                  jlmin = jlptr + 1
                  Ijl(k) = jlmin
                  if ( luk/=0 ) then
                     jlptr = jlptr + luk
                     if ( jlptr>Jlmax ) then
                        !  ** error.. insufficient storage for jl
                        Flag = 3*N + k
                        return
                     else
                        qm = Q(np1)
                        do j = jlmin, jlptr
                           qm = Q(qm)
                           Jl(j) = qm
                        enddo
                     endif
                  endif
               else
               !   ******  if so, jl can be compressed  ********************************
                  irll = Irl(lasti)
                  Ijl(k) = irll + 1
                  if ( Jl(irll)/=k ) Ijl(k) = Ijl(k) - 1
               endif
            endif LUKZ
            Irl(k) = Ijl(k)
            Il(k+1) = Il(k) + luk
            !
            !   ******  initialize q for computing kth row of u  ********************
            Q(np1) = np1
            luk = -1
            !   ******  by filling in kth row of reordered a  ***********************
            rk = R(k)
            jmin = Ira(k)
            jmax = Ia(rk+1) - 1
            if ( jmin<=jmax ) then
               do j = jmin, jmax
                  vj = Ic(Ja(j))
                  qm = np1
                  do
                     m = qm
                     qm = Q(m)
                     if ( qm>=vj ) then
                        if ( qm==vj )then
                           !  ** error.. duplicate entry in a
                           Flag = 2*N + rk
                           return
                     endif
                        luk = luk + 1
                        Q(m) = vj
                        Q(vj) = qm
                        exit
                     endif
                  enddo
               enddo
            endif
            !   ******  link through jrl,  ******************************************
            lastid = 0
            lasti = 0
            Iju(k) = juptr
            i = k
            i1 = Jrl(k)
            do
               i = i1
               if ( i==0 ) then
               !   ******  update jrl(k) and irl(k)  ***********************************
                  if ( Il(k+1)>Il(k) ) then
                     j = Jl(Irl(k))
                     Jrl(k) = Jrl(j)
                     Jrl(j) = k
                  endif
                  !   ******  lasti is the longest row merged into the kth  ***************
                  !   ******  see if it equals the entire kth row  ************************
                  qm = Q(np1)
                  if ( qm/=k )then
                     !  ** error.. null pivot
                     Flag = 5*N + k
                     return
               endif
                  LUKY: if ( luk/=0 ) then
                     if ( lastid/=luk ) then
                     !   ******  if not, see if kth row can overlap the previous one  ********
                        KTHROW: if ( jumin<=juptr ) then
                           qm = Q(qm)
                           do j = jumin, juptr
                              if ( Ju(j)<qm ) cycle
                              if ( Ju(j)==qm ) then
                                 Iju(k) = j
                                 do i = j, juptr
                                    if ( Ju(i)/=qm ) exit KTHROW
                                    qm = Q(qm)
                                    if ( qm>N ) exit LUKY
                                 enddo
                                 juptr = j - 1
                              endif
                              exit
                           enddo
                        endif KTHROW
                        !   ******  move row indices from q to ju, update vectors  **************
                        jumin = juptr + 1
                        Iju(k) = jumin
                        if ( luk/=0 ) then
                           juptr = juptr + luk
                           if ( juptr>Jumax ) then
                           !  ** error.. insufficient storage for ju
                              Flag = 6*N + k
                              return
                           else
                              qm = Q(np1)
                              do j = jumin, juptr
                                 qm = Q(qm)
                                 Ju(j) = qm
                              enddo
                           endif
                        endif
                     else
                        !   ******  if so, ju can be compressed  ********************************
                        irul = Iru(lasti)
                        Iju(k) = irul + 1
                        if ( Ju(irul)/=k ) Iju(k) = Iju(k) - 1
                     endif
                  endif LUKY
                  Iru(k) = Iju(k)
                  Iu(k+1) = Iu(k) + luk
               !
                  !   ******  update iru, jru  ********************************************
                  i = k
                  exit LOOP
               else
                  i1 = Jrl(i)
                  qm = np1
                  jmin = Iru(i)
                  jmax = Iju(i) + Iu(i+1) - Iu(i) - 1
                  long = jmax - jmin
                  if ( long>=0 ) then
                     jtmp = Ju(jmin)
                     if ( jtmp/=k ) then
                     !   ******  update irl and jrl, *****************************************
                        long = long + 1
                        cend = Ijl(i) + Il(i+1) - Il(i)
                        Irl(i) = Irl(i) + 1
                        if ( Irl(i)<cend ) then
                           j = Jl(Irl(i))
                           Jrl(i) = Jrl(j)
                           Jrl(j) = i
                        endif
                     endif
                     if ( lastid<long ) then
                        lasti = i
                        lastid = long
                     endif
                     !   ******  and merge the corresponding rows into the kth row  **********
                     do j = jmin, jmax
                        vj = Ju(j)
                        do
                           m = qm
                           qm = Q(m)
                           if ( qm>=vj ) then
                              if ( qm/=vj ) then
                                 luk = luk + 1
                                 Q(m) = vj
                                 Q(vj) = qm
                                 qm = vj
                              endif
                              exit
                           endif
                        enddo
                     enddo
                  endif
               endif
            enddo
         else
            qm = np1
            jmin = Irl(i)
            jmax = Ijl(i) + Il(i+1) - Il(i) - 1
            long = jmax - jmin
            if ( long>=0 ) then
               jtmp = Jl(jmin)
               if ( jtmp/=k ) long = long + 1
               if ( jtmp==k ) R(i) = -R(i)
               if ( lastid<long ) then
                  lasti = i
                  lastid = long
               endif
               !   ******  and merge the corresponding columns into the kth column  ****
               do j = jmin, jmax
                  vj = Jl(j)
                  do
                     m = qm
                     qm = Q(m)
                     if ( qm>=vj ) then
                        if ( qm/=vj ) then
                           luk = luk + 1
                           Q(m) = vj
                           Q(vj) = qm
                           qm = vj
                        endif
                        exit
                     endif
                  enddo
               enddo
            endif
         endif
      enddo LOOP
      BACK: do
         i1 = Jru(i)
         if ( R(i)<0 ) then
            R(i) = -R(i)
         else
            rend = Iju(i) + Iu(i+1) - Iu(i)
            if ( Iru(i)<rend ) then
               j = Ju(Iru(i))
               Jru(i) = Jru(j)
               Jru(j) = i
            endif
         endif
         i = i1
         if ( i==0 ) then
         !
         !   ******  update ira, jra, irac  **************************************
            i = Irac(k)
            if ( i==0 ) cycle ALL
         else
            Iru(i) = Iru(i) + 1
         cycle BACK
         endif
         exit BACK
      enddo BACK
      INFINITE: do
         i1 = Jra(i)
         Ira(i) = Ira(i) + 1
         if ( Ira(i)<Ia(R(i)+1) ) then
            irai = Ira(i)
            jairai = Ic(Ja(irai))
            if ( jairai<=i ) then
               Jra(i) = Irac(jairai)
               Irac(jairai) = i
            endif
         endif
         i = i1
         if ( i==0 ) exit INFINITE
      enddo INFINITE
   enddo ALL
   !
   Ijl(N) = jlptr
   Iju(N) = juptr
   Flag = 0
end subroutine nsfc