c  m3xxx98.f    sept 1998
c
c           poles of Pade of function involving sqrt(P(z))
c
      PROGRAM M3xxx98
      parameter(ndim=200)
      INTEGER twom,m,i,nmax
      complex X(20),A(ndim),B2(ndim),b(ndim),root(ndim)
      double precision hr(ndim,ndim),hi(ndim,ndim),wr(ndim),wi(ndim)
      READ *,m,nmax
      twom=2*m
c  2m zeros of P:
      do i=1,twom
       read *,xr,yr
       x(i)=cmplx(xr,yr)
      enddo
       CALL recsqr(m,X,A,B2,NMAX)
       amaxb=0
       do i=1,nmax
        b(i)=sqrt(b2(i))
        amaxb=max(amaxb, abs(b(i)) )
       enddo
       print '(a1,''psset{unit=1cm}'')',char(92)
       print '(a1,''begin{pspicture}(17,2.5)'')',char(92)
       do i=1,nmax
       print '(a1,''psdots[dotsize=0.15]('',f7.3,'','',f7.3,'')'')',
     &  char(92),i*17.0/nmax,0.2+2*real(b(i))/amaxb
       enddo
       print '(a1,''end{pspicture}'')',char(92)
       PRINT '(10F8.4)',(A(I),I=1,5)
       PRINT '(10F8.4)',(B(I),I=1,5)
c     call trilr(nmax,ntrue,b,a,root)
c     print *,ntrue,' found'
c      PRINT '(10F8.4)',(root(I),I=1,ntrue)
      do i=1,nmax
      do j=1,nmax
       hr(i,j)=0
       hi(i,j)=0
      enddo
      hr(i,i)=a(i)
      hi(i,i)=aimag( a(i) )
      enddo
      hr(1,2)=b(1)
      hi(1,2)=aimag( b(1) )
      do i=2,nmax-1
       hr(i,i-1)=b(i-1)
       hi(i,i-1)=aimag( b(i-1) )
       hr(i,i+1)=b(i)
       hi(i,i+1)=aimag( b(i) )
      enddo
      hr(nmax,nmax-1)=b(nmax-1)
      hi(nmax,nmax-1)=aimag( b(nmax-1) )
      low=1
      igh=nmax
      call comlr(ndim,nmax,low,igh,hr,hi,wr,wi,ierr)
      if(ierr.ne.0)print *,' err: ',ierr
      print '(10f8.4)',(wr(i),i=1,nmax)
      print '(10f8.4)',(wi(i),i=1,nmax)
      amaxx=wr(1)
      aminx=wr(1)
      amaxy=wi(1)
      aminy=wi(1)
      do i=2,nmax
        amaxx=max(amaxx, wr(i) )
        aminx=min(aminx, wr(i) )
        amaxy=max(amaxy, wi(i) )
        aminy=min(aminy, wi(i) )
      enddo
       print '(a1,''psset{unit=1cm}'')',char(92)
       print '(a1,''begin{pspicture}(17,5)'')',char(92)
       do i=1,twom
       print '(a1,''psdots[dotstyle=x,dotsize=0.25]
     & ('',f7.3,'','',f7.3,'')'')',
     &  char(92),16.0*(real(x(i))-aminx)/(amaxx-aminx),
     &           5.0* (aimag(x(i))-aminy)/(amaxy-aminy)
       enddo
       do i=1,nmax
       print '(a1,''psdots[dotsize=0.15]('',f7.3,'','',f7.3,'')'')',
     &  char(92),16.0*(wr(i)-aminx)/(amaxx-aminx),
     &           5.0* (wi(i)-aminy)/(amaxy-aminy)
       enddo
       print '(a1,''end{pspicture}'')',char(92)
      STOP
      END
      SUBROUTINE recsqr(P,X,A,B,NMAX)
