← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/integrate/odepack/sro.f

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/integrate/odepack/sro.f
1
 
      subroutine sro
2
 
     *     (n, ip, ia,ja,a, q, r, dflag)
3
 
clll. optimize
4
 
c***********************************************************************
5
 
c  sro -- symmetric reordering of sparse symmetric matrix
6
 
c***********************************************************************
7
 
c
8
 
c  description
9
 
c
10
 
c    the nonzero entries of the matrix m are assumed to be stored
11
 
c    symmetrically in (ia,ja,a) format (i.e., not both m(i,j) and m(j,i)
12
 
c    are stored if i ne j).
13
 
c
14
 
c    sro does not rearrange the order of the rows, but does move
15
 
c    nonzeroes from one row to another to ensure that if m(i,j) will be
16
 
c    in the upper triangle of m with respect to the new ordering, then
17
 
c    m(i,j) is stored in row i (and thus m(j,i) is not stored),  whereas
18
 
c    if m(i,j) will be in the strict lower triangle of m, then m(j,i) is
19
 
c    stored in row j (and thus m(i,j) is not stored).
20
 
c
21
 
c
22
 
c  additional parameters
23
 
c
24
 
c    q     - integer one-dimensional work array.  dimension = n
25
 
c
26
 
c    r     - integer one-dimensional work array.  dimension = number of
27
 
c            nonzero entries in the upper triangle of m
28
 
c
29
 
c    dflag - logical variable.  if dflag = .true., then store nonzero
30
 
c            diagonal elements at the beginning of the row
31
 
c
32
 
c-----------------------------------------------------------------------
33
 
c
34
 
      integer  ip(1),  ia(1), ja(1),  q(1), r(1)
35
 
      double precision  a(1),  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