share/odepack/fortran/nnfc.f

   1       subroutine nnfc
   2      *     (n, r,c,ic, ia,ja,a, z, b,
   3      *      lmax,il,jl,ijl,l, d, umax,iu,ju,iju,u,
   4      *      row, tmp, irl,jrl, flag)
   5 c*** subroutine nnfc
   6 c*** numerical ldu-factorization of sparse nonsymmetric matrix and
   7 c      solution of system of linear equations (compressed pointer
   8 c      storage)
   9 c
  10 c
  11 c       input variables..  n, r, c, ic, ia, ja, a, b,
  12 c                          il, jl, ijl, lmax, iu, ju, iju, umax
  13 c       output variables.. z, l, d, u, flag
  14 c
  15 c       parameters used internally..
  16 c nia   - irl,  - vectors used to find the rows of  l.  at the kth step
  17 c nia   - jrl       of the factorization,  jrl(k)  points to the head
  18 c       -           of a linked list in  jrl  of column indices j
  19 c       -           such j .lt. k and  l(k,j)  is nonzero.  zero
  20 c       -           indicates the end of the list.  irl(j)  (j.lt.k)
  21 c       -           points to the smallest i such that i .ge. k and
  22 c       -           l(i,j)  is nonzero.
  23 c       -           size of each = n.
  24 c fia   - row   - holds intermediate values in calculation of  u and l.
  25 c       -           size = n.
  26 c fia   - tmp   - holds new right-hand side  b*  for solution of the
  27 c       -           equation ux = b*.
  28 c       -           size = n.
  29 c
  30 c  internal variables..
  31 c    jmin, jmax - indices of the first and last positions in a row to
  32 c      be examined.
  33 c    sum - used in calculating  tmp.
  34 c
  35       integer rk,umax
  36       integer  r(*), c(*), ic(*), ia(*), ja(*), il(*), jl(*), ijl(*)
  37       integer  iu(*), ju(*), iju(*), irl(*), jrl(*), flag
  38 c     real  a(*), l(*), d(*), u(*), z(*), b(*), row(*)
  39 c     real tmp(*), lki, sum, dk
  40       double precision  a(*), l(*), d(*), u(*), z(*), b(*), row(*)
  41       double precision  tmp(*), lki, sum, dk
  42 c
  43 c  ******  initialize pointers and test storage  ***********************
  44       if(il(n+1)-1 .gt. lmax) go to 104
  45       if(iu(n+1)-1 .gt. umax) go to 107
  46       do 1 k=1,n
  47         irl(k) = il(k)
  48         jrl(k) = 0
  49    1    continue
  50 c
  51 c  ******  for each row  ***********************************************
  52       do 19 k=1,n
  53 c  ******  reverse jrl and zero row where kth row of l will fill in  ***
  54         row(k) = 0
  55         i1 = 0
  56         if (jrl(k) .eq. 0) go to 3
  57         i = jrl(k)
  58    2    i2 = jrl(i)
  59         jrl(i) = i1
  60         i1 = i
  61         row(i) = 0
  62         i = i2
  63         if (i .ne. 0) go to 2
  64 c  ******  set row to zero where u will fill in  ***********************
  65    3    jmin = iju(k)
  66         jmax = jmin + iu(k+1) - iu(k) - 1
  67         if (jmin .gt. jmax) go to 5
  68         do 4 j=jmin,jmax
  69    4      row(ju(j)) = 0
  70 c  ******  place kth row of a in row  **********************************
  71    5    rk = r(k)
  72         jmin = ia(rk)
  73         jmax = ia(rk+1) - 1
  74         do 6 j=jmin,jmax
  75           row(ic(ja(j))) = a(j)
  76    6      continue
  77 c  ******  initialize sum, and link through jrl  ***********************
  78         sum = b(rk)
  79         i = i1
  80         if (i .eq. 0) go to 10
  81 c  ******  assign the kth row of l and adjust row, sum  ****************
  82    7      lki = -row(i)
  83 c  ******  if l is not required, then comment out the following line  **
  84           l(irl(i)) = -lki
  85           sum = sum + lki * tmp(i)
  86           jmin = iu(i)
  87           jmax = iu(i+1) - 1
  88           if (jmin .gt. jmax) go to 9
  89           mu = iju(i) - jmin
  90           do 8 j=jmin,jmax
  91    8        row(ju(mu+j)) = row(ju(mu+j)) + lki * u(j)
  92    9      i = jrl(i)
  93           if (i .ne. 0) go to 7
  94 c
  95 c  ******  assign kth row of u and diagonal d, set tmp(k)  *************
  96   10    if (row(k) .eq. 0.0d0) go to 108
  97         dk = 1.0d0 / row(k)
  98         d(k) = dk
  99         tmp(k) = sum * dk
 100         if (k .eq. n) go to 19
 101         jmin = iu(k)
 102         jmax = iu(k+1) - 1
 103         if (jmin .gt. jmax)  go to 12
 104         mu = iju(k) - jmin
 105         do 11 j=jmin,jmax
 106   11      u(j) = row(ju(mu+j)) * dk
 107   12    continue
 108 c
 109 c  ******  update irl and jrl, keeping jrl in decreasing order  ********
 110         i = i1
 111         if (i .eq. 0) go to 18
 112   14    irl(i) = irl(i) + 1
 113         i1 = jrl(i)
 114         if (irl(i) .ge. il(i+1)) go to 17
 115         ijlb = irl(i) - il(i) + ijl(i)
 116         j = jl(ijlb)
 117   15    if (i .gt. jrl(j)) go to 16
 118           j = jrl(j)
 119           go to 15
 120   16    jrl(i) = jrl(j)
 121         jrl(j) = i
 122   17    i = i1
 123         if (i .ne. 0) go to 14
 124   18    if (irl(k) .ge. il(k+1)) go to 19
 125         j = jl(ijl(k))
 126         jrl(k) = jrl(j)
 127         jrl(j) = k
 128   19    continue
 129 c
 130 c  ******  solve  ux = tmp  by back substitution  **********************
 131       k = n
 132       do 22 i=1,n
 133         sum =  tmp(k)
 134         jmin = iu(k)
 135         jmax = iu(k+1) - 1
 136         if (jmin .gt. jmax)  go to 21
 137         mu = iju(k) - jmin
 138         do 20 j=jmin,jmax
 139   20      sum = sum - u(j) * tmp(ju(mu+j))
 140   21    tmp(k) =  sum
 141         z(c(k)) =  sum
 142   22    k = k-1
 143       flag = 0
 144       return
 145 c
 146 c ** error.. insufficient storage for l
 147  104  flag = 4*n + 1
 148       return
 149 c ** error.. insufficient storage for u
 150  107  flag = 7*n + 1
 151       return
 152 c ** error.. zero pivot
 153  108  flag = 8*n + k
 154       return
 155       end