C  recurrence coeff  for
C                    [x ,x ] [x ,x ] ... [x    ,x  ]
C                      1  2    3  4        2p-1  2p
C
C  A partir de p, x (vecteur des 2p xi) et nmax, la routine
C  donne A et B, vecteurs des coefficients de recurrence
      INTEGER P,N,NMAX
      complex X(1),A(1),B(1),
     &       R(300),S(300),PTS(300),R0(300),S0(300),
     &       LP1(300),LP2(300),XPTS(300),LPP(20),LPP1,
     &       SOMX,X0,EPS,EPS1,PB,R10,B0,A0,RL10,PI
      SOMX=0
      DO 11 I=1,P
   11 PTS(I)=(X(2*I-1)+X(2*I))/2
      A0=0
      SOMX=0
      DO 115 I=1,P
       SOMX=SOMX+X(2*I-1)+X(2*I)
       IF(I.LT.P)A0=A0-X(2*I)
       R(I)=1
       XPTS(I)=1
       LP1(I)=1
       LP2(I)=0
       DO 1151 J=1,P
        IF(J.LT.P)R(I)=R(I)*( PTS(I)-X(2*J) )
        XPTS(I)=XPTS(I)*(PTS(I)-X(2*J-1))*(PTS(I)-X(2*J))
        IF(J.EQ.I)GOTO 1151
        LP1(I)=LP1(I)/(PTS(I)-PTS(J))
        LP2(I)=LP2(I)-PTS(J)
 1151  CONTINUE
 115   LP2(I)=LP2(I)*LP1(I)
      A0=A0+SOMX/2
      DO 14 I=1,P
      S(I)=R(I)*(PTS(I)-A0)
      R0(I)=R(I)
  14  S0(I)=S(I)
      A(1)=A0
C
C  Iterations
C
      DO 2 IT=1,NMAX
       ID=4
       IF(IT.EQ.1)ID=2
       DO 21 I=1,P
   21  R(I)=(S(I)**2-XPTS(I))/(ID*R(I))
       B0=0
       A0=0
       DO 22 I=1,P
       A0=A0+R(I)*LP2(I)
   22  B0=B0+R(I)*LP1(I)
       A0=A0/B0+SOMX/2
       DO 23 I=1,P
       R(I)=R(I)/B0
   23  S(I)=2*(PTS(I)-A0)*R(I)-S(I)
       A(IT+1)=A0
       B(IT)=B0
 2    CONTINUE
      RETURN
      END
      subroutine comlr(nm,n,low,igh,hr,hi,wr,wi,ierr)
c
cReceived: from CUNYVM.BITNET by Pythag.ucl.ac.be (Mailer R2.04) with BSMTP id
c 2883; Mon, 13 Nov 89 09:32:44 +0100
cReceived: from CUNYVM by CUNYVM.BITNET (Mailer R2.03B) with BSMTP id 2335; Mon,
c 13 Nov 89 02:56:23 EST
cReceived: from research.att.com by CUNYVM.CUNY.EDU (IBM VM SMTP R1.2.1MX) with
c TCP; Mon, 13 Nov 89 02:56:16 EST
cFrom: netlibd@research.att.com
cDate: Mon, 13 Nov 89 02:57 EST
cTo: annu76%buclln11.bitnet@cunyvm.cuny.edu
cSubject: send comlr from eispack
c
cCareful! Anything free comes with no guarantee.
c*** from netlib, Mon Nov 13 02:57:26 EST 1989 ***
      integer i,j,l,m,n,en,ll,mm,nm,igh,im1,itn,its,low,mp1,enm1,ierr
      double precision hr(nm,n),hi(nm,n),wr(n),wi(n)
      double precision si,sr,ti,tr,xi,xr,yi,yr,zzi,zzr,tst1,tst2
c
c     this subroutine is a translation of the algol procedure comlr,
c     num. math. 12, 369-376(1968) by martin and wilkinson.
c     handbook for auto. comp., vol.ii-linear algebra, 396-403(1971).
c
c     this subroutine finds the eigenvalues of a complex
c     upper hessenberg matrix by the modified lr method.
c
c     on input
c
c        nm must be set to the row dimension of two-dimensional
c          array parameters as declared in the calling program
c          dimension statement.
c
c        n is the order of the matrix.
c
c        low and igh are integers determined by the balancing
c          subroutine  cbal.  if  cbal  has not been used,
c          set low=1, igh=n.
c
c        hr and hi contain the real and imaginary parts,
c          respectively, of the complex upper hessenberg matrix.
c          their lower triangles below the subdiagonal contain the
c          multipliers which were used in the reduction by  comhes,
c          if performed.
c
c     on output
c
c        the upper hessenberg portions of hr and hi have been
c          destroyed.  therefore, they must be saved before
c          calling  comlr  if subsequent calculation of
c          eigenvectors is to be performed.
c
c        wr and wi contain the real and imaginary parts,
c          respectively, of the eigenvalues.  if an error
c          exit is made, the eigenvalues should be correct
c          for indices ierr+1,...,n.
c
c        ierr is set to
c          zero       for normal return,
c          j          if the limit of 30*n iterations is exhausted
c                     while the j-th eigenvalue is being sought.
c
c     calls cdiv for complex division.
c     calls csroot for complex square root.
c
c     questions and comments should be directed to burton s. garbow,
c     mathematics and computer science div, argonne national laboratory
c
c     this version dated august 1983.
c
c     ------------------------------------------------------------------
c
      ierr = 0
