← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/interpolate/fitpack/fpgrpa.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/interpolate/fitpack/fpgrpa.f
1
 
      subroutine fpgrpa(ifsu,ifsv,ifbu,ifbv,idim,ipar,u,mu,v,mv,z,mz,
2
 
     * tu,nu,tv,nv,p,c,nc,fp,fpu,fpv,mm,mvnu,spu,spv,right,q,au,au1,
3
 
     * av,av1,bu,bv,nru,nrv)
4
 
c  ..
5
 
c  ..scalar arguments..
6
 
      real*8 p,fp
7
 
      integer ifsu,ifsv,ifbu,ifbv,idim,mu,mv,mz,nu,nv,nc,mm,mvnu
8
 
c  ..array arguments..
9
 
      real*8 u(mu),v(mv),z(mz*idim),tu(nu),tv(nv),c(nc*idim),fpu(nu),
10
 
     * fpv(nv),spu(mu,4),spv(mv,4),right(mm*idim),q(mvnu),au(nu,5),
11
 
     * au1(nu,4),av(nv,5),av1(nv,4),bu(nu,5),bv(nv,5)
12
 
      integer ipar(2),nru(mu),nrv(mv)
13
 
c  ..local scalars..
14
 
      real*8 arg,fac,term,one,half,value
15
 
      integer i,id,ii,it,iz,i1,i2,j,jz,k,k1,k2,l,l1,l2,mvv,k0,muu,
16
 
     * ncof,nroldu,nroldv,number,nmd,numu,numu1,numv,numv1,nuu,nvv,
17
 
     * nu4,nu7,nu8,nv4,nv7,nv8
18
 
c  ..local arrays..
19
 
      real*8 h(5)
20
 
c  ..subroutine references..
21
 
c    fpback,fpbspl,fpdisc,fpbacp,fptrnp,fptrpe
22
 
c  ..
23
 
c  let
24
 
c               |   (spu)    |            |   (spv)    |
25
 
c        (au) = | ---------- |     (av) = | ---------- |
26
 
c               | (1/p) (bu) |            | (1/p) (bv) |
27
 
c
28
 
c                                | z  ' 0 |
29
 
c                            q = | ------ |
30
 
c                                | 0  ' 0 |
31
 
c
32
 
c  with c      : the (nu-4) x (nv-4) matrix which contains the b-spline
33
 
c                coefficients.
34
 
c       z      : the mu x mv matrix which contains the function values.
35
 
c       spu,spv: the mu x (nu-4), resp. mv x (nv-4) observation matrices
36
 
c                according to the least-squares problems in the u-,resp.
37
 
c                v-direction.
38
 
c       bu,bv  : the (nu-7) x (nu-4),resp. (nv-7) x (nv-4) matrices
39
 
c                containing the discontinuity jumps of the derivatives
40
 
c                of the b-splines in the u-,resp.v-variable at the knots
41
 
c  the b-spline coefficients of the smoothing spline are then calculated
42
 
c  as the least-squares solution of the following over-determined linear
43
 
c  system of equations
44
 
c
45
 
c    (1)  (av) c (au)' = q
46
 
c
47
 
c  subject to the constraints
48
 
c
49
 
c    (2)  c(nu-3+i,j) = c(i,j), i=1,2,3 ; j=1,2,...,nv-4
50
 
c            if(ipar(1).ne.0)
51
 
c
52
 
c    (3)  c(i,nv-3+j) = c(i,j), j=1,2,3 ; i=1,2,...,nu-4
53
 
c            if(ipar(2).ne.0)
54
 
c
55
 
c  set constants
56
 
      one = 1
57
 
      half = 0.5
58
 
c  initialization
59
 
      nu4 = nu-4
60
 
      nu7 = nu-7
61
 
      nu8 = nu-8
62
 
      nv4 = nv-4
63
 
      nv7 = nv-7
64
 
      nv8 = nv-8
65
 
      muu = mu
66
 
      if(ipar(1).ne.0) muu = mu-1
67
 
      mvv = mv
68
 
      if(ipar(2).ne.0) mvv = mv-1
69
 
c  it depends on the value of the flags ifsu,ifsv,ifbu  and ibvand
70
 
c  on the value of p whether the matrices (spu), (spv), (bu) and (bv)
71
 
c  still must be determined.
72
 
      if(ifsu.ne.0) go to 50
73
 
c  calculate the non-zero elements of the matrix (spu) which is the ob-
74
 
c  servation matrix according to the least-squares spline approximation
75
 
c  problem in the u-direction.
76
 
      l = 4
77
 
      l1 = 5
78
 
      number = 0
79
 
      do 40 it=1,muu
80
 
        arg = u(it)
81
 
  10    if(arg.lt.tu(l1) .or. l.eq.nu4) go to 20
82
 
        l = l1
83
 
        l1 = l+1
84
 
        number = number+1
85
 
        go to 10
86
 
  20    call fpbspl(tu,nu,3,arg,l,h)
87
 
        do 30 i=1,4
88
 
          spu(it,i) = h(i)
89
 
  30    continue
90
 
        nru(it) = number
91
 
  40  continue
