cdrv Subroutine

subroutine cdrv(N, R, C, Ic, Ia, Ja, A, B, Z, Nsp, Isp, Rsp, Esp, Path, Flag)

NAME

cdrv(3f) - [M_odepack ] driver for solving sparse non-symmetric systems of linear equations

DESCRIPTION

driver for subroutines for solving sparse nonsymmetric systems of linear equations (compressed pointer storage)

PARAMETERS

class abbreviations are--
   n - integer variable
   f - real variable
   v - supplies a value to the driver
   r - returns a result from the driver
   i - used internally by the driver
   a - array

class - parameter

the nonzero entries of the coefficient matrix m are stored row-by-row in the array a. to identify the individual nonzero entries in each row, we need to know in which column each entry lies. the column indices which correspond to the nonzero entries of m are stored in the array ja. i.e., if a(k) = m(i,j), then ja(k) = j. in addition, we need to know where each row starts and how long it is. the index positions in ja and a where the rows of m begin are stored in the array ia. i.e., if m(i,j) is the first nonzero entry (stored) in the i-th row and a(k) = m(i,j), then ia(i) = k. moreover, the index in ja and a of the first location following the last element in the last row is stored in ia(n+1). thus, the number of entries in the i-th row is given by ia(i+1) - ia(i), the nonzero entries of the i-th row are stored consecutively in

            a(ia(i)),  a(ia(i)+1),  ..., a(ia(i+1)-1),

and the corresponding column indices are stored consecutively in

            ja(ia(i)), ja(ia(i)+1), ..., ja(ia(i+1)-1).

for example, the 5 by 5 matrix

                ( 1. 0. 2. 0. 0.)
                ( 0. 3. 0. 0. 0.)
            m = ( 0. 4. 5. 6. 0.)
                ( 0. 0. 0. 7. 0.)
                ( 0. 0. 0. 8. 9.)

would be stored as

               - 1  2  3  4  5  6  7  8  9
            ---+--------------------------
            ia - 1  3  4  7  8 10
            ja - 1  3  2  2  3  4  4  4  5
             a - 1. 2. 3. 4. 5. 6. 7. 8. 9.         .

     nv    - n     - number of variables/equations.
     fva   - a     - nonzero entries of the coefficient matrix m, stored
           -           by rows.
           -           size = number of nonzero entries in m.
     nva   - ia    - pointers to delimit the rows in a.
           -           size = n+1.
     nva   - ja    - column numbers corresponding to the elements of a.
           -           size = size of a.
     fva   - b     - right-hand side b.  b and z can the same array.
           -           size = n.
     fra   - z     - solution x.  b and z can be the same array.
           -           size = n.

the rows and columns of the original matrix m can be reordered (e.g., to reduce fillin or ensure numerical stability) before calling the driver. if no reordering is done, then set r(i) = c(i) = ic(i) = i for i=1,…,n. the solution z is returned in the original order.

if the columns have been reordered (i.e., c(i).ne.i for some i), then the driver will call a subroutine (nroc) which rearranges each row of ja and a, leaving the rows in the original order, but placing the elements of each row in increasing order with respect to the new ordering. if path.ne.1, then nroc is assumed to have been called already.

 nva   - r     - ordering of the rows of m.
       -           size = n.
 nva   - c     - ordering of the columns of m.
       -           size = n.
 nva   - ic    - inverse of the ordering of the columns of m.  i.e.,
       -           ic(c(i)) = i  for i=1,...,n.
       -           size = n.

the solution of the system of linear equations is divided into three stages –

      nsfc -- the matrix m is processed symbolically to determine where
               fillin will occur during the numeric factorization.
      nnfc -- the matrix m is factored numerically into the product ldu
               of a unit lower triangular matrix l, a diagonal matrix
               d, and a unit upper triangular matrix u, and the system
               mx = b  is solved.
      nnsc -- the linear system  mx = b  is solved using the ldu
  or           factorization from nnfc.
      nntc -- the transposed linear system  mt x = b  is solved using
               the ldu factorization from nnf.

for several systems whose coefficient matrices have the same nonzero structure, nsfc need be done only once (for the first system). then nnfc is done once for each additional system. for several systems with the same coefficient matrix, nsfc and nnfc need be done only once (for the first system). then nnsc or nntc is done once for each additional right-hand side.

 nv    - path  - path specification.  values and their meanings are --
       -           1  perform nroc, nsfc, and nnfc.
       -           2  perform nnfc only  (nsfc is assumed to have been
       -               done in a manner compatible with the storage
       -               allocation used in the driver).
       -           3  perform nnsc only  (nsfc and nnfc are assumed to
       -               have been done in a manner compatible with the
       -               storage allocation used in the driver).
       -           4  perform nntc only  (nsfc and nnfc are assumed to
       -               have been done in a manner compatible with the
       -               storage allocation used in the driver).
       -           5  perform nroc and nsfc.

various errors are detected by the driver and the individual subroutines.

 nr    - flag  - error flag.  values and their meanings are --
       -             0     no errors detected
       -             n+k   null row in a  --  row = k
       -            2n+k   duplicate entry in a  --  row = k
       -            3n+k   insufficient storage in nsfc  --  row = k
       -            4n+1   insufficient storage in nnfc
       -            5n+k   null pivot  --  row = k
       -            6n+k   insufficient storage in nsfc  --  row = k
       -            7n+1   insufficient storage in nnfc
       -            8n+k   zero pivot  --  row = k
       -           10n+1   insufficient storage in cdrv
       -           11n+1   illegal path specification