c     .......... store roots isolated by cbal ..........
      do 200 i = 1, n
         if (i .ge. low .and. i .le. igh) go to 200
         wr(i) = hr(i,i)
         wi(i) = hi(i,i)
  200 continue
c
      en = igh
      tr = 0.0d0
      ti = 0.0d0
      itn = 30*n
c     .......... search for next eigenvalue ..........
  220 if (en .lt. low) go to 1001
      its = 0
      enm1 = en - 1
c     .......... look for single small sub-diagonal element
c                for l=en step -1 until low d0 -- ..........
  240 do 260 ll = low, en
         l = en + low - ll
         if (l .eq. low) go to 300
         tst1 = dabs(hr(l-1,l-1)) + dabs(hi(l-1,l-1))
     x            + dabs(hr(l,l)) + dabs(hi(l,l))
         tst2 = tst1 + dabs(hr(l,l-1)) + dabs(hi(l,l-1))
         if (tst2 .eq. tst1) go to 300
  260 continue
c     .......... form shift ..........
  300 if (l .eq. en) go to 660
      if (itn .eq. 0) go to 1000
      if (its .eq. 10 .or. its .eq. 20) go to 320
      sr = hr(en,en)
      si = hi(en,en)
      xr = hr(enm1,en) * hr(en,enm1) - hi(enm1,en) * hi(en,enm1)
      xi = hr(enm1,en) * hi(en,enm1) + hi(enm1,en) * hr(en,enm1)
      if (xr .eq. 0.0d0 .and. xi .eq. 0.0d0) go to 340
      yr = (hr(enm1,enm1) - sr) / 2.0d0
      yi = (hi(enm1,enm1) - si) / 2.0d0
      call csroot(yr**2-yi**2+xr,2.0d0*yr*yi+xi,zzr,zzi)
      if (yr * zzr + yi * zzi .ge. 0.0d0) go to 310
      zzr = -zzr
      zzi = -zzi
  310 call cdiv(xr,xi,yr+zzr,yi+zzi,xr,xi)
      sr = sr - xr
      si = si - xi
      go to 340
c     .......... form exceptional shift ..........
  320 sr = dabs(hr(en,enm1)) + dabs(hr(enm1,en-2))
      si = dabs(hi(en,enm1)) + dabs(hi(enm1,en-2))
c
  340 do 360 i = low, en
         hr(i,i) = hr(i,i) - sr
         hi(i,i) = hi(i,i) - si
  360 continue
c
      tr = tr + sr
      ti = ti + si
      its = its + 1
      itn = itn - 1
c     .......... look for two consecutive small
c                sub-diagonal elements ..........
      xr = dabs(hr(enm1,enm1)) + dabs(hi(enm1,enm1))
      yr = dabs(hr(en,enm1)) + dabs(hi(en,enm1))
      zzr = dabs(hr(en,en)) + dabs(hi(en,en))
c     .......... for m=en-1 step -1 until l do -- ..........
      do 380 mm = l, enm1
         m = enm1 + l - mm
         if (m .eq. l) go to 420
         yi = yr
         yr = dabs(hr(m,m-1)) + dabs(hi(m,m-1))
         xi = zzr
         zzr = xr
         xr = dabs(hr(m-1,m-1)) + dabs(hi(m-1,m-1))
         tst1 = zzr / yi * (zzr + xr + xi)
         tst2 = tst1 + yr
         if (tst2 .eq. tst1) go to 420
  380 continue
c     .......... triangular decomposition h=l*r ..........
  420 mp1 = m + 1
