1      subroutine tn_fitting_coeff(n,m,q,f,c,r)
2*
3* $Id$
4*
5      implicit none
6#include "errquit.fh"
7      integer n, m
8      double precision q(0:n,1:m), f(1:m), c(0:n), r
9c
10c     Given the matrix from tn_fitting_matrix return
11c     the fitting coeffs for a given function f() and
12c     the residual at the fitting points
13c
14      integer i, j
15      integer maxn
16      parameter (maxn=64)
17      double precision t(0:maxn), x, sum
18c
19      if (n .gt. maxn) call errquit('tn_f_c: n>maxn ', n, UNKNOWN_ERR)
20c
21      do i = 0,n
22         c(i) = 0.0d0
23      enddo
24c
25      do j = 1,m
26         do i = 0,n
27            c(i) = c(i) + q(i,j)*f(j)
28         enddo
29      enddo
30c
31      r = 0.0d0
32      do j = 1, m
33         x = -1.0d0 + 2.0d0*(j-1)/dble(m-1)
34         call tn_evaluate(n,x,t)
35         sum = 0.0d0
36         do i = 0,n
37            sum = sum + c(i)*t(i)
38         enddo
39         r = r + (f(j)-sum)**2
40      enddo
41c
42      r = sqrt(r)
43c
44      end
45      subroutine tn_collate_matrix(nx,xlo,xhi,nptx,x,order,q,lq)
46      implicit none
47#include "errquit.fh"
48      integer nx, nptx, lq, order
49      double precision xlo, xhi, x(nptx), q(lq,nx)
50c
51c     Return in q(nx,nptx) the matrix that uses Chebyshev interpolation
52c     of the requested order to interpolate from nx points uniformly
53c     distributed in [xlo,xhi] to the list of points x(i), i=1,nptx
54c
55      integer maxorder, maxn
56      parameter (maxorder=35,maxn=64)
57      double precision qq((maxorder+1)*maxn),t(0:maxorder)
58      double precision xx
59      integer i,j,k,ind
60c
61      if (nx .gt. maxn) call errquit('tn_collate: nx>maxn',nx,
62     &       UNKNOWN_ERR)
63      if (order.gt.maxorder) call errquit('tn_collate: order',order,
64     &       UNKNOWN_ERR)
65c
66      call tn_fitting_matrix(order,nx,qq)
67      do j = 1, nx
68         do i = 1, nptx
69            q(i,j) = 0.0d0
70         enddo
71      enddo
72      do i = 1, nptx
73         xx = 2.0d0*(x(i)-xlo)/(xhi-xlo) - 1.0d0
74         call tn_evaluate(order,xx,t)
75         ind = 1
76         do j = 1, nx
77            do k = 0,order
78               q(i,j) = q(i,j) + t(k)*qq(ind+k)
79            enddo
80            ind = ind + order + 1
81         enddo
82      enddo
83c
84      end
85      subroutine tn_fitting_matrix(n,m,q)
86      implicit none
87#include "errquit.fh"
88      integer n,m
89      double precision q(0:n,1:m)
90c
91c     Return in q the matrix that forms a least-squares fit of
92c     m points uniformly distributed on [-1,1] to a Chebyshev
93c     expansion of order n.
94c
95c     Q = (At.A)^-1.A where Aij = Tj(xi), xi = -1 + 2*(i-1)/(m-1)
96c
97c     Then, given a function f() evaluated at the m uniform points
98c     a fit is computed as
99c
100c     ci = sum(j) Qij.f(xj) , or c = Qf
101c
102c     f(x) ~= sum(i) Ti(x)*ci
103c
104c     This calls the LAPACK routine dgels to use QR or a related
105c     method ... DON'T use the normal equations ... I tried it
106c     and it is too unstable for the present purpose.
107c
108      integer maxn, maxm, lwork
109      parameter (maxn=64,maxm=100,lwork=10000)
110      double precision at((maxn+1)*maxm),
111     $     b((maxm+1)*(maxm+1)),h, x, work(lwork)
112      integer i, j, info
113c
114      if (n.gt.maxn .or. m.gt.maxm) call errquit('tn_f_m: n or m ',
115     &       UNKNOWN_ERR,
116     $     n*10000+m)
117c
118      h = 2.0d0/(m-1)
119      do i = 1, m
120         x = -1.0d0 + h*(i-1)
121         call tn_evaluate(n,x,at(1+(i-1)*(n+1)))
122      enddo
123c
124      call dfill(m*m,0.0d0,b,1)
125      call dfill(m,1.0d0,b,m+1)
126      call dgels('t',n+1,m,m,at,n+1,b,m,work,lwork,info)
127      if (info .ne. 0) call errquit('tn_f_m: dgels failed ', info,
128     &       UNKNOWN_ERR)
129c
130      do i = 1, m
131         do j = 0,n
132            q(j,i) = b(1 + j + (i-1)*m)
133         enddo
134      enddo
135c
136      end
137      subroutine tn_evaluate(n,x,t)
138      implicit none
139      double precision x
140      integer n
141      double precision t(0:n)
142c
143c     Compute the Chebyshev polynomials defined on -1,1 at the
144c     point x.  n must be at least 1.
145c
146      integer j
147c
148      t(0)=1.0d0
149      t(1)=x
150      if (n.gt.1) then
151         do j=1,n-1
152            t(j+1) = 2.0d0*x*t(j) - t(j-1)
153         enddo
154      endif
155c
156      end
157      subroutine tn_interp_2d(nx, ny, xlo, xhi, ylo, yhi, nptx, npty,
158     $     x, y, f, ldx, ff, ldxff, order)
159      implicit none
160#include "errquit.fh"
161c
162c     Given a discretization of the area (xlo:xhi,ylo:yhi)
163c     in f(1:nx,1:ny) return in ff(1:nptx,1:npty)
164c     the values of f() interpolated onto the coordinates
165c     (x(i),y(j)) i=1,nptx, j=1,npty using a least-squares
166c     Chebyshev approximation of the given order in each dimension.
167c
168c     The input array f() is DESTROYED and must be dimensioned
169c     ldx>=max(nx,ptx)
170c
171      integer nx, ny, nptx, npty, ldx, ldxff, order
172      double precision xlo,xhi,ylo,yhi
173      double precision f(ldx, *), ff(ldxff,*)
174      double precision x(nptx), y(npty)
175c
176      integer i, j, l
177      integer maxnpt, maxn
178      parameter (maxnpt=100, maxn=100)
179      double precision c(maxnpt,maxn), fij, tmp(maxnpt)
180c
181      if (nx.gt.ldx) call errquit('tn_i_2d: nx>ldx',nx*10000+ldx,
182     &       UNKNOWN_ERR)
183      if (nptx.gt.ldx) call errquit('tn_i_2d: nptx>ldx',nptx*10000+ldx,
184     &       UNKNOWN_ERR)
185      if (nx.gt.maxn .or. nptx.gt.maxnpt) call errquit
186     $     ('tn_i_2d: nx or nptx', 10000*nx+nptx,
187     &       UNKNOWN_ERR)
188      if (ny.gt.maxn .or. npty.gt.maxnpt) call errquit
189     $     ('tn_i_2d: ny or npty', 10000*ny+npty,
190     &       UNKNOWN_ERR)
191c
192      do l = 1, nptx
193         if (x(l).lt.xlo .or. x(l).gt.xhi) call errquit
194     $        ('tn_interp_2d: x is out of range',l,
195     &       UNKNOWN_ERR)
196      enddo
197      do l = 1, npty
198         if (y(l).lt.ylo .or. y(l).gt.yhi) call errquit
199     $        ('tn_interp_2d: y is out of range',l,
200     &       UNKNOWN_ERR)
201      enddo
202c
203      call tn_collate_matrix(nx, xlo, xhi, nptx, x, order, c, maxnpt)
204      do j = 1, ny
205         do l = 1, nptx
206            tmp(l) = 0.0d0
207         end do
208         do i = 1, nx
209            fij = f(i,j)
210            if (abs(fij) .gt. 0.0d0) then
211               do l = 1, nptx
212                  tmp(l) = tmp(l)+c(l,i)*fij
213               end do
214            end if
215         end do
216         do l = 1, nptx
217            f(l,j) = tmp(l)
218         end do
219      end do
220c
221      call tn_collate_matrix(ny, ylo, yhi, npty, y, order, c, maxnpt)
222      do i = 1, nptx
223         do l = 1, npty
224            tmp(l) = 0.0d0
225         end do
226         do j = 1, ny
227            fij = f(i,j)
228            if (abs(fij) .gt. 0.0d0) then
229               do l = 1, npty
230                  tmp(l) = tmp(l) + c(l,j)*fij
231               end do
232            end if
233         end do
234         do l = 1, npty
235            ff(i,l) = tmp(l)
236         end do
237      end do
238c
239      end
240      double precision function tn_interp_3d_point(
241     $     g, nx, ny, nz, hx, hy, hz,
242     $     x, y, z,
243     $     n, order, qq)
244      implicit none
245#include "errquit.fh"
246c
247c     Given a function tablulated on a grid with given spacing return
248c     the value interpolated at (x,y,z) within the cube
249c     (0:hx*(nx-1),0:hy*(ny-1),0:hz*(nz-1))
250c     a Chebyshev fit of the given order to the given number of points
251c     in all dimensions.
252c
253c     Note that the exterior points on the grid are ommited so
254c     the grid is numbered rather oddly.  See solver.F.
255c
256c     The input array of points is PRESERVED.
257c
258      integer nx, ny, nz, n, order
259      double precision x, y, z, g(nx,ny,nz), hx, hy, hz, qq(*)
260c
261      integer i, j, k, ind, ilo, jlo, klo, nhalf, jtop
262      integer maxn, maxorder
263      parameter (maxn=63, maxorder=53)
264      double precision cx(maxn), cy(maxn), cz(maxn), sumi, sumj, sum,
265     $     sumi0, sumi1, sumi2, sumi3
266      double precision xx, yy, zz, twornm1,
267     $     tx(0:maxorder), ty(0:maxorder), tz(0:maxorder)
268c
269      if (n.gt.maxn) call errquit('tn_i_3d_pt: n', n, UNKNOWN_ERR)
270      if (order.gt.maxorder) call errquit('tn_i_3d_pt:order',order,
271     &       UNKNOWN_ERR)
272c
273      nhalf = n/2
274      twornm1 = 2.0d0/dble(n-1)
275      xx = x/hx                 ! Coords within whole volume
276      yy = y/hy
277      zz = z/hz
278      ilo = max(1, int(xx)-nhalf) ! Corner of interpolation volume
279      ilo = min(ilo,nx-n+1)
280      jlo = max(1, int(yy)-nhalf)
281      jlo = min(jlo,ny-n+1)
282      klo = max(1, int(zz)-nhalf)
283      klo = min(klo,nz-n+1)
284      xx = (xx-dble(ilo))*twornm1 - 1.0d0 ! Rescale to [-1,1]
285      yy = (yy-dble(jlo))*twornm1 - 1.0d0
286      zz = (zz-dble(klo))*twornm1 - 1.0d0
287c
288      tx(0)=1.0d0
289      ty(0)=1.0d0
290      tz(0)=1.0d0
291      tx(1)=xx
292      ty(1)=yy
293      tz(1)=zz
294      do j=1,order-1
295         tx(j+1) = 2.0d0*xx*tx(j) - tx(j-1)
296         ty(j+1) = 2.0d0*yy*ty(j) - ty(j-1)
297         tz(j+1) = 2.0d0*zz*tz(j) - tz(j-1)
298      enddo
299c
300      ind = 1
301      do j = 1, n
302         cx(j) = 0.0d0
303         cy(j) = 0.0d0
304         cz(j) = 0.0d0
305         do k = 0,order
306            cx(j) = cx(j) + tx(k)*qq(ind+k)
307            cy(j) = cy(j) + ty(k)*qq(ind+k)
308            cz(j) = cz(j) + tz(k)*qq(ind+k)
309         enddo
310         ind = ind + order + 1
311      enddo
312c
313      ilo = ilo - 1
314      jlo = jlo - 1
315      klo = klo - 1
316      sum = 0.0d0
317c
318c     Manually unroll the j loop into the i loop
319c
320      jtop = (n/4)*4
321      do k = 1, n
322         sumj = 0.0d0
323         do j = 1, jtop, 4
324            sumi0 = 0.0d0
325            sumi1 = 0.0d0
326            sumi2 = 0.0d0
327            sumi3 = 0.0d0
328            do i = 1, n
329               sumi0 = sumi0 + cx(i)*g(i+ilo,j+jlo  ,k+klo)
330               sumi1 = sumi1 + cx(i)*g(i+ilo,j+jlo+1,k+klo)
331               sumi2 = sumi2 + cx(i)*g(i+ilo,j+jlo+2,k+klo)
332               sumi3 = sumi3 + cx(i)*g(i+ilo,j+jlo+3,k+klo)
333            end do
334            sumj = sumj + cy(j)*sumi0 + cy(j+1)*sumi1 + cy(j+2)*sumi2 +
335     $           cy(j+3)*sumi3
336         end do
337         do j = jtop+1,n
338            sumi = 0.0d0
339            do i = 1, n
340               sumi = sumi + cx(i)*g(i+ilo,j+jlo,k+klo)
341            end do
342            sumj = sumj + cy(j)*sumi
343         end do
344         sum = sum + cz(k)*sumj
345      end do
346c$$$      do k = 1, n
347c$$$         sumj = 0.0d0
348c$$$         do j = 1, n
349c$$$            sumi = 0.0d0
350c$$$            do i = 1, n
351c$$$               sumi = sumi + cx(i)*g(i+ilo,j+jlo,k+klo)
352c$$$            end do
353c$$$            sumj = sumj + cy(j)*sumi
354c$$$         end do
355c$$$         sum = sum + cz(k)*sumj
356c$$$      end do
357c
358      tn_interp_3d_point = sum
359c
360      end
361      subroutine tn_lsq_fit(m,nfunc,order,f)
362      implicit none
363#include "errquit.fh"
364      integer m, nfunc, order
365      double precision f(m,nfunc)
366c
367c     Input is a set of nfunc functions on a uniform grid (1:m)
368c     tablulated in f(1:m,1:nfunc).
369c
370c     Return in f(1:order+1,1:nfunc) a least-squares fit to
371c     Chebyshev polynomials up to the specified order.
372c
373c     order+1 <= m
374c
375      integer maxorder, maxm, lwork
376      parameter (maxorder=35,maxm=63,lwork=10000)
377      double precision at((maxorder+1)*maxm), h, x, work(lwork)
378      integer i, info
379c
380      if (order+1.gt.m) call errquit('order+1>m',order,
381     &       UNKNOWN_ERR)
382      if (order.gt.maxorder) call errquit('order ', order,
383     &       UNKNOWN_ERR)
384      if (m .gt. maxm) call errquit('m', m,
385     &       UNKNOWN_ERR)
386c
387      h = 2.0d0/(m-1)
388      do i = 1, m
389         x = -1.0d0 + h*dble(i-1)
390         call tn_evaluate(order,x,at(1+(i-1)*(order+1)))
391      enddo
392c
393      call dgels('t',order+1,m,nfunc,at,order+1,f,m,work,lwork,info)
394      if (info .ne. 0) call errquit('info ', info,
395     &       UNKNOWN_ERR)
396c
397*      write(6,*) ' solution from dgels'
398*      call doutput(f,1,order+1,1,nfunc,order+1,nfunc,1)
399c
400      end
401      subroutine tn_lsq_fit_cube(m,order,f)
402      implicit none
403      integer m, order
404      double precision f(m,m,m)
405c
406c     Input in f(1:m,1:m,1:m) is a function tabulated in a cube.
407c
408c     Return in f(1:order+1,1:order+1,1:order+1) a least squares
409c     fit to Chebychev polynomials of the given order.
410c     Also return in r the residual.
411c
412      integer maxorder, maxm
413      parameter (maxorder=33, maxm=63)
414      double precision g(maxm*maxm*maxm)
415      integer i, j, k, op1, ind
416c
417      op1 = order+1
418c
419c     Transform the first dimension
420c
421      call tn_lsq_fit(m,m*m,order,f)
422c
423c     Rotate the dimensions and do the next
424c
425      call rotate_dims(f,m,m,g,m,op1,op1,m,m)
426      call tn_lsq_fit(m,(op1)*m,order,g)
427c
428c     And the last one
429c
430      call rotate_dims(g,m,op1,f,m,op1,op1,op1,m)
431      call tn_lsq_fit(m,(op1)**2,order,f)
432c
433c     Rotate again and copy back
434c
435      call rotate_dims(f,m,op1,g,op1,op1,op1,op1,op1)
436      ind = 1
437      do k = 1, op1
438         do j = 1, op1
439            do i = 1, op1
440               f(i,j,k) = g(ind)
441               ind = ind + 1
442            enddo
443         enddo
444      enddo
445c
446      end
447      subroutine rotate_dims(f, ldf1, ldf2, g, ldg1, ldg2,
448     $     ni, nj, nk)
449      implicit none
450      integer ldf1, ldf2, ldg1, ldg2, ni, nj, nk
451      double precision f(ldf1,ldf2,*), g(ldg1,ldg2,*)
452c
453c     Copy f into g permuting the dimensions so that g(k,i,j) = f(i,j,k)
454c
455      integer i, j, k
456c
457      do k = 1, nk
458         do j = 1, nj
459            do i = 1, ni
460               g(k,i,j) = f(i,j,k)
461            enddo
462         enddo
463      enddo
464c
465      end
466      double precision function tn_cube_eval(g,ldg1,ldg2,order,x,y,z)
467      implicit none
468      integer ldg1, ldg2, order
469      double precision g(ldg1,ldg2,*), x, y, z
470c
471c     Given in g a Chebychev fit of given order of a function defined
472c     in a cube [-1,1] compute the value at the interior point x,y,z
473c
474      integer maxorder
475      parameter (maxorder = 33)
476      double precision tx(maxorder+1), ty(maxorder+1), tz(maxorder+1)
477      double precision sum, sumi, sumj
478      integer i, j, k, op1
479c
480      tx(1)=1.0d0
481      ty(1)=1.0d0
482      tz(1)=1.0d0
483      tx(2)=x
484      ty(2)=y
485      tz(2)=z
486      do j=2,order
487         tx(j+1) = 2.0d0*x*tx(j) - tx(j-1)
488         ty(j+1) = 2.0d0*y*ty(j) - ty(j-1)
489         tz(j+1) = 2.0d0*z*tz(j) - tz(j-1)
490      enddo
491*      call tn_evaluate(order,x,tx)
492*      call tn_evaluate(order,y,ty)
493*      call tn_evaluate(order,z,tz)
494c
495      op1 = order+1
496      sum = 0.0d0
497      do k = 1, op1
498         sumj = 0.0d0
499         do j = 1, op1
500            sumi = 0.0d0
501            do i = 1, op1
502               sumi = sumi + tx(i)*g(i,j,k)
503            end do
504            sumj = sumj + ty(j)*sumi
505         end do
506         sum = sum + tz(k)*sumj
507      end do
508c
509      tn_cube_eval = sum
510c
511      end
512