*DECK CDRV
subroutine cdrv
* (n, r,c,ic, ia,ja,a, b, z, nsp,isp,rsp,esp, path, flag)
c*** subroutine cdrv
c*** driver for subroutines for solving sparse nonsymmetric systems of
c linear equations (compressed pointer storage)
c
c
c parameters
c class abbreviations are--
c n - integer variable
c f - real variable
c v - supplies a value to the driver
c r - returns a result from the driver
c i - used internally by the driver
c a - array
c
c class - parameter
c ------+----------
c -
c the nonzero entries of the coefficient matrix m are stored
c row-by-row in the array a. to identify the individual nonzero
c entries in each row, we need to know in which column each entry
c lies. the column indices which correspond to the nonzero entries
c of m are stored in the array ja. i.e., if a(k) = m(i,j), then
c ja(k) = j. in addition, we need to know where each row starts and
c how long it is. the index positions in ja and a where the rows of
c m begin are stored in the array ia. i.e., if m(i,j) is the first
c nonzero entry (stored) in the i-th row and a(k) = m(i,j), then
c ia(i) = k. moreover, the index in ja and a of the first location
c following the last element in the last row is stored in ia(n+1).
c thus, the number of entries in the i-th row is given by
c ia(i+1) - ia(i), the nonzero entries of the i-th row are stored
c consecutively in
c a(ia(i)), a(ia(i)+1), ..., a(ia(i+1)-1),
c and the corresponding column indices are stored consecutively in
c ja(ia(i)), ja(ia(i)+1), ..., ja(ia(i+1)-1).
c for example, the 5 by 5 matrix
c ( 1. 0. 2. 0. 0.)
c ( 0. 3. 0. 0. 0.)
c m = ( 0. 4. 5. 6. 0.)
c ( 0. 0. 0. 7. 0.)
c ( 0. 0. 0. 8. 9.)
c would be stored as
c - 1 2 3 4 5 6 7 8 9
c ---+--------------------------
c ia - 1 3 4 7 8 10
c ja - 1 3 2 2 3 4 4 4 5
c a - 1. 2. 3. 4. 5. 6. 7. 8. 9. .
c
c nv - n - number of variables/equations.
c fva - a - nonzero entries of the coefficient matrix m, stored
c - by rows.
c - size = number of nonzero entries in m.
c nva - ia - pointers to delimit the rows in a.
c - size = n+1.
c nva - ja - column numbers corresponding to the elements of a.
c - size = size of a.
c fva - b - right-hand side b. b and z can the same array.
c - size = n.
c fra - z - solution x. b and z can be the same array.
c - size = n.
c
c the rows and columns of the original matrix m can be
c reordered (e.g., to reduce fillin or ensure numerical stability)
c before calling the driver. if no reordering is done, then set
c r(i) = c(i) = ic(i) = i for i=1,...,n. the solution z is returned
c in the original order.
c if the columns have been reordered (i.e., c(i).ne.i for some
c i), then the driver will call a subroutine (nroc) which rearranges
c each row of ja and a, leaving the rows in the original order, but
c placing the elements of each row in increasing order with respect
c to the new ordering. if path.ne.1, then nroc is assumed to have
c been called already.
c
c nva - r - ordering of the rows of m.
c - size = n.
c nva - c - ordering of the columns of m.
c - size = n.
c nva - ic - inverse of the ordering of the columns of m. i.e.,
c - ic(c(i)) = i for i=1,...,n.
c - size = n.
c
c the solution of the system of linear equations is divided into
c three stages --
c nsfc -- the matrix m is processed symbolically to determine where
c fillin will occur during the numeric factorization.
c nnfc -- the matrix m is factored numerically into the product ldu
c of a unit lower triangular matrix l, a diagonal matrix
c d, and a unit upper triangular matrix u, and the system
c mx = b is solved.
c nnsc -- the linear system mx = b is solved using the ldu
c or factorization from nnfc.
c nntc -- the transposed linear system mt x = b is solved using
c the ldu factorization from nnf.
c for several systems whose coefficient matrices have the same
c nonzero structure, nsfc need be done only once (for the first
c system). then nnfc is done once for each additional system. for
c several systems with the same coefficient matrix, nsfc and nnfc
c need be done only once (for the first system). then nnsc or nntc
c is done once for each additional right-hand side.
c
c nv - path - path specification. values and their meanings are --
c - 1 perform nroc, nsfc, and nnfc.
c - 2 perform nnfc only (nsfc is assumed to have been
c - done in a manner compatible with the storage
c - allocation used in the driver).
c - 3 perform nnsc only (nsfc and nnfc are assumed to
c - have been done in a manner compatible with the
c - storage allocation used in the driver).
c - 4 perform nntc only (nsfc and nnfc are assumed to
c - have been done in a manner compatible with the
c - storage allocation used in the driver).
c - 5 perform nroc and nsfc.
c
c various errors are detected by the driver and the individual
c subroutines.
c
c nr - flag - error flag. values and their meanings are --
c - 0 no errors detected
c - n+k null row in a -- row = k
c - 2n+k duplicate entry in a -- row = k
c - 3n+k insufficient storage in nsfc -- row = k
c - 4n+1 insufficient storage in nnfc
c - 5n+k null pivot -- row = k
c - 6n+k insufficient storage in nsfc -- row = k
c - 7n+1 insufficient storage in nnfc
c - 8n+k zero pivot -- row = k
c - 10n+1 insufficient storage in cdrv
c - 11n+1 illegal path specification
c
c working storage is needed for the factored form of the matrix
c m plus various temporary vectors. the arrays isp and rsp should be
c equivalenced. integer storage is allocated from the beginning of
c isp and real storage from the end of rsp.
c
c nv - nsp - declared dimension of rsp. nsp generally must
c - be larger than 8n+2 + 2k (where k = (number of
c - nonzero entries in m)).
c nvira - isp - integer working storage divided up into various arrays
c - needed by the subroutines. isp and rsp should be
c - equivalenced.
c - size = lratio*nsp.
c fvira - rsp - real working storage divided up into various arrays
c - needed by the subroutines. isp and rsp should be
c - equivalenced.
c - size = nsp.
c nr - esp - if sufficient storage was available to perform the
c - symbolic factorization (nsfc), then esp is set to
c - the amount of excess storage provided (negative if
c - insufficient storage was available to perform the
c - numeric factorization (nnfc)).
c
c
c conversion to double precision
c
c to convert these routines for double precision arrays..
c (1) use the double precision declarations in place of the real
c declarations in each subprogram, as given in comment cards.
c (2) change the data-loaded value of the integer lratio
c in subroutine cdrv, as indicated below.
c (3) change e0 to d0 in the constants in statement number 10
c in subroutine nnfc and the line following that.
c
integer r(*), c(*), ic(*), ia(*), ja(*), isp(*), esp, path,
* flag, d, u, q, row, tmp, ar, umax
c real a(*), b(*), z(*), rsp(*)
double precision a(*), b(*), z(*), rsp(*)
c
c set lratio equal to the ratio between the length of floating point
c and integer array data. e. g., lratio = 1 for (real, integer),
c lratio = 2 for (double precision, integer)
c
data lratio/2/
c
if (path.lt.1 .or. 5.lt.path) go to 111
c******initialize and divide up temporary storage *******************
il = 1
ijl = il + (n+1)
iu = ijl + n
iju = iu + (n+1)
irl = iju + n
jrl = irl + n
jl = jrl + n
c
c ****** reorder a if necessary, call nsfc if flag is set ***********
if ((path-1) * (path-5) .ne. 0) go to 5
max = (lratio*nsp + 1 - jl) - (n+1) - 5*n
jlmax = max/2
q = jl + jlmax
ira = q + (n+1)
jra = ira + n
irac = jra + n
iru = irac + n
jru = iru + n
jutmp = jru + n
jumax = lratio*nsp + 1 - jutmp
esp = max/lratio
if (jlmax.le.0 .or. jumax.le.0) go to 110
c
do 1 i=1,n
if (c(i).ne.i) go to 2
1 continue
go to 3
2 ar = nsp + 1 - n
call nroc
* (n, ic, ia,ja,a, isp(il), rsp(ar), isp(iu), flag)
if (flag.ne.0) go to 100
c
3 call nsfc
* (n, r, ic, ia,ja,
* jlmax, isp(il), isp(jl), isp(ijl),
* jumax, isp(iu), isp(jutmp), isp(iju),
* isp(q), isp(ira), isp(jra), isp(irac),
* isp(irl), isp(jrl), isp(iru), isp(jru), flag)
if(flag .ne. 0) go to 100
c ****** move ju next to jl *****************************************
jlmax = isp(ijl+n-1)
ju = jl + jlmax
jumax = isp(iju+n-1)
if (jumax.le.0) go to 5
do 4 j=1,jumax
4 isp(ju+j-1) = isp(jutmp+j-1)
c
c ****** call remaining subroutines *********************************
5 jlmax = isp(ijl+n-1)
ju = jl + jlmax
jumax = isp(iju+n-1)
l = (ju + jumax - 2 + lratio) / lratio + 1
lmax = isp(il+n) - 1
d = l + lmax
u = d + n
row = nsp + 1 - n
tmp = row - n
umax = tmp - u
esp = umax - (isp(iu+n) - 1)
c
if ((path-1) * (path-2) .ne. 0) go to 6
if (umax.lt.0) go to 110
call nnfc
* (n, r, c, ic, ia, ja, a, z, b,
* lmax, isp(il), isp(jl), isp(ijl), rsp(l), rsp(d),
* umax, isp(iu), isp(ju), isp(iju), rsp(u),
* rsp(row), rsp(tmp), isp(irl), isp(jrl), flag)
if(flag .ne. 0) go to 100
c
6 if ((path-3) .ne. 0) go to 7
call nnsc
* (n, r, c, isp(il), isp(jl), isp(ijl), rsp(l),
* rsp(d), isp(iu), isp(ju), isp(iju), rsp(u),
* z, b, rsp(tmp))
c
7 if ((path-4) .ne. 0) go to 8
call nntc
* (n, r, c, isp(il), isp(jl), isp(ijl), rsp(l),
* rsp(d), isp(iu), isp(ju), isp(iju), rsp(u),
* z, b, rsp(tmp))
8 return
c
c ** error.. error detected in nroc, nsfc, nnfc, or nnsc
100 return
c ** error.. insufficient storage
110 flag = 10*n + 1
return
c ** error.. illegal path specification
111 flag = 11*n + 1
return
end