c
      do 520 i = mp1, en
         im1 = i - 1
         xr = hr(im1,im1)
         xi = hi(im1,im1)
         yr = hr(i,im1)
         yi = hi(i,im1)
         if (dabs(xr) + dabs(xi) .ge. dabs(yr) + dabs(yi)) go to 460
c     .......... interchange rows of hr and hi ..........
         do 440 j = im1, en
            zzr = hr(im1,j)
            hr(im1,j) = hr(i,j)
            hr(i,j) = zzr
            zzi = hi(im1,j)
            hi(im1,j) = hi(i,j)
            hi(i,j) = zzi
  440    continue
c
         call cdiv(xr,xi,yr,yi,zzr,zzi)
         wr(i) = 1.0d0
         go to 480
  460    call cdiv(yr,yi,xr,xi,zzr,zzi)
         wr(i) = -1.0d0
  480    hr(i,im1) = zzr
         hi(i,im1) = zzi
c
         do 500 j = i, en
            hr(i,j) = hr(i,j) - zzr * hr(im1,j) + zzi * hi(im1,j)
            hi(i,j) = hi(i,j) - zzr * hi(im1,j) - zzi * hr(im1,j)
  500    continue
c
  520 continue
c     .......... composition r*l=h ..........
      do 640 j = mp1, en
         xr = hr(j,j-1)
         xi = hi(j,j-1)
         hr(j,j-1) = 0.0d0
         hi(j,j-1) = 0.0d0
c     .......... interchange columns of hr and hi,
c                if necessary ..........
         if (wr(j) .le. 0.0d0) go to 580
c
         do 540 i = l, j
            zzr = hr(i,j-1)
            hr(i,j-1) = hr(i,j)
            hr(i,j) = zzr
            zzi = hi(i,j-1)
            hi(i,j-1) = hi(i,j)
            hi(i,j) = zzi
  540    continue
c
  580    do 600 i = l, j
            hr(i,j-1) = hr(i,j-1) + xr * hr(i,j) - xi * hi(i,j)
            hi(i,j-1) = hi(i,j-1) + xr * hi(i,j) + xi * hr(i,j)
  600    continue
c
  640 continue
c
      go to 240
c     .......... a root found ..........
  660 wr(en) = hr(en,en) + tr
      wi(en) = hi(en,en) + ti
      en = enm1
      go to 220
c     .......... set error -- all eigenvalues have not
c                converged after 30*n iterations ..........
 1000 ierr = en
 1001 return
      end
      subroutine csroot(xr,xi,yr,yi)
      double precision xr,xi,yr,yi
c
c     (yr,yi) = complex dsqrt(xr,xi)
c     branch chosen so that yr .ge. 0.0 and sign(yi) .eq. sign(xi)
c
      double precision s,tr,ti,pythag
      tr = xr
      ti = xi
      s = dsqrt(0.5d0*(pythag(tr,ti) + dabs(tr)))
      if (tr .ge. 0.0d0) yr = s
      if (ti .lt. 0.0d0) s = -s
      if (tr .le. 0.0d0) yi = s
      if (tr .lt. 0.0d0) yr = 0.5d0*(ti/yi)
      if (tr .gt. 0.0d0) yi = 0.5d0*(ti/yr)
      return
      end
      subroutine cdiv(ar,ai,br,bi,cr,ci)
      double precision ar,ai,br,bi,cr,ci
c
c     complex division, (cr,ci) = (ar,ai)/(br,bi)
c
      double precision s,ars,ais,brs,bis
      s = dabs(br) + dabs(bi)
      ars = ar/s
      ais = ai/s
      brs = br/s
      bis = bi/s
      s = brs**2 + bis**2
      cr = (ars*brs + ais*bis)/s
      ci = (ais*brs - ars*bis)/s
      return
      end
      double precision function pythag(a,b)
      double precision a,b
c
c     finds dsqrt(a**2+b**2) without overflow or destructive underflow
c
      double precision p,r,s,t,u
      p = dmax1(dabs(a),dabs(b))
      if (p .eq. 0.0d0) go to 20
      r = (dmin1(dabs(a),dabs(b))/p)**2
   10 continue
         t = 4.0d0 + r
         if (t .eq. 4.0d0) go to 20
         s = r/t
         u = 1.0d0 + 2.0d0*s
         p = u*p
         r = (s/u)**2 * r
      go to 10
   20 pythag = p
      return
      end