92
 
      ifsu = 1
93
 
c  calculate the non-zero elements of the matrix (spv) which is the ob-
94
 
c  servation matrix according to the least-squares spline approximation
95
 
c  problem in the v-direction.
96
 
  50  if(ifsv.ne.0) go to 100
97
 
      l = 4
98
 
      l1 = 5
99
 
      number = 0
100
 
      do 90 it=1,mvv
101
 
        arg = v(it)
102
 
  60    if(arg.lt.tv(l1) .or. l.eq.nv4) go to 70
103
 
        l = l1
104
 
        l1 = l+1
105
 
        number = number+1
106
 
        go to 60
107
 
  70    call fpbspl(tv,nv,3,arg,l,h)
108
 
        do 80 i=1,4
109
 
          spv(it,i) = h(i)
110
 
  80    continue
111
 
        nrv(it) = number
112
 
  90  continue
113
 
      ifsv = 1
114
 
 100  if(p.le.0.) go to  150
115
 
c  calculate the non-zero elements of the matrix (bu).
116
 
      if(ifbu.ne.0 .or. nu8.eq.0) go to 110
117
 
      call fpdisc(tu,nu,5,bu,nu)
118
 
      ifbu = 1
119
 
c  calculate the non-zero elements of the matrix (bv).
120
 
 110  if(ifbv.ne.0 .or. nv8.eq.0) go to 150
121
 
      call fpdisc(tv,nv,5,bv,nv)
122
 
      ifbv = 1
123
 
c  substituting (2)  and (3) into (1), we obtain the overdetermined
124
 
c  system
125
 
c         (4)  (avv) (cr) (auu)' = (qq)
126
 
c  from which the nuu*nvv remaining coefficients
127
 
c         c(i,j) , i=1,...,nu-4-3*ipar(1) ; j=1,...,nv-4-3*ipar(2) ,
128
 
c  the elements of (cr), are then determined in the least-squares sense.
129
 
c  we first determine the matrices (auu) and (qq). then we reduce the
130
 
c  matrix (auu) to upper triangular form (ru) using givens rotations.
131
 
c  we apply the same transformations to the rows of matrix qq to obtain
132
 
c  the (mv) x nuu matrix g.
133
 
c  we store matrix (ru) into au (and au1 if ipar(1)=1) and g into q.
134
 
 150  if(ipar(1).ne.0) go to 160
135
 
      nuu = nu4
136
 
      call fptrnp(mu,mv,idim,nu,nru,spu,p,bu,z,au,q,right)
137
 
      go to 180
138
 
 160  nuu = nu7
139
 
      call fptrpe(mu,mv,idim,nu,nru,spu,p,bu,z,au,au1,q,right)
140
 
c  we determine the matrix (avv) and then we reduce this matrix to
141
 
c  upper triangular form (rv) using givens rotations.
142
 
c  we apply the same transformations to the columns of matrix
143
 
c  g to obtain the (nvv) x (nuu) matrix h.
144
 
c  we store matrix (rv) into av (and av1 if ipar(2)=1) and h into c.
145
 
 180  if(ipar(2).ne.0) go to 190
146
 
      nvv = nv4
147
 
      call fptrnp(mv,nuu,idim,nv,nrv,spv,p,bv,q,av,c,right)
148
 
      go to 200
149
 
 190  nvv = nv7
150
 
      call fptrpe(mv,nuu,idim,nv,nrv,spv,p,bv,q,av,av1,c,right)
151
 
c  backward substitution to obtain the b-spline coefficients as the
152
 
c  solution of the linear system    (rv) (cr) (ru)' = h.
153
 
c  first step: solve the system  (rv) (c1) = h.
154
 
 200  ncof = nuu*nvv
155
 
      k = 1
156
 
      if(ipar(2).ne.0) go to 240
157
 
      do 220 ii=1,idim
158
 
      do 220 i=1,nuu
159
 
         call fpback(av,c(k),nvv,5,c(k),nv)
160
 
         k = k+nvv
161
 
 220  continue
162
 
      go to 300
163
 
 240  do 260 ii=1,idim
164
 
      do 260 i=1,nuu
165
 
         call fpbacp(av,av1,c(k),nvv,4,c(k),5,nv)
166
 
         k = k+nvv
167
 
 260  continue
168
 
c  second step: solve the system  (cr) (ru)' = (c1).
169
 
 300  if(ipar(1).ne.0) go to 400
170
 
      do 360 ii=1,idim
171
 
      k = (ii-1)*ncof
172
 
      do 360 j=1,nvv
173
 
        k = k+1
174
 
        l = k
175
 
        do 320 i=1,nuu
176
 
          right(i) = c(l)
177
 
          l = l+nvv
178
 
 320    continue
179
 
        call fpback(au,right,nuu,5,right,nu)
180
 
        l = k
181
 
        do 340 i=1,nuu
182
 
           c(l) = right(i)
183
 
           l = l+nvv
184
 
 340    continue
185
 
 360  continue
186
 
      go to 500
187
 
 400  do 460 ii=1,idim
188
 
      k = (ii-1)*ncof
189
 
      do 460 j=1,nvv
190
 
        k = k+1
191
 
        l = k
192
 
        do 420 i=1,nuu
193
 
          right(i) = c(l)
194
 
          l = l+nvv
195
 
 420    continue
196
 
        call fpbacp(au,au1,right,nuu,4,right,5,nu)
197
 
        l = k
198
 
        do 440 i=1,nuu
199
 
           c(l) = right(i)
200
 
           l = l+nvv
201
 
 440    continue
202
 
 460  continue
203
 
c  calculate from the conditions (2)-(3), the remaining b-spline
204
 
c  coefficients.
205
 
 500  if(ipar(2).eq.0) go to 600
206
 
      i = 0
207
 
      j = 0
208
 
      do 560 id=1,idim
209
 
      do 560 l=1,nuu
210
 
         ii = i
211
 
         do 520 k=1,nvv
212
 
            i = i+1
213
 
            j = j+1
214
 
            q(i) = c(j)
215
 
 520     continue
216
 
         do 540 k=1,3
217
 
            ii = ii+1
218
 
            i = i+1
219
 
            q(i) = q(ii)
220
 
 540     continue
221
 
 560  continue
222
 
      ncof = nv4*nuu
223
 
      nmd = ncof*idim
224
 
      do 580 i=1,nmd
225
 
         c(i) = q(i)
226
 
 580  continue
227
 
 600  if(ipar(1).eq.0) go to 700
228
 
      i = 0
229
 
      j = 0
230
 
      n33 = 3*nv4
231
 
      do 660 id=1,idim
232
 
         ii = i
233
 
         do 620 k=1,ncof
234
 
            i = i+1
235
 
            j = j+1
236
 
            q(i) = c(j)
237
 
 620     continue
238
 
         do 640 k=1,n33
239
 
            ii = ii+1
240
 
            i = i+1
241
 
            q(i) = q(ii)
242
 
 640     continue
243
 
 660  continue
244
 
      ncof = nv4*nu4
245
 
      nmd = ncof*idim
246
 
      do 680 i=1,nmd
247
 
         c(i) = q(i)
248
 
 680  continue
249
 
c  calculate the quantities
250
 
c    res(i,j) = (z(i,j) - s(u(i),v(j)))**2 , i=1,2,..,mu;j=1,2,..,mv
251
 
c    fp = sumi=1,mu(sumj=1,mv(res(i,j)))
252
 
c    fpu(r) = sum''i(sumj=1,mv(res(i,j))) , r=1,2,...,nu-7
253
 
c                  tu(r+3) <= u(i) <= tu(r+4)
254
 
c    fpv(r) = sumi=1,mu(sum''j(res(i,j))) , r=1,2,...,nv-7
255
 
c                  tv(r+3) <= v(j) <= tv(r+4)
256
 
 700  fp = 0.
257
 
      do 720 i=1,nu
258
 
        fpu(i) = 0.
259
 
 720  continue
260
 
      do 740 i=1,nv
261
 
        fpv(i) = 0.
262
 
 740  continue
263
 
      nroldu = 0
264
 
c  main loop for the different grid points.
265
 
      do 860 i1=1,muu
266
 
        numu = nru(i1)
267
 
        numu1 = numu+1
268
 
        nroldv = 0
269
 
        iz = (i1-1)*mv
270
 
        do 840 i2=1,mvv
271
 
          numv = nrv(i2)
272
 
          numv1 = numv+1
273
 
          iz = iz+1
274
 
c  evaluate s(u,v) at the current grid point by making the sum of the
275
 
c  cross products of the non-zero b-splines at (u,v), multiplied with
276
 
c  the appropiate b-spline coefficients.
277
 
          term = 0.
278
 
          k0 = numu*nv4+numv
279
 
          jz = iz
280
 
          do 800 id=1,idim
281
 
          k1 = k0
282
 
          value = 0.
283
 
          do 780 l1=1,4
284
 
            k2 = k1
285
 
            fac = spu(i1,l1)
286
 
            do 760 l2=1,4
287
 
              k2 = k2+1
288
 
              value = value+fac*spv(i2,l2)*c(k2)
289
 
 760        continue
290
 
            k1 = k1+nv4
291
 
 780      continue
292
 
c  calculate the squared residual at the current grid point.
293
 
          term = term+(z(jz)-value)**2
294
 
          jz = jz+mz
295
 
          k0 = k0+ncof
296
 
 800      continue
297
 
c  adjust the different parameters.
298
 
          fp = fp+term
299
 
          fpu(numu1) = fpu(numu1)+term
300
 
          fpv(numv1) = fpv(numv1)+term
301
 
          fac = term*half
302
 
          if(numv.eq.nroldv) go to 820
303
 
          fpv(numv1) = fpv(numv1)-fac
304
 
          fpv(numv) = fpv(numv)+fac
305
 
 820      nroldv = numv
306
 
          if(numu.eq.nroldu) go to 840
307
 
          fpu(numu1) = fpu(numu1)-fac
308
 
          fpu(numu) = fpu(numu)+fac
309
 
 840    continue
310
 
        nroldu = numu
311
 
 860  continue
312
 
      return
313
 
      end