share/odepack/fortran/odrv.f

   1 *DECK ODRV
   2       subroutine odrv
   3      *     (n, ia,ja,a, p,ip, nsp,isp, path, flag)
   4 c                                                                 5/2/83
   5 c***********************************************************************
   6 c  odrv -- driver for sparse matrix reordering routines
   7 c***********************************************************************
   8 c
   9 c  description
  10 c
  11 c    odrv finds a minimum degree ordering of the rows and columns
  12 c    of a matrix m stored in (ia,ja,a) format (see below).  for the
  13 c    reordered matrix, the work and storage required to perform
  14 c    gaussian elimination is (usually) significantly less.
  15 c
  16 c    note.. odrv and its subordinate routines have been modified to
  17 c    compute orderings for general matrices, not necessarily having any
  18 c    symmetry.  the miminum degree ordering is computed for the
  19 c    structure of the symmetric matrix  m + m-transpose.
  20 c    modifications to the original odrv module have been made in
  21 c    the coding in subroutine mdi, and in the initial comments in
  22 c    subroutines odrv and md.
  23 c
  24 c    if only the nonzero entries in the upper triangle of m are being
  25 c    stored, then odrv symmetrically reorders (ia,ja,a), (optionally)
  26 c    with the diagonal entries placed first in each row.  this is to
  27 c    ensure that if m(i,j) will be in the upper triangle of m with
  28 c    respect to the new ordering, then m(i,j) is stored in row i (and
  29 c    thus m(j,i) is not stored),  whereas if m(i,j) will be in the
  30 c    strict lower triangle of m, then m(j,i) is stored in row j (and
  31 c    thus m(i,j) is not stored).
  32 c
  33 c
  34 c  storage of sparse matrices
  35 c
  36 c    the nonzero entries of the matrix m are stored row-by-row in the
  37 c    array a.  to identify the individual nonzero entries in each row,
  38 c    we need to know in which column each entry lies.  these column
  39 c    indices are stored in the array ja.  i.e., if  a(k) = m(i,j),  then
  40 c    ja(k) = j.  to identify the individual rows, we need to know where
  41 c    each row starts.  these row pointers are stored in the array ia.
  42 c    i.e., if m(i,j) is the first nonzero entry (stored) in the i-th row
  43 c    and  a(k) = m(i,j),  then  ia(i) = k.  moreover, ia(n+1) points to
  44 c    the first location following the last element in the last row.
  45 c    thus, the number of entries in the i-th row is  ia(i+1) - ia(i),
  46 c    the nonzero entries in the i-th row are stored consecutively in
  47 c
  48 c            a(ia(i)),  a(ia(i)+1),  ..., a(ia(i+1)-1),
  49 c
  50 c    and the corresponding column indices are stored consecutively in
  51 c
  52 c            ja(ia(i)), ja(ia(i)+1), ..., ja(ia(i+1)-1).
  53 c
  54 c    when the coefficient matrix is symmetric, only the nonzero entries
  55 c    in the upper triangle need be stored.  for example, the matrix
  56 c
  57 c             ( 1  0  2  3  0 )
  58 c             ( 0  4  0  0  0 )
  59 c         m = ( 2  0  5  6  0 )
  60 c             ( 3  0  6  7  8 )
  61 c             ( 0  0  0  8  9 )
  62 c
  63 c    could be stored as
  64 c
  65 c            - 1  2  3  4  5  6  7  8  9 10 11 12 13
  66 c         ---+--------------------------------------
  67 c         ia - 1  4  5  8 12 14
  68 c         ja - 1  3  4  2  1  3  4  1  3  4  5  4  5
  69 c          a - 1  2  3  4  2  5  6  3  6  7  8  8  9
  70 c
  71 c    or (symmetrically) as
  72 c
  73 c            - 1  2  3  4  5  6  7  8  9
  74 c         ---+--------------------------
  75 c         ia - 1  4  5  7  9 10
  76 c         ja - 1  3  4  2  3  4  4  5  5
  77 c          a - 1  2  3  4  5  6  7  8  9          .
  78 c
  79 c
  80 c  parameters
  81 c
  82 c    n    - order of the matrix
  83 c
  84 c    ia   - integer one-dimensional array containing pointers to delimit
  85 c           rows in ja and a.  dimension = n+1
  86 c
  87 c    ja   - integer one-dimensional array containing the column indices
  88 c           corresponding to the elements of a.  dimension = number of
  89 c           nonzero entries in (the upper triangle of) m
  90 c
  91 c    a    - real one-dimensional array containing the nonzero entries in
  92 c           (the upper triangle of) m, stored by rows.  dimension =
  93 c           number of nonzero entries in (the upper triangle of) m
  94 c
  95 c    p    - integer one-dimensional array used to return the permutation
  96 c           of the rows and columns of m corresponding to the minimum
  97 c           degree ordering.  dimension = n
  98 c
  99 c    ip   - integer one-dimensional array used to return the inverse of
 100 c           the permutation returned in p.  dimension = n
 101 c
 102 c    nsp  - declared dimension of the one-dimensional array isp.  nsp
 103 c           must be at least  3n+4k,  where k is the number of nonzeroes
 104 c           in the strict upper triangle of m
 105 c
 106 c    isp  - integer one-dimensional array used for working storage.
 107 c           dimension = nsp
 108 c
 109 c    path - integer path specification.  values and their meanings are -
 110 c             1  find minimum degree ordering only
 111 c             2  find minimum degree ordering and reorder symmetrically
 112 c                  stored matrix (used when only the nonzero entries in
 113 c                  the upper triangle of m are being stored)
 114 c             3  reorder symmetrically stored matrix as specified by
 115 c                  input permutation (used when an ordering has already
 116 c                  been determined and only the nonzero entries in the
 117 c                  upper triangle of m are being stored)
 118 c             4  same as 2 but put diagonal entries at start of each row
 119 c             5  same as 3 but put diagonal entries at start of each row
 120 c
 121 c    flag - integer error flag.  values and their meanings are -
 122 c               0    no errors detected
 123 c              9n+k  insufficient storage in md
 124 c             10n+1  insufficient storage in odrv
 125 c             11n+1  illegal path specification
 126 c
 127 c
 128 c  conversion from real to double precision
 129 c
 130 c    change the real declarations in odrv and sro to double precision
 131 c    declarations.
 132 c
 133 c-----------------------------------------------------------------------
 134 c
 135       integer  ia(*), ja(*),  p(*), ip(*),  isp(*),  path,  flag,
 136      *   v, l, head,  tmp, q
 137 c...  real  a(*)
 138       double precision  a(*)
 139       logical  dflag
 140 c
 141 c----initialize error flag and validate path specification
 142       flag = 0
 143       if (path.lt.1 .or. 5.lt.path)  go to 111
 144 c
 145 c----allocate storage and find minimum degree ordering
 146       if ((path-1) * (path-2) * (path-4) .ne. 0)  go to 1
 147         max = (nsp-n)/2
 148         v    = 1
 149         l    = v     +  max
 150         head = l     +  max
 151         next = head  +  n
 152         if (max.lt.n)  go to 110
 153 c
 154         call  md
 155      *     (n, ia,ja, max,isp(v),isp(l), isp(head),p,ip, isp(v), flag)
 156         if (flag.ne.0)  go to 100
 157 c
 158 c----allocate storage and symmetrically reorder matrix
 159    1  if ((path-2) * (path-3) * (path-4) * (path-5) .ne. 0)  go to 2
 160         tmp = (nsp+1) -      n
 161         q   = tmp     - (ia(n+1)-1)
 162         if (q.lt.1)  go to 110
 163 c
 164         dflag = path.eq.4 .or. path.eq.5
 165         call sro
 166      *     (n,  ip,  ia, ja, a,  isp(tmp),  isp(q),  dflag)
 167 c
 168    2  return
 169 c
 170 c ** error -- error detected in md
 171  100  return
 172 c ** error -- insufficient storage
 173  110  flag = 10*n + 1
 174       return
 175 c ** error -- illegal path specified
 176  111  flag = 11*n + 1
 177       return
 178       end