share/odepack/fortran/nsfc.f

   1       subroutine nsfc
   2      *      (n, r, ic, ia,ja, jlmax,il,jl,ijl, jumax,iu,ju,iju,
   3      *       q, ira,jra, irac, irl,jrl, iru,jru, flag)
   4 c*** subroutine nsfc
   5 c*** symbolic ldu-factorization of nonsymmetric sparse matrix
   6 c      (compressed pointer storage)
   7 c
   8 c
   9 c       input variables.. n, r, ic, ia, ja, jlmax, jumax.
  10 c       output variables.. il, jl, ijl, iu, ju, iju, flag.
  11 c
  12 c       parameters used internally..
  13 c nia   - q     - suppose  m*  is the result of reordering  m.  if
  14 c       -           processing of the ith row of  m*  (hence the ith
  15 c       -           row of  u) is being done,  q(j)  is initially
  16 c       -           nonzero if  m*(i,j) is nonzero (j.ge.i).  since
  17 c       -           values need not be stored, each entry points to the
  18 c       -           next nonzero and  q(n+1)  points to the first.  n+1
  19 c       -           indicates the end of the list.  for example, if n=9
  20 c       -           and the 5th row of  m*  is
  21 c       -              0 x x 0 x 0 0 x 0
  22 c       -           then  q  will initially be
  23 c       -              a a a a 8 a a 10 5           (a - arbitrary).
  24 c       -           as the algorithm proceeds, other elements of  q
  25 c       -           are inserted in the list because of fillin.
  26 c       -           q  is used in an analogous manner to compute the
  27 c       -           ith column of  l.
  28 c       -           size = n+1.
  29 c nia   - ira,  - vectors used to find the columns of  m.  at the kth
  30 c nia   - jra,      step of the factorization,  irac(k)  points to the
  31 c nia   - irac      head of a linked list in  jra  of row indices i
  32 c       -           such that i .ge. k and  m(i,k)  is nonzero.  zero
  33 c       -           indicates the end of the list.  ira(i)  (i.ge.k)
  34 c       -           points to the smallest j such that j .ge. k and
  35 c       -           m(i,j)  is nonzero.
  36 c       -           size of each = n.
  37 c nia   - irl,  - vectors used to find the rows of  l.  at the kth step
  38 c nia   - jrl       of the factorization,  jrl(k)  points to the head
  39 c       -           of a linked list in  jrl  of column indices j
  40 c       -           such j .lt. k and  l(k,j)  is nonzero.  zero
  41 c       -           indicates the end of the list.  irl(j)  (j.lt.k)
  42 c       -           points to the smallest i such that i .ge. k and
  43 c       -           l(i,j)  is nonzero.
  44 c       -           size of each = n.
  45 c nia   - iru,  - vectors used in a manner analogous to  irl and jrl
  46 c nia   - jru       to find the columns of  u.
  47 c       -           size of each = n.
  48 c
  49 c  internal variables..
  50 c    jlptr - points to the last position used in  jl.
  51 c    juptr - points to the last position used in  ju.
  52 c    jmin,jmax - are the indices in  a or u  of the first and last
  53 c                elements to be examined in a given row.
  54 c                for example,  jmin=ia(k), jmax=ia(k+1)-1.
  55 c
  56       integer cend, qm, rend, rk, vj
  57       integer ia(*), ja(*), ira(*), jra(*), il(*), jl(*), ijl(*)
  58       integer iu(*), ju(*), iju(*), irl(*), jrl(*), iru(*), jru(*)
  59       integer r(*), ic(*), q(*), irac(*), flag
  60 c
  61 c  ******  initialize pointers  ****************************************
  62       np1 = n + 1
  63       jlmin = 1
  64       jlptr = 0
  65       il(1) = 1
  66       jumin = 1
  67       juptr = 0
  68       iu(1) = 1
  69       do 1 k=1,n
  70         irac(k) = 0
  71         jra(k) = 0
  72         jrl(k) = 0
  73    1    jru(k) = 0
  74 c  ******  initialize column pointers for a  ***************************
  75       do 2 k=1,n
  76         rk = r(k)
  77         iak = ia(rk)
  78         if (iak .ge. ia(rk+1))  go to 101
  79         jaiak = ic(ja(iak))
  80         if (jaiak .gt. k)  go to 105
  81         jra(k) = irac(jaiak)
  82         irac(jaiak) = k
  83    2    ira(k) = iak
  84 c
  85 c  ******  for each column of l and row of u  **************************
  86       do 41 k=1,n
  87 c
  88 c  ******  initialize q for computing kth column of l  *****************
  89         q(np1) = np1
  90         luk = -1
  91 c  ******  by filling in kth column of a  ******************************
  92         vj = irac(k)
  93         if (vj .eq. 0)  go to 5
  94    3      qm = np1
  95    4      m = qm
  96           qm =  q(m)
  97           if (qm .lt. vj)  go to 4
  98           if (qm .eq. vj)  go to 102
  99             luk = luk + 1
 100             q(m) = vj
 101             q(vj) = qm
 102             vj = jra(vj)
 103             if (vj .ne. 0)  go to 3
 104 c  ******  link through jru  *******************************************
 105    5    lastid = 0
 106         lasti = 0
 107         ijl(k) = jlptr
 108         i = k
 109    6      i = jru(i)
 110           if (i .eq. 0)  go to 10
 111           qm = np1
 112           jmin = irl(i)
 113           jmax = ijl(i) + il(i+1) - il(i) - 1
 114           long = jmax - jmin
 115           if (long .lt. 0)  go to 6
 116           jtmp = jl(jmin)
 117           if (jtmp .ne. k)  long = long + 1
 118           if (jtmp .eq. k)  r(i) = -r(i)
 119           if (lastid .ge. long)  go to 7
 120             lasti = i
 121             lastid = long
 122 c  ******  and merge the corresponding columns into the kth column  ****
 123    7      do 9 j=jmin,jmax
 124             vj = jl(j)
 125    8        m = qm
 126             qm = q(m)
 127             if (qm .lt. vj)  go to 8
 128             if (qm .eq. vj)  go to 9
 129               luk = luk + 1
 130               q(m) = vj
 131               q(vj) = qm
 132               qm = vj
 133    9        continue
 134             go to 6
 135 c  ******  lasti is the longest column merged into the kth  ************
 136 c  ******  see if it equals the entire kth column  *********************
 137   10    qm = q(np1)
 138         if (qm .ne. k)  go to 105
 139         if (luk .eq. 0)  go to 17
 140         if (lastid .ne. luk)  go to 11
 141 c  ******  if so, jl can be compressed  ********************************
 142         irll = irl(lasti)
 143         ijl(k) = irll + 1
 144         if (jl(irll) .ne. k)  ijl(k) = ijl(k) - 1
 145         go to 17
 146 c  ******  if not, see if kth column can overlap the previous one  *****
 147   11    if (jlmin .gt. jlptr)  go to 15
 148         qm = q(qm)
 149         do 12 j=jlmin,jlptr
 150           if (jl(j) - qm)  12, 13, 15
 151   12      continue
 152         go to 15
 153   13    ijl(k) = j
 154         do 14 i=j,jlptr
 155           if (jl(i) .ne. qm)  go to 15
 156           qm = q(qm)
 157           if (qm .gt. n)  go to 17
 158   14      continue
 159         jlptr = j - 1
 160 c  ******  move column indices from q to jl, update vectors  ***********
 161   15    jlmin = jlptr + 1
 162         ijl(k) = jlmin
 163         if (luk .eq. 0)  go to 17
 164         jlptr = jlptr + luk
 165         if (jlptr .gt. jlmax)  go to 103
 166           qm = q(np1)
 167           do 16 j=jlmin,jlptr
 168             qm = q(qm)
 169   16        jl(j) = qm
 170   17    irl(k) = ijl(k)
 171         il(k+1) = il(k) + luk
 172 c
 173 c  ******  initialize q for computing kth row of u  ********************
 174         q(np1) = np1
 175         luk = -1
 176 c  ******  by filling in kth row of reordered a  ***********************
 177         rk = r(k)
 178         jmin = ira(k)
 179         jmax = ia(rk+1) - 1
 180         if (jmin .gt. jmax)  go to 20
 181         do 19 j=jmin,jmax
 182           vj = ic(ja(j))
 183           qm = np1
 184   18      m = qm
 185           qm = q(m)
 186           if (qm .lt. vj)  go to 18
 187           if (qm .eq. vj)  go to 102
 188             luk = luk + 1
 189             q(m) = vj
 190             q(vj) = qm
 191   19      continue
 192 c  ******  link through jrl,  ******************************************
 193   20    lastid = 0
 194         lasti = 0
 195         iju(k) = juptr
 196         i = k
 197         i1 = jrl(k)
 198   21      i = i1
 199           if (i .eq. 0)  go to 26
 200           i1 = jrl(i)
 201           qm = np1
 202           jmin = iru(i)
 203           jmax = iju(i) + iu(i+1) - iu(i) - 1
 204           long = jmax - jmin
 205           if (long .lt. 0)  go to 21
 206           jtmp = ju(jmin)
 207           if (jtmp .eq. k)  go to 22
 208 c  ******  update irl and jrl, *****************************************
 209             long = long + 1
 210             cend = ijl(i) + il(i+1) - il(i)
 211             irl(i) = irl(i) + 1
 212             if (irl(i) .ge. cend)  go to 22
 213               j = jl(irl(i))
 214               jrl(i) = jrl(j)
 215               jrl(j) = i
 216   22      if (lastid .ge. long)  go to 23
 217             lasti = i
 218             lastid = long
 219 c  ******  and merge the corresponding rows into the kth row  **********
 220   23      do 25 j=jmin,jmax
 221             vj = ju(j)
 222   24        m = qm
 223             qm = q(m)
 224             if (qm .lt. vj)  go to 24
 225             if (qm .eq. vj)  go to 25
 226               luk = luk + 1
 227               q(m) = vj
 228               q(vj) = qm
 229               qm = vj
 230   25        continue
 231           go to 21
 232 c  ******  update jrl(k) and irl(k)  ***********************************
 233   26    if (il(k+1) .le. il(k))  go to 27
 234           j = jl(irl(k))
 235           jrl(k) = jrl(j)
 236           jrl(j) = k
 237 c  ******  lasti is the longest row merged into the kth  ***************
 238 c  ******  see if it equals the entire kth row  ************************
 239   27    qm = q(np1)
 240         if (qm .ne. k)  go to 105
 241         if (luk .eq. 0)  go to 34
 242         if (lastid .ne. luk)  go to 28
 243 c  ******  if so, ju can be compressed  ********************************
 244         irul = iru(lasti)
 245         iju(k) = irul + 1
 246         if (ju(irul) .ne. k)  iju(k) = iju(k) - 1
 247         go to 34
 248 c  ******  if not, see if kth row can overlap the previous one  ********
 249   28    if (jumin .gt. juptr)  go to 32
 250         qm = q(qm)
 251         do 29 j=jumin,juptr
 252           if (ju(j) - qm)  29, 30, 32
 253   29      continue
 254         go to 32
 255   30    iju(k) = j
 256         do 31 i=j,juptr
 257           if (ju(i) .ne. qm)  go to 32
 258           qm = q(qm)
 259           if (qm .gt. n)  go to 34
 260   31      continue
 261         juptr = j - 1
 262 c  ******  move row indices from q to ju, update vectors  **************
 263   32    jumin = juptr + 1
 264         iju(k) = jumin
 265         if (luk .eq. 0)  go to 34
 266         juptr = juptr + luk
 267         if (juptr .gt. jumax)  go to 106
 268           qm = q(np1)
 269           do 33 j=jumin,juptr
 270             qm = q(qm)
 271   33        ju(j) = qm
 272   34    iru(k) = iju(k)
 273         iu(k+1) = iu(k) + luk
 274 c
 275 c  ******  update iru, jru  ********************************************
 276         i = k
 277   35      i1 = jru(i)
 278           if (r(i) .lt. 0)  go to 36
 279           rend = iju(i) + iu(i+1) - iu(i)
 280           if (iru(i) .ge. rend)  go to 37
 281             j = ju(iru(i))
 282             jru(i) = jru(j)
 283             jru(j) = i
 284             go to 37
 285   36      r(i) = -r(i)
 286   37      i = i1
 287           if (i .eq. 0)  go to 38
 288           iru(i) = iru(i) + 1
 289           go to 35
 290 c
 291 c  ******  update ira, jra, irac  **************************************
 292   38    i = irac(k)
 293         if (i .eq. 0)  go to 41
 294   39      i1 = jra(i)
 295           ira(i) = ira(i) + 1
 296           if (ira(i) .ge. ia(r(i)+1))  go to 40
 297           irai = ira(i)
 298           jairai = ic(ja(irai))
 299           if (jairai .gt. i)  go to 40
 300           jra(i) = irac(jairai)
 301           irac(jairai) = i
 302   40      i = i1
 303           if (i .ne. 0)  go to 39
 304   41    continue
 305 c
 306       ijl(n) = jlptr
 307       iju(n) = juptr
 308       flag = 0
 309       return
 310 c
 311 c ** error.. null row in a
 312  101  flag = n + rk
 313       return
 314 c ** error.. duplicate entry in a
 315  102  flag = 2*n + rk
 316       return
 317 c ** error.. insufficient storage for jl
 318  103  flag = 3*n + k
 319       return
 320 c ** error.. null pivot
 321  105  flag = 5*n + k
 322       return
 323 c ** error.. insufficient storage for ju
 324  106  flag = 6*n + k
 325       return
 326       end