working storage is needed for the factored form of the matrix m plus various temporary vectors. the arrays isp and rsp should be equivalenced. integer storage is allocated from the beginning of isp and real storage from the end of rsp.

        nv    - nsp   - declared dimension of rsp.  nsp generally must
                        be larger than  8n+2 + 2k  (where  k = (number of
                        nonzero entries in m)).
        nvira - isp   - integer working storage divided up into various arrays
                        needed by the subroutines.  isp and rsp should be
                        equivalenced.
                             size = lratio*nsp.
        fvira - rsp   - real working storage divided up into various arrays
                        needed by the subroutines.  isp and rsp should be
                        equivalenced.
                             size = nsp.
        nr    - esp   - if sufficient storage was available to perform the
                        symbolic factorization (nsfc), then esp is set to
                        the amount of excess storage provided (negative if
                        insufficient storage was available to perform the
                        numeric factorization (nnfc)).

Arguments

Type IntentOptional Attributes Name
integer :: N
integer :: R(*)
integer :: C(*)
integer :: Ic(*)
integer :: Ia(*)
integer :: Ja(*)
real(kind=dp) :: A(*)
real(kind=dp) :: B(*)
real(kind=dp) :: Z(*)
integer, intent(in) :: Nsp
integer, intent(inout) :: Isp(*)
real(kind=dp) :: Rsp(*)
integer, intent(out) :: Esp
integer, intent(in) :: Path
integer, intent(inout) :: Flag

Calls

proc~~cdrv~2~~CallsGraph proc~cdrv~2 cdrv.inc::cdrv nnfc nnfc proc~cdrv~2->nnfc nnsc nnsc proc~cdrv~2->nnsc nntc nntc proc~cdrv~2->nntc nroc nroc proc~cdrv~2->nroc nsfc nsfc proc~cdrv~2->nsfc

Variables

Type Visibility Attributes Name Initial
integer, public :: ar
integer, public :: d
integer, public :: i
integer, public :: ijl
integer, public :: iju
integer, public :: il
integer, public :: ira
integer, public :: irac
integer, public :: irl
integer, public :: iru
integer, public :: iu
integer, public :: j
integer, public :: jl
integer, public :: jlmax
integer, public :: jra
integer, public :: jrl
integer, public :: jru
integer, public :: ju
integer, public :: jumax
integer, public :: jutmp
integer, public :: l
integer, public :: lmax
integer, public, save :: lratio
integer, public :: max
integer, public :: q
integer, public :: row
integer, public :: tmp
integer, public :: u
integer, public :: umax

Source Code

subroutine cdrv(N,R,C,Ic,Ia,Ja,A,B,Z,Nsp,Isp,Rsp,Esp,Path,Flag)
!
integer               :: N
integer               :: R(*)
integer               :: C(*)
integer               :: Ic(*)
integer               :: Ia(*)
integer               :: Ja(*)
real(kind=dp)         :: A(*)
real(kind=dp)         :: B(*)
real(kind=dp)         :: Z(*)
integer,intent(inout) :: Isp(*)
integer,intent(in)    :: Nsp
real(kind=dp)         :: Rsp(*)
integer,intent(out)   :: Esp
integer,intent(in)    :: Path
integer,intent(inout) :: Flag
!
integer :: ar, d, i, ijl, iju, il, ira, irac, irl, iru, iu, j, jl, jlmax, jra, jrl, jru, ju, jumax, jutmp, l,  &
         & lmax, max, q, row, tmp, u, umax
integer , save :: lratio
!
!   set lratio equal to the ratio between the length of floating point
!   and integer array data.  e. g., lratio = 1 for (real, integer),
!   lratio = 2 for (double precision, integer)
!
data lratio/2/
!
if ( Path<1 .or. 5<Path ) then
!  ** error.. illegal path specification
   Flag = 11*N + 1
   return
else
! ### ***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
!
!   ******  reorder a if necessary, call nsfc if flag is set  ***********
   if ( (Path-1)*(Path-5)==0 ) then
      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<=0 .or. jumax<=0 )then
         !  ** error.. insufficient storage
         Flag = 10*N + 1
         return
      endif
!
      do i = 1 , N
         if ( C(i)/=i ) then
            ar = Nsp + 1 - N
            call nroc(N,Ic,Ia,Ja,A,Isp(il),Rsp(ar),Isp(iu),Flag)
            if ( Flag==0 ) exit
               !  ** error.. error detected in nroc, nsfc, nnfc, or nnsc
               return
         endif
      enddo
!
      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/=0 )then
         !  ** error.. error detected in nroc, nsfc, nnfc, or nnsc
         return
      endif
!   ******  move ju next to jl  *****************************************
      jlmax = Isp(ijl+N-1)
      ju = jl + jlmax
      jumax = Isp(iju+N-1)
      if ( jumax>0 ) then
         do j = 1 , jumax
            Isp(ju+j-1) = Isp(jutmp+j-1)
         enddo
      endif
   endif
!
!   ******  call remaining subroutines  *********************************
   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)
!
   if ( (Path-1)*(Path-2)==0 ) then
      if ( umax<0 ) then
         !  ** error.. insufficient storage
         Flag = 10*N + 1
         return
      endif
      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/=0 ) then
         !  ** error.. error detected in nroc, nsfc, nnfc, or nnsc
         return
      endif
   endif
!
   if ( (Path-3)==0 ) 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))
!
   if ( (Path-4)==0 ) 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))
   return
endif
!
end subroutine cdrv