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