← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/integrate/odepack/nnfc.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/nnfc.f
1
 
      subroutine nnfc
2
 
     *     (n, r,c,ic, ia,ja,a, z, b,
3
 
     *      lmax,il,jl,ijl,l, d, umax,iu,ju,iju,u,
4
 
     *      row, tmp, irl,jrl, flag)
5
 
clll. optimize
6
 
c*** subroutine nnfc
7
 
c*** numerical ldu-factorization of sparse nonsymmetric matrix and
8
 
c      solution of system of linear equations (compressed pointer
9
 
c      storage)
10
 
c
11
 
c
12
 
c       input variables..  n, r, c, ic, ia, ja, a, b,
13
 
c                          il, jl, ijl, lmax, iu, ju, iju, umax
14
 
c       output variables.. z, l, d, u, flag
15
 
c
16
 
c       parameters used internally..
17
 
c nia   - irl,  - vectors used to find the rows of  l.  at the kth step
18
 
c nia   - jrl       of the factorization,  jrl(k)  points to the head
19
 
c       -           of a linked list in  jrl  of column indices j
20
 
c       -           such j .lt. k and  l(k,j)  is nonzero.  zero
21
 
c       -           indicates the end of the list.  irl(j)  (j.lt.k)
22
 
c       -           points to the smallest i such that i .ge. k and
23
 
c       -           l(i,j)  is nonzero.
24
 
c       -           size of each = n.
25
 
c fia   - row   - holds intermediate values in calculation of  u and l.
26
 
c       -           size = n.
27
 
c fia   - tmp   - holds new right-hand side  b*  for solution of the
28
 
c       -           equation ux = b*.
29
 
c       -           size = n.
30
 
c
31
 
c  internal variables..
32
 
c    jmin, jmax - indices of the first and last positions in a row to
33
 
c      be examined.
34
 
c    sum - used in calculating  tmp.
35
 
c
36
 
      integer rk,umax
37
 
      integer  r(1), c(1), ic(1), ia(1), ja(1), il(1), jl(1), ijl(1)
38
 
      integer  iu(1), ju(1), iju(1), irl(1), jrl(1), flag
39
 
      double precision  a(1), l(1), d(1), u(1), z(1), b(1), row(1)
40
 
      double precision  tmp(1), lki, sum, dk
41
 
c
42
 
c  ******  initialize pointers and test storage  ***********************
43
 
      if(il(n+1)-1 .gt. lmax) go to 104
44
 
      if(iu(n+1)-1 .gt. umax) go to 107
45
 
      do 1 k=1,n
46
 
        irl(k) = il(k)
47
 
        jrl(k) = 0
48
 
   1    continue
49
 
c
50
 
c  ******  for each row  ***********************************************
51
 
      do 19 k=1,n
52
 
c  ******  reverse jrl and zero row where kth row of l will fill in  ***
53
 
        row(k) = 0
54
 
        i1 = 0
55
 
        if (jrl(k) .eq. 0) go to 3
56
 
        i = jrl(k)
57
 
   2    i2 = jrl(i)
58
 
        jrl(i) = i1
59
 
        i1 = i
60
 
        row(i) = 0
61
 
        i = i2
62
 
        if (i .ne. 0) go to 2
63
 
c  ******  set row to zero where u will fill in  ***********************
64
 
   3    jmin = iju(k)
65
 
        jmax = jmin + iu(k+1) - iu(k) - 1
66
 
        if (jmin .gt. jmax) go to 5
67
 
        do 4 j=jmin,jmax
68
 
   4      row(ju(j)) = 0
69
 
c  ******  place kth row of a in row  **********************************
70
 
   5    rk = r(k)
71
 
        jmin = ia(rk)
72
 
        jmax = ia(rk+1) - 1
73
 
        do 6 j=jmin,jmax
74
 
          row(ic(ja(j))) = a(j)
75
 
   6      continue
76
 
c  ******  initialize sum, and link through jrl  ***********************
77
 
        sum = b(rk)
78
 
        i = i1
79
 
        if (i .eq. 0) go to 10
80
 
c  ******  assign the kth row of l and adjust row, sum  ****************
81
 
   7      lki = -row(i)
82
 
c  ******  if l is not required, then comment out the following line  **
83
 
          l(irl(i)) = -lki
84
 
          sum = sum + lki * tmp(i)
85
 
          jmin = iu(i)
86
 
          jmax = iu(i+1) - 1
87
 
          if (jmin .gt. jmax) go to 9
88
 
          mu = iju(i) - jmin
89
 
          do 8 j=jmin,jmax
90
 
   8        row(ju(mu+j)) = row(ju(mu+j)) + lki * u(j)
91
 
   9      i = jrl(i)
92
 
          if (i .ne. 0) go to 7
93
 
c
94
 
c  ******  assign kth row of u and diagonal d, set tmp(k)  *************
95
 
  10    if (row(k) .eq. 0.0d0) go to 108
96
 
        dk = 1.0d0 / row(k)
97
 
        d(k) = dk
98
 
        tmp(k) = sum * dk
99
 
        if (k .eq. n) go to 19
100
 
        jmin = iu(k)
101
 
        jmax = iu(k+1) - 1
102
 
        if (jmin .gt. jmax)  go to 12
103
 
        mu = iju(k) - jmin
104
 
        do 11 j=jmin,jmax
105
 
  11      u(j) = row(ju(mu+j)) * dk
106
 
  12    continue
107
 
c
108
 
c  ******  update irl and jrl, keeping jrl in decreasing order  ********
109
 
        i = i1
110
 
        if (i .eq. 0) go to 18
111
 
  14    irl(i) = irl(i) + 1
112
 
        i1 = jrl(i)
113
 
        if (irl(i) .ge. il(i+1)) go to 17
114
 
        ijlb = irl(i) - il(i) + ijl(i)
115
 
        j = jl(ijlb)
116
 
  15    if (i .gt. jrl(j)) go to 16
117
 
          j = jrl(j)
118
 
          go to 15
119
 
  16    jrl(i) = jrl(j)
120
 
        jrl(j) = i
121
 
  17    i = i1
122
 
        if (i .ne. 0) go to 14
123
 
  18    if (irl(k) .ge. il(k+1)) go to 19
124
 
        j = jl(ijl(k))
125
 
        jrl(k) = jrl(j)
126
 
        jrl(j) = k
127
 
  19    continue
128
 
c
129
 
c  ******  solve  ux = tmp  by back substitution  **********************
130
 
      k = n
131
 
      do 22 i=1,n
132
 
        sum =  tmp(k)
133
 
        jmin = iu(k)
134
 
        jmax = iu(k+1) - 1
135
 
        if (jmin .gt. jmax)  go to 21
136
 
        mu = iju(k) - jmin
137
 
        do 20 j=jmin,jmax
138
 
  20      sum = sum - u(j) * tmp(ju(mu+j))
139
 
  21    tmp(k) =  sum
140
 
        z(c(k)) =  sum
141
 
  22    k = k-1
142
 
      flag = 0
143
 
      return
144
 
c
145
 
c ** error.. insufficient storage for l
146
 
 104  flag = 4*n + 1
147
 
      return
148
 
c ** error.. insufficient storage for u
149
 
 107  flag = 7*n + 1
150
 
      return
151
 
c ** error.. zero pivot
152
 
 108  flag = 8*n + k
153
 
      return
154
 
      end