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