share/odepack/fortran/sro.f

   1       subroutine sro
   2      *     (n, ip, ia,ja,a, q, r, dflag)
   3 c***********************************************************************
   4 c  sro -- symmetric reordering of sparse symmetric matrix
   5 c***********************************************************************
   6 c
   7 c  description
   8 c
   9 c    the nonzero entries of the matrix m are assumed to be stored
  10 c    symmetrically in (ia,ja,a) format (i.e., not both m(i,j) and m(j,i)
  11 c    are stored if i ne j).
  12 c
  13 c    sro does not rearrange the order of the rows, but does move
  14 c    nonzeroes from one row to another to ensure that if m(i,j) will be
  15 c    in the upper triangle of m with respect to the new ordering, then
  16 c    m(i,j) is stored in row i (and thus m(j,i) is not stored),  whereas
  17 c    if m(i,j) will be in the strict lower triangle of m, then m(j,i) is
  18 c    stored in row j (and thus m(i,j) is not stored).
  19 c
  20 c
  21 c  additional parameters
  22 c
  23 c    q     - integer one-dimensional work array.  dimension = n
  24 c
  25 c    r     - integer one-dimensional work array.  dimension = number of
  26 c            nonzero entries in the upper triangle of m
  27 c
  28 c    dflag - logical variable.  if dflag = .true., then store nonzero
  29 c            diagonal elements at the beginning of the row
  30 c
  31 c-----------------------------------------------------------------------
  32 c
  33       integer  ip(*),  ia(*), ja(*),  q(*), r(*)
  34 c...  real  a(*),  ak
  35       double precision  a(*),  ak
  36       logical  dflag
  37 c
  38 c
  39 c--phase 1 -- find row in which to store each nonzero
  40 c----initialize count of nonzeroes to be stored in each row
  41       do 1 i=1,n
  42   1     q(i) = 0
  43 c
  44 c----for each nonzero element a(j)
  45       do 3 i=1,n
  46         jmin = ia(i)
  47         jmax = ia(i+1) - 1
  48         if (jmin.gt.jmax)  go to 3
  49         do 2 j=jmin,jmax
  50 c
  51 c--------find row (=r(j)) and column (=ja(j)) in which to store a(j) ...
  52           k = ja(j)
  53           if (ip(k).lt.ip(i))  ja(j) = i
  54           if (ip(k).ge.ip(i))  k = i
  55           r(j) = k
  56 c
  57 c--------... and increment count of nonzeroes (=q(r(j)) in that row
  58   2       q(k) = q(k) + 1
  59   3     continue
  60 c
  61 c
  62 c--phase 2 -- find new ia and permutation to apply to (ja,a)
  63 c----determine pointers to delimit rows in permuted (ja,a)
  64       do 4 i=1,n
  65         ia(i+1) = ia(i) + q(i)
  66   4     q(i) = ia(i+1)
  67 c
  68 c----determine where each (ja(j),a(j)) is stored in permuted (ja,a)
  69 c----for each nonzero element (in reverse order)
  70       ilast = 0
  71       jmin = ia(1)
  72       jmax = ia(n+1) - 1
  73       j = jmax
  74       do 6 jdummy=jmin,jmax
  75         i = r(j)
  76         if (.not.dflag .or. ja(j).ne.i .or. i.eq.ilast)  go to 5
  77 c
  78 c------if dflag, then put diagonal nonzero at beginning of row
  79           r(j) = ia(i)
  80           ilast = i
  81           go to 6
  82 c
  83 c------put (off-diagonal) nonzero in last unused location in row
  84   5       q(i) = q(i) - 1
  85           r(j) = q(i)
  86 c
  87   6     j = j-1
  88 c
  89 c
  90 c--phase 3 -- permute (ja,a) to upper triangular form (wrt new ordering)
  91       do 8 j=jmin,jmax
  92   7     if (r(j).eq.j)  go to 8
  93           k = r(j)
  94           r(j) = r(k)
  95           r(k) = k
  96           jak = ja(k)
  97           ja(k) = ja(j)
  98           ja(j) = jak
  99           ak = a(k)
 100           a(k) = a(j)
 101           a(j) = ak
 102           go to 7
 103   8     continue
 104 c
 105       return
 106       end