1      subroutine fpsurf(iopt,m,x,y,z,w,xb,xe,yb,ye,kxx,kyy,s,nxest,
2     * nyest,eta,tol,maxit,nmax,km1,km2,ib1,ib3,nc,intest,nrest,
3     * nx0,tx,ny0,ty,c,fp,fp0,fpint,coord,f,ff,a,q,bx,by,spx,spy,h,
4     * index,nummer,wrk,lwrk,ier)
5c  ..
6c  ..scalar arguments..
7      real*8 xb,xe,yb,ye,s,eta,tol,fp,fp0
8      integer iopt,m,kxx,kyy,nxest,nyest,maxit,nmax,km1,km2,ib1,ib3,
9     * nc,intest,nrest,nx0,ny0,lwrk,ier
10c  ..array arguments..
11      real*8 x(m),y(m),z(m),w(m),tx(nmax),ty(nmax),c(nc),fpint(intest),
12     * coord(intest),f(nc),ff(nc),a(nc,ib1),q(nc,ib3),bx(nmax,km2),
13     * by(nmax,km2),spx(m,km1),spy(m,km1),h(ib3),wrk(lwrk)
14      integer index(nrest),nummer(m)
15c  ..local scalars..
16      real*8 acc,arg,cos,dmax,fac1,fac2,fpmax,fpms,f1,f2,f3,hxi,p,pinv,
17     * piv,p1,p2,p3,sigma,sin,sq,store,wi,x0,x1,y0,y1,zi,eps,
18     * rn,one,con1,con9,con4,half,ten
19      integer i,iband,iband1,iband3,iband4,ibb,ichang,ich1,ich3,ii,
20     * in,irot,iter,i1,i2,i3,j,jrot,jxy,j1,kx,kx1,kx2,ky,ky1,ky2,l,
21     * la,lf,lh,lwest,lx,ly,l1,l2,n,ncof,nk1x,nk1y,nminx,nminy,nreg,
22     * nrint,num,num1,nx,nxe,nxx,ny,nye,nyy,n1,rank
23c  ..local arrays..
24      real*8 hx(6),hy(6)
25c  ..function references..
26      real*8 abs,fprati,sqrt
27      integer min0
28c  ..subroutine references..
29c    fpback,fpbspl,fpgivs,fpdisc,fporde,fprank,fprota
30c  ..
31c  set constants
32      one = 0.1e+01
33      con1 = 0.1e0
34      con9 = 0.9e0
35      con4 = 0.4e-01
36      half = 0.5e0
37      ten = 0.1e+02
38cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
39c part 1: determination of the number of knots and their position.     c
40c ****************************************************************     c
41c given a set of knots we compute the least-squares spline sinf(x,y),  c
42c and the corresponding weighted sum of squared residuals fp=f(p=inf). c
43c if iopt=-1  sinf(x,y) is the requested approximation.                c
44c if iopt=0 or iopt=1 we check whether we can accept the knots:        c
45c   if fp <=s we will continue with the current set of knots.          c
46c   if fp > s we will increase the number of knots and compute the     c
47c      corresponding least-squares spline until finally  fp<=s.        c
48c the initial choice of knots depends on the value of s and iopt.      c
49c   if iopt=0 we first compute the least-squares polynomial of degree  c
50c     kx in x and ky in y; nx=nminx=2*kx+2 and ny=nminy=2*ky+2.        c
51c     fp0=f(0) denotes the corresponding weighted sum of squared       c
52c     residuals                                                        c
53c   if iopt=1 we start with the knots found at the last call of the    c
54c     routine, except for the case that s>=fp0; then we can compute    c
55c     the least-squares polynomial directly.                           c
56c eventually the independent variables x and y (and the corresponding  c
57c parameters) will be switched if this can reduce the bandwidth of the c
58c system to be solved.                                                 c
59cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
60c  ichang denotes whether(1) or not(-1) the directions have been inter-
61c  changed.
62      ichang = -1
63      x0 = xb
64      x1 = xe
65      y0 = yb
66      y1 = ye
67      kx = kxx
68      ky = kyy
69      kx1 = kx+1
70      ky1 = ky+1
71      nxe = nxest
72      nye = nyest
73      eps = sqrt(eta)
74      if(iopt.lt.0) go to 20
75c  calculation of acc, the absolute tolerance for the root of f(p)=s.
76      acc = tol*s
77      if(iopt.eq.0) go to 10
78      if(fp0.gt.s) go to 20
79c  initialization for the least-squares polynomial.
80  10  nminx = 2*kx1
81      nminy = 2*ky1
82      nx = nminx
83      ny = nminy
84      ier = -2
85      go to 30
86  20  nx = nx0
87      ny = ny0
88c  main loop for the different sets of knots. m is a save upper bound
89c  for the number of trials.
90  30  do 420 iter=1,m
91c  find the position of the additional knots which are needed for the
92c  b-spline representation of s(x,y).
93        l = nx
94        do 40 i=1,kx1
95          tx(i) = x0
96          tx(l) = x1
97          l = l-1
98  40    continue
99        l = ny
100        do 50 i=1,ky1
101          ty(i) = y0
102          ty(l) = y1
103          l = l-1
104  50    continue
105c  find nrint, the total number of knot intervals and nreg, the number
106c  of panels in which the approximation domain is subdivided by the
107c  intersection of knots.
108        nxx = nx-2*kx1+1
109        nyy = ny-2*ky1+1
110        nrint = nxx+nyy
111        nreg = nxx*nyy
112c  find the bandwidth of the observation matrix a.
113c  if necessary, interchange the variables x and y, in order to obtain
114c  a minimal bandwidth.
115        iband1 = kx*(ny-ky1)+ky
116        l = ky*(nx-kx1)+kx
117        if(iband1.le.l) go to 130
118        iband1 = l
119        ichang = -ichang
120        do 60 i=1,m
121          store = x(i)
122          x(i) = y(i)
123          y(i) = store
124  60    continue
125        store = x0
126        x0 = y0
127        y0 = store
128        store = x1
129        x1 = y1
130        y1 = store
131        n = min0(nx,ny)
132        do 70 i=1,n
133          store = tx(i)
134          tx(i) = ty(i)
135          ty(i) = store
136  70    continue
137        n1 = n+1
138        if (nx.lt.ny) go to 80
139        if (nx.eq.ny) go to 120
140        go to 100
141  80    do 90 i=n1,ny
142          tx(i) = ty(i)
143  90    continue
144        go to 120
145 100    do 110 i=n1,nx
146          ty(i) = tx(i)
147 110    continue
148 120    l = nx
149        nx = ny
150        ny = l
151        l = nxe
152        nxe = nye
153        nye = l
154        l = nxx
155        nxx = nyy
156        nyy = l
157        l = kx
158        kx = ky
159        ky = l
160        kx1 = kx+1
161        ky1 = ky+1
162 130    iband = iband1+1
163c  arrange the data points according to the panel they belong to.
164        call fporde(x,y,m,kx,ky,tx,nx,ty,ny,nummer,index,nreg)
165c  find ncof, the number of b-spline coefficients.
166        nk1x = nx-kx1
167        nk1y = ny-ky1
168        ncof = nk1x*nk1y
169c  initialize the observation matrix a.
170        do 140 i=1,ncof
171          f(i) = 0.
172          do 140 j=1,iband
173            a(i,j) = 0.
174 140    continue
175c  initialize the sum of squared residuals.
176        fp = 0.
177c  fetch the data points in the new order. main loop for the
178c  different panels.
179        do 250 num=1,nreg
180c  fix certain constants for the current panel; jrot records the column
181c  number of the first non-zero element in a row of the observation
182c  matrix according to a data point of the panel.
183          num1 = num-1
184          lx = num1/nyy
185          l1 = lx+kx1
186          ly = num1-lx*nyy
187          l2 = ly+ky1
188          jrot = lx*nk1y+ly
189c  test whether there are still data points in the panel.
190          in = index(num)
191 150      if(in.eq.0) go to 250
192c  fetch a new data point.
193          wi = w(in)
194          zi = z(in)*wi
195c  evaluate for the x-direction, the (kx+1) non-zero b-splines at x(in).
196          call fpbspl(tx,nx,kx,x(in),l1,hx)
197c  evaluate for the y-direction, the (ky+1) non-zero b-splines at y(in).
198          call fpbspl(ty,ny,ky,y(in),l2,hy)
199c  store the value of these b-splines in spx and spy respectively.
200          do 160 i=1,kx1
201            spx(in,i) = hx(i)
202 160      continue
203          do 170 i=1,ky1
204            spy(in,i) = hy(i)
205 170      continue
206c  initialize the new row of observation matrix.
207          do 180 i=1,iband
208            h(i) = 0.
209 180      continue
210c  calculate the non-zero elements of the new row by making the cross
211c  products of the non-zero b-splines in x- and y-direction.
212          i1 = 0
213          do 200 i=1,kx1
214            hxi = hx(i)
215            j1 = i1
216            do 190 j=1,ky1
217              j1 = j1+1
218              h(j1) = hxi*hy(j)*wi
219 190        continue
220            i1 = i1+nk1y
221 200      continue
222c  rotate the row into triangle by givens transformations .
223          irot = jrot
224          do 220 i=1,iband
225            irot = irot+1
226            piv = h(i)
227            if(piv.eq.0.) go to 220
228c  calculate the parameters of the givens transformation.
229            call fpgivs(piv,a(irot,1),cos,sin)
230c  apply that transformation to the right hand side.
231            call fprota(cos,sin,zi,f(irot))
232            if(i.eq.iband) go to 230
233c  apply that transformation to the left hand side.
234            i2 = 1
235            i3 = i+1
236            do 210 j=i3,iband
237              i2 = i2+1
238              call fprota(cos,sin,h(j),a(irot,i2))
239 210        continue
240 220      continue
241c  add the contribution of the row to the sum of squares of residual
242c  right hand sides.
243 230      fp = fp+zi**2
244c  find the number of the next data point in the panel.
245          in = nummer(in)
246          go to 150
247 250    continue
248c  find dmax, the maximum value for the diagonal elements in the reduced
249c  triangle.
250        dmax = 0.
251        do 260 i=1,ncof
252          if(a(i,1).le.dmax) go to 260
253          dmax = a(i,1)
254 260    continue
255c  check whether the observation matrix is rank deficient.
256        sigma = eps*dmax
257        do 270 i=1,ncof
258          if(a(i,1).le.sigma) go to 280
259 270    continue
260c  backward substitution in case of full rank.
261        call fpback(a,f,ncof,iband,c,nc)
262        rank = ncof
263        do 275 i=1,ncof
264          q(i,1) = a(i,1)/dmax
265 275    continue
266        go to 300
267c  in case of rank deficiency, find the minimum norm solution.
268c  check whether there is sufficient working space
269 280    lwest = ncof*iband+ncof+iband
270        if(lwrk.lt.lwest) go to 780
271        do 290 i=1,ncof
272          ff(i) = f(i)
273          do 290 j=1,iband
274            q(i,j) = a(i,j)
275 290    continue
276        lf =1
277        lh = lf+ncof
278        la = lh+iband
279        call fprank(q,ff,ncof,iband,nc,sigma,c,sq,rank,wrk(la),
280     *    wrk(lf),wrk(lh))
281        do 295 i=1,ncof
282          q(i,1) = q(i,1)/dmax
283 295    continue
284c  add to the sum of squared residuals, the contribution of reducing
285c  the rank.
286        fp = fp+sq
287 300    if(ier.eq.(-2)) fp0 = fp
288c  test whether the least-squares spline is an acceptable solution.
289        if(iopt.lt.0) go to 820
290        fpms = fp-s
291        if(abs(fpms).le.acc) then
292          if (fp.le.0) go to 815
293          go to 820
294        endif
295c  test whether we can accept the choice of knots.
296        if(fpms.lt.0.) go to 430
297c  test whether we cannot further increase the number of knots.
298        if(ncof.gt.m) go to 790
299        ier = 0
300c  search where to add a new knot.
301c  find for each interval the sum of squared residuals fpint for the
302c  data points having the coordinate belonging to that knot interval.
303c  calculate also coord which is the same sum, weighted by the position
304c  of the data points considered.
305        do 320 i=1,nrint
306          fpint(i) = 0.
307          coord(i) = 0.
308 320    continue
309        do 360 num=1,nreg
310          num1 = num-1
311          lx = num1/nyy
312          l1 = lx+1
313          ly = num1-lx*nyy
314          l2 = ly+1+nxx
315          jrot = lx*nk1y+ly
316          in = index(num)
317 330      if(in.eq.0) go to 360
318          store = 0.
319          i1 = jrot
320          do 350 i=1,kx1
321            hxi = spx(in,i)
322            j1 = i1
323            do 340 j=1,ky1
324              j1 = j1+1
325              store = store+hxi*spy(in,j)*c(j1)
326 340        continue
327            i1 = i1+nk1y
328 350      continue
329          store = (w(in)*(z(in)-store))**2
330          fpint(l1) = fpint(l1)+store
331          coord(l1) = coord(l1)+store*x(in)
332          fpint(l2) = fpint(l2)+store
333          coord(l2) = coord(l2)+store*y(in)
334          in = nummer(in)
335          go to 330
336 360    continue
337c  find the interval for which fpint is maximal on the condition that
338c  there still can be added a knot.
339 370    l = 0
340        fpmax = 0.
341        l1 = 1
342        l2 = nrint
343        if(nx.eq.nxe) l1 = nxx+1
344        if(ny.eq.nye) l2 = nxx
345        if(l1.gt.l2) go to 810
346        do 380 i=l1,l2
347          if(fpmax.ge.fpint(i)) go to 380
348          l = i
349          fpmax = fpint(i)
350 380    continue
351c  test whether we cannot further increase the number of knots.
352        if(l.eq.0) go to 785
353c  calculate the position of the new knot.
354        arg = coord(l)/fpint(l)
355c  test in what direction the new knot is going to be added.
356        if(l.gt.nxx) go to 400
357c  addition in the x-direction.
358        jxy = l+kx1
359        fpint(l) = 0.
360        fac1 = tx(jxy)-arg
361        fac2 = arg-tx(jxy-1)
362        if(fac1.gt.(ten*fac2) .or. fac2.gt.(ten*fac1)) go to 370
363        j = nx
364        do 390 i=jxy,nx
365          tx(j+1) = tx(j)
366          j = j-1
367 390    continue
368        tx(jxy) = arg
369        nx = nx+1
370        go to 420
371c  addition in the y-direction.
372 400    jxy = l+ky1-nxx
373        fpint(l) = 0.
374        fac1 = ty(jxy)-arg
375        fac2 = arg-ty(jxy-1)
376        if(fac1.gt.(ten*fac2) .or. fac2.gt.(ten*fac1)) go to 370
377        j = ny
378        do 410 i=jxy,ny
379          ty(j+1) = ty(j)
380          j = j-1
381 410    continue
382        ty(jxy) = arg
383        ny = ny+1
384c  restart the computations with the new set of knots.
385 420  continue
386c  test whether the least-squares polynomial is a solution of our
387c  approximation problem.
388 430  if(ier.eq.(-2)) go to 830
389cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
390c part 2: determination of the smoothing spline sp(x,y)                c
391c *****************************************************                c
392c we have determined the number of knots and their position. we now    c
393c compute the b-spline coefficients of the smoothing spline sp(x,y).   c
394c the observation matrix a is extended by the rows of a matrix,        c
395c expressing that sp(x,y) must be a polynomial of degree kx in x and   c
396c ky in y. the corresponding weights of these additional rows are set  c
397c to 1./p.  iteratively we than have to determine the value of p       c
398c such that f(p)=sum((w(i)*(z(i)-sp(x(i),y(i))))**2) be = s.           c
399c we already know that the least-squares polynomial corresponds to     c
400c p=0  and that the least-squares spline corresponds to p=infinity.    c
401c the iteration process which is proposed here makes use of rational   c
402c interpolation. since f(p) is a convex and strictly decreasing        c
403c function of p, it can be approximated by a rational function r(p)=   c
404c (u*p+v)/(p+w). three values of p(p1,p2,p3) with corresponding values c
405c of f(p) (f1=f(p1)-s,f2=f(p2)-s,f3=f(p3)-s) are used to calculate the c
406c new value of p such that r(p)=s. convergence is guaranteed by taking c
407c f1 > 0 and f3 < 0.                                                   c
408cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
409      kx2 = kx1+1
410c  test whether there are interior knots in the x-direction.
411      if(nk1x.eq.kx1) go to 440
412c  evaluate the discotinuity jumps of the kx-th order derivative of
413c  the b-splines at the knots tx(l),l=kx+2,...,nx-kx-1.
414      call fpdisc(tx,nx,kx2,bx,nmax)
415 440  ky2 = ky1 + 1
416c  test whether there are interior knots in the y-direction.
417      if(nk1y.eq.ky1) go to 450
418c  evaluate the discontinuity jumps of the ky-th order derivative of
419c  the b-splines at the knots ty(l),l=ky+2,...,ny-ky-1.
420      call fpdisc(ty,ny,ky2,by,nmax)
421c  initial value for p.
422 450  p1 = 0.
423      f1 = fp0-s
424      p3 = -one
425      f3 = fpms
426      p = 0.
427      do 460 i=1,ncof
428        p = p+a(i,1)
429 460  continue
430      rn = ncof
431      p = rn/p
432c  find the bandwidth of the extended observation matrix.
433      iband3 = kx1*nk1y
434      iband4 = iband3 +1
435      ich1 = 0
436      ich3 = 0
437c  iteration process to find the root of f(p)=s.
438      do 770 iter=1,maxit
439        pinv = one/p
440c  store the triangularized observation matrix into q.
441        do 480 i=1,ncof
442          ff(i) = f(i)
443          do 470 j=1,iband
444            q(i,j) = a(i,j)
445 470      continue
446          ibb = iband+1
447          do 480 j=ibb,iband4
448            q(i,j) = 0.
449 480    continue
450        if(nk1y.eq.ky1) go to 560
451c  extend the observation matrix with the rows of a matrix, expressing
452c  that for x=cst. sp(x,y) must be a polynomial in y of degree ky.
453        do 550 i=ky2,nk1y
454          ii = i-ky1
455          do 550 j=1,nk1x
456c  initialize the new row.
457            do 490 l=1,iband
458              h(l) = 0.
459 490        continue
460c  fill in the non-zero elements of the row. jrot records the column
461c  number of the first non-zero element in the row.
462            do 500 l=1,ky2
463              h(l) = by(ii,l)*pinv
464 500        continue
465            zi = 0.
466            jrot = (j-1)*nk1y+ii
467c  rotate the new row into triangle by givens transformations without
468c  square roots.
469            do 540 irot=jrot,ncof
470              piv = h(1)
471              i2 = min0(iband1,ncof-irot)
472              if(piv.eq.0.) then
473                if (i2.le.0) go to 550
474                go to 520
475              endif
476c  calculate the parameters of the givens transformation.
477              call fpgivs(piv,q(irot,1),cos,sin)
478c  apply that givens transformation to the right hand side.
479              call fprota(cos,sin,zi,ff(irot))
480              if(i2.eq.0) go to 550
481c  apply that givens transformation to the left hand side.
482              do 510 l=1,i2
483                l1 = l+1
484                call fprota(cos,sin,h(l1),q(irot,l1))
485 510          continue
486 520          do 530 l=1,i2
487                h(l) = h(l+1)
488 530          continue
489              h(i2+1) = 0.
490 540        continue
491 550    continue
492 560    if(nk1x.eq.kx1) go to 640
493c  extend the observation matrix with the rows of a matrix expressing
494c  that for y=cst. sp(x,y) must be a polynomial in x of degree kx.
495        do 630 i=kx2,nk1x
496          ii = i-kx1
497          do 630 j=1,nk1y
498c  initialize the new row
499            do 570 l=1,iband4
500              h(l) = 0.
501 570        continue
502c  fill in the non-zero elements of the row. jrot records the column
503c  number of the first non-zero element in the row.
504            j1 = 1
505            do 580 l=1,kx2
506              h(j1) = bx(ii,l)*pinv
507              j1 = j1+nk1y
508 580        continue
509            zi = 0.
510            jrot = (i-kx2)*nk1y+j
511c  rotate the new row into triangle by givens transformations .
512            do 620 irot=jrot,ncof
513              piv = h(1)
514              i2 = min0(iband3,ncof-irot)
515              if(piv.eq.0.) then
516                if (i2.le.0) go to 630
517                go to 600
518              endif
519c  calculate the parameters of the givens transformation.
520              call fpgivs(piv,q(irot,1),cos,sin)
521c  apply that givens transformation to the right hand side.
522              call fprota(cos,sin,zi,ff(irot))
523              if(i2.eq.0) go to 630
524c  apply that givens transformation to the left hand side.
525              do 590 l=1,i2
526                l1 = l+1
527                call fprota(cos,sin,h(l1),q(irot,l1))
528 590          continue
529 600          do 610 l=1,i2
530                h(l) = h(l+1)
531 610          continue
532              h(i2+1) = 0.
533 620        continue
534 630    continue
535c  find dmax, the maximum value for the diagonal elements in the
536c  reduced triangle.
537 640    dmax = 0.
538        do 650 i=1,ncof
539          if(q(i,1).le.dmax) go to 650
540          dmax = q(i,1)
541 650    continue
542c  check whether the matrix is rank deficient.
543        sigma = eps*dmax
544        do 660 i=1,ncof
545          if(q(i,1).le.sigma) go to 670
546 660    continue
547c  backward substitution in case of full rank.
548        call fpback(q,ff,ncof,iband4,c,nc)
549        rank = ncof
550        go to 675
551c  in case of rank deficiency, find the minimum norm solution.
552 670    lwest = ncof*iband4+ncof+iband4
553        if(lwrk.lt.lwest) go to 780
554        lf = 1
555        lh = lf+ncof
556        la = lh+iband4
557        call fprank(q,ff,ncof,iband4,nc,sigma,c,sq,rank,wrk(la),
558     *   wrk(lf),wrk(lh))
559 675    do 680 i=1,ncof
560          q(i,1) = q(i,1)/dmax
561 680    continue
562c  compute f(p).
563        fp = 0.
564        do 720 num = 1,nreg
565          num1 = num-1
566          lx = num1/nyy
567          ly = num1-lx*nyy
568          jrot = lx*nk1y+ly
569          in = index(num)
570 690      if(in.eq.0) go to 720
571          store = 0.
572          i1 = jrot
573          do 710 i=1,kx1
574            hxi = spx(in,i)
575            j1 = i1
576            do 700 j=1,ky1
577              j1 = j1+1
578              store = store+hxi*spy(in,j)*c(j1)
579 700        continue
580            i1 = i1+nk1y
581 710      continue
582          fp = fp+(w(in)*(z(in)-store))**2
583          in = nummer(in)
584          go to 690
585 720    continue
586c  test whether the approximation sp(x,y) is an acceptable solution.
587        fpms = fp-s
588        if(abs(fpms).le.acc) go to 820
589c  test whether the maximum allowable number of iterations has been
590c  reached.
591        if(iter.eq.maxit) go to 795
592c  carry out one more step of the iteration process.
593        p2 = p
594        f2 = fpms
595        if(ich3.ne.0) go to 740
596        if((f2-f3).gt.acc) go to 730
597c  our initial choice of p is too large.
598        p3 = p2
599        f3 = f2
600        p = p*con4
601        if(p.le.p1) p = p1*con9 + p2*con1
602        go to 770
603 730    if(f2.lt.0.) ich3 = 1
604 740    if(ich1.ne.0) go to 760
605        if((f1-f2).gt.acc) go to 750
606c  our initial choice of p is too small
607        p1 = p2
608        f1 = f2
609        p = p/con4
610        if(p3.lt.0.) go to 770
611        if(p.ge.p3) p = p2*con1 + p3*con9
612        go to 770
613 750    if(f2.gt.0.) ich1 = 1
614c  test whether the iteration process proceeds as theoretically
615c  expected.
616 760    if(f2.ge.f1 .or. f2.le.f3) go to 800
617c  find the new value of p.
618        p = fprati(p1,f1,p2,f2,p3,f3)
619 770  continue
620c  error codes and messages.
621 780  ier = lwest
622      go to 830
623 785  ier = 5
624      go to 830
625 790  ier = 4
626      go to 830
627 795  ier = 3
628      go to 830
629 800  ier = 2
630      go to 830
631 810  ier = 1
632      go to 830
633 815  ier = -1
634      fp = 0.
635 820  if(ncof.ne.rank) ier = -rank
636c  test whether x and y are in the original order.
637 830  if(ichang.lt.0) go to 930
638c  if not, interchange x and y once more.
639      l1 = 1
640      do 840 i=1,nk1x
641        l2 = i
642        do 840 j=1,nk1y
643          f(l2) = c(l1)
644          l1 = l1+1
645          l2 = l2+nk1x
646 840  continue
647      do 850 i=1,ncof
648        c(i) = f(i)
649 850  continue
650      do 860 i=1,m
651        store = x(i)
652        x(i) = y(i)
653        y(i) = store
654 860  continue
655      n = min0(nx,ny)
656      do 870 i=1,n
657        store = tx(i)
658        tx(i) = ty(i)
659        ty(i) = store
660 870  continue
661      n1 = n+1
662      if (nx.lt.ny) go to 880
663      if (nx.eq.ny) go to 920
664      go to 900
665 880  do 890 i=n1,ny
666        tx(i) = ty(i)
667 890  continue
668      go to 920
669 900  do 910 i=n1,nx
670        ty(i) = tx(i)
671 910  continue
672 920  l = nx
673      nx = ny
674      ny = l
675 930  if(iopt.lt.0) go to 940
676      nx0 = nx
677      ny0 = ny
678 940  return
679      end
680
681