program lor990 c c Exemple of construction of recurrence coefficients c for a Lorentzian weight function c c w(x) = x^p (1+x^2/a)^(-k-1) on (0,c) c c c In the paper, p=2, a=kappa, b=kappa+1 parameter(ndim=2000) integer k,nc ,j,jr double precision ka,aux,aux1,aux2,sk,c,a,b,rr,ri,rro,rio,al2,ep dimension a(ndim),b(ndim),rr(ndim),ri(ndim),rro(ndim),rio(ndim) complex z,sz c read p a k c precision number of coeff. c integer c open(unit=1,file='loren99.d',status='old') read(1,*)al2,ka,k,c,ep,nc print '('' x** '',f5.2,''(1+x**2/'',f6.3,'')**('',i6,'') '')', & al2,ka,-k-1 print '('' on (0,'',f7.2,'')'')',c print '(i5,'' recurrence coefficients'')',nc c call lorentz(al2,ka,k,c,ep,nc,ndim,a,b,rr,ri,rro,rio,ier) print *,' ier subroutine= ',ier,' (ndim=',ndim c do 4 j=1,nc+1 4 print '(i5,1x,2f22.15)',j-1,a(j),b(j) c c Test reconstruction of weight function at x=1 with c Haydock & Nex' continued fraction z=1.0 c Termination of continued fraction b/(z-a-b/(z-a- ...)) sz=(2*z-c)/c sz=sz+sqrt(sz*sz-1) if(abs(sz).gt.1)sz=1/sz sz=c*sz/4 aux=sz aux1=imag(sz) do 43 jr1=1,nc jr=nc+1-jr1 c 23 rho(jr)=b(jr)/(z-a(jr) -rho(jr+1) ) aux2=(1-a(jr)-aux)**2+aux1**2 aux=b(jr)*(1-a(jr)-aux)/aux2 43 aux1=b(jr)*aux1/aux2 print *,' check at x=1:' print '(i5,1x,2f22.15)',nc,1/(1.0d0+1.0d0/ka)**(k+1), & aux1/3.14159265358979d0 c c nodes and weights c do 6 j=1,nc-1 rr(j)=a(j) rro(j)=sqrt( b(j+1) ) 6 ri(j+1)=0 ri(1)=1 rr(nc)=a(nc) call gausq2(nc, rr, rro, ri, ep,ierr) if(ierr.ne.0)print *,' err gaussq2=',ierr do 65 j=1,nc ri(j)=b(1)*ri(j)*ri(j) 65 print '(2f20.12)',rr(j),ri(j) c check (1+x2/a)**(k+1) aux=0 do 67 j=1,nc 67 aux=aux+ri(j)* (1+rr(j)**2/ka)**(k+1) print *,' check quadrature ',aux, c**(al2+1)/(al2+1) stop end c subroutine lorentz(al2,ka,k,c,ep,nc,ndim,a,b,rr,ri,rro,rio,ier) c Construction of recurrence coefficients c for a Lorentzian weight function c w(x) = x^p (1+x^2/a)^(-k-1) on (0,c) c c monic pol.: Q_{n+1}(x) = (x-\alpha_n)Q_n(x) - \beta_n Q_{n-1}(x) c c input: c al2 = p scalar real > -1 c ka = a scalar real > 0 c k scalar integer >= 0 c c scalar real > 0 c ep = precision scalar real c nc = number of wanted coefficients scalar integer > 0 c ndim = available vector dimension c c output: c a vector real alpha : alpha(0) = a(1), etc. c b vector real beta : beta(1) = b(2), etc. c b(1) = beta(0) = total weight c c ier scalar integer = largest index encountered in c vectors. Caution if close to ndim! c c rr, ri, rro,rio : auxiliary vectors c integer k,nc,nx0,nx ,nn,i,j,jr,jr1 double precision ka,aux,aux1,aux2,sk,c,a,b,rr,ri,rro,rio,al2,ep dimension a(ndim),b(ndim),rr(ndim),ri(ndim),rro(ndim),rio(ndim) complex z,sz sk=sqrt(ka) ssk=sk z=cmplx(0.0, ssk ) c print *,' z= ',z c Termination of continued fraction b/(z-a-b/(z-a- ...)) z=(2*z-c)/c sz=z+sqrt(z*z-1) if(abs(sz).gt.1)sz=1/sz nx0=1+0.5*log(ep)/log( abs(sz) ) sz=c*sz/4 c print *,' terminator = ',sz, ' overshot per kappa=',nx0 nx=nc+(k+1)*nx0 1 if(nx.gt.ndim-1)nx=ndim-1 ier=nx c coeff. de x^al2 sur (0,c) : a(1)=c*(al2+1)/(al2+2) do 21 i=1,nx a(i+1)=c*( i*(i+al2+1)+al2*(al2+1)/2 )/ & ( 2*(i+al2/2)*(i+1+al2/2) ) 21 b(i+1)=c*c*i*i*(i+al2)*(i+al2)/ & ((2*i+al2-1)*(2*i+al2)*(2*i+al2)*(2*i+al2+1)) b(1)=c**(al2+1) /(al2+1) c k+1 divisions of the weight by 1+x^2/k : nn=nx-1 do 2 j=1,k+1 c print '('' division '',i3,i5,1p,2e20.12)',j,nn,a(1),b(1) c print '(18x, 1p,2e20.12)', a(2),b(2) c print '(18x, 1p,2e20.12)', a(3),b(3) rr(nn+1)=sz ri(nn+1)=imag(sz) rro(nn-4)=rr(nn+1) rio(nn-4)=ri(nn+1) do 23 jr1=1,nn jr=nn+1-jr1 c 23 rho(jr)=b(jr)/(z-a(jr) -rho(jr+1) ) aux=(a(jr)+rr(jr+1))**2+(sk-ri(jr+1))**2 rr(jr)=-b(jr)*(a(jr)+rr(jr+1))/aux 23 ri(jr)=b(jr)*(ri(jr+1)-sk)/aux nn=nn-5 do 231 jr1=1,nn jr=nn+1-jr1 aux=(a(jr)+rro(jr+1))**2+(sk-rio(jr+1))**2 rro(jr)=-b(jr)*(a(jr)+rro(jr+1))/aux 231 rio(jr)=b(jr)*(rio(jr+1)-sk)/aux do 232 jr=1,nn nn2=jr ab=abs(rr(jr))+abs(ri(jr)) if(abs(rro(jr)-rr(jr))+abs(rio(jr)-ri(jr)).gt.ep*ab)goto 235 232 continue 235 nn=nn2 aux=a(1) a(1)=sk*rr(1)/ri(1) b(1)=-ri(1)*sk do 24 jr=2,nn aux1=a(jr) a(jr)=aux+sk*( rr(jr)/ri(jr)-rr(jr-1)/ri(jr-1) ) aux=aux1 24 b(jr)=-(-ri(jr-1)+sk )*(rr(jr-1)**2+ri(jr-1)**2)*ri(jr) & /ri(jr-1)**2 b(2)=-sk*(rr(1)**2+ri(1)**2)*ri(2)/ri(1)**2 if(nn.lt.nc)then c print *,' reprise avec un index max plus grand ' if(nx.eq.ndim-1)return nx=nx+20 goto 1 endif 2 continue return end c c subroutine gausq2(n, d, e, z, ep,ierr) c c this subroutine is a translation of an algol procedure, c num. math. 12, 377-383(1968) by martin and wilkinson, c as modified in num. math. 15, 450(1970) by dubrulle. c handbook for auto. comp., vol.ii-linear algebra, 241-248(1971). c this is a modified version of the 'eispack' routine imtql2. c c this subroutine finds the eigenvalues and first components of the c eigenvectors of a symmetric tridiagonal matrix by the implicit ql c method. c c on input: c c n is the order of the matrix; c c d contains the diagonal elements of the input matrix; c c e contains the subdiagonal elements of the input matrix c in its first n-1 positions. e(n) is arbitrary; c c z contains the first row of the identity matrix. c c on output: c c d contains the eigenvalues in ascending order. if an c error exit is made, the eigenvalues are correct but c unordered for indices 1, 2, ..., ierr-1; c c e has been destroyed; c c z contains the first components of the orthonormal eigenvectors c of the symmetric tridiagonal matrix. if an error exit is c made, z contains the eigenvectors associated with the stored c eigenvalues; c c ierr is set to c zero for normal return, c j if the j-th eigenvalue has not been c determined after 30 iterations. c c ------------------------------------------------------------------ c integer i, j, k, l, m, n, ii, mml, ierr real*8 d(n), e(n), z(n), b, c, f, g, p, r, s, ep real*8 dsqrt, dabs, dsign, d1mach c ierr = 0 if (n .eq. 1) go to 1001 c e(n) = 0.0d0 do 240 l = 1, n j = 0 c :::::::::: look for small sub-diagonal element :::::::::: 105 do 110 m = l, n if (m .eq. n) go to 120 if (dabs(e(m)) .le. ep * (dabs(d(m)) + dabs(d(m+1)))) x go to 120 110 continue c 120 p = d(l) if (m .eq. l) go to 240 if (j .eq. 30) go to 1000 j = j + 1 c :::::::::: form shift :::::::::: g = (d(l+1) - p) / (2.0d0 * e(l)) r = dsqrt(g*g+1.0d0) g = d(m) - p + e(l) / (g + dsign(r, g)) s = 1.0d0 c = 1.0d0 p = 0.0d0 mml = m - l c c :::::::::: for i=m-1 step -1 until l do -- :::::::::: do 200 ii = 1, mml i = m - ii f = s * e(i) b = c * e(i) if (dabs(f) .lt. dabs(g)) go to 150 c = g / f r = dsqrt(c*c+1.0d0) e(i+1) = f * r s = 1.0d0 / r c = c * s go to 160 150 s = f / g r = dsqrt(s*s+1.0d0) e(i+1) = g * r c = 1.0d0 / r s = s * c 160 g = d(i+1) - p r = (d(i) - g) * s + 2.0d0 * c * b p = s * r d(i+1) = g + p g = c * r - b c :::::::::: form first component of vector :::::::::: f = z(i+1) z(i+1) = s * z(i) + c * f 200 z(i) = c * z(i) - s * f c d(l) = d(l) - p e(l) = g e(m) = 0.0d0 go to 105 240 continue c c :::::::::: order eigenvalues and eigenvectors :::::::::: do 300 ii = 2, n i = ii - 1 k = i p = d(i) c do 260 j = ii, n if (d(j) .ge. p) go to 260 k = j p = d(j) 260 continue c if (k .eq. i) go to 300 d(k) = d(i) d(i) = p p = z(i) z(i) = z(k) z(k) = p 300 continue c go to 1001 c :::::::::: set error -- no convergence to an c eigenvalue after 30 iterations :::::::::: 1000 ierr = l 1001 return c :::::::::: last card of gausq2 :::::::::: end 