libks.f90

I[1].M.Smith所著的《有限元方法编程》第三版Fortran程序

module libks 
contains
  subroutine lancz1(n,el,er,acc,leig,lx,lalfa,lp,itape,iflag,u,v, &
                     eig,jeig,neig,x,del,nu,alfa,beta)            
 implicit none
! use the lanczos algorithm to find the spectrum of a symmetric matrix
!     or that part of the spectrum that lies in a given interval. this
!     version returns information for later calculation of eigenvectors
!     by lancz2.
!  standard fortran 66 (a verified pfort subroutine)
!  lousy portabilty(ims) ;; converted to f90 (kind=2)
      integer ::   n,leig,lx,lalfa,iflag,lp,itape,neig
      integer ::   nu(lx), jeig(2,leig)
      real::   el,er,acc,u(0:n),v(0:n), eig(leig),x(lx),del(lx), &
                        alfa(lalfa),beta(lalfa)
!
! n must be set to the matrix order. it is not altered.
! el,er must be set to indicate the range within which the spectrum
!     is wanted. if el.ge.er then it is assumed that the whole spectrum
!     is wanted.they are not altered.
! acc must be set to the required precision, relative to the largest
!     eigenvalue of a. if it is very small or negative then as much
!     accuracy as the precision allows will be found. it is not altered.
! leig must be set to length of array eig. it must be as large as the
!     number of distinct eigenvalues in the range (el,er). it is not
!     altered.
! lx must be set to the length of arrays x,del,nu. a value three times
!     number of distinct eigenvalues in the range (el,er) usually
!     suffices. it is not altered.
! lalfa must be set to the length of arrays alfa and beta. it limits
!      the number of lanczos steps possible. it is not altered.
! lp must be set to the unit number for diagnostic messages. if lp.le.0
!     the messages are suppressed. it is not altered.
! itape must be set to the unit number of a sequential file if
!     eigenvectors are wanted from lancz2 or non-positive if they are
!     not wanted. it is not altered.
! iflag must be set prior to the first call for a particular matrix to
!       -1 if the user does not want to specify a start vector, or
!       -2 if the user has placed a required start vector in v.
!     it should not otherwise be changed.
!     on return it has the value
!       0 on successful completion
!       1 on a intermediate return
!       2,3,...,7 on error conditions
! u,v hold the working vectors of the lanczos recurrence. on a return
!     with iflag=1 the user must add to u the vector a*v without
!     altering v.
! eig need not be set by the user. on any return it contains neig
!     computed eigenvalues, stored in increasing order.
! jeig need not be set by the user. on return jeig(1,i) contains the
!     lanczos step at which eig(i) was accepted and jeig(2,i) contains
!     the matching point for the recurrences to get the eigenvector.
! neig need not be set by the user. on any return it contains the number
!     of eigenvalues in eig.
! x is used for a work array.  x(1).lt.x(2).lt.x(3).lt. ... are points
!     approximating eigenvalues of a.
! del is used as a workarray.  del(i),i=1,2,3,... contain the last pivot
!     in the lu factorization of t-lamda*i, where t is the lanczos
!     tridiagonal matrix, i is the identity matrix and lamda=x(i).
! nu is used as a workarray.  iabs(nu(i))-1,i=1,2,3,... contain the
!     number of negative pivots in the lu factorization of t-lamda*i,
!     lamda=x(i). fixed intervals x(i),x(i+1)  about converged
!     eigenvalues are flagged by negative values of nu(i).
! alfa and beta are used as workarrays holding the lanczos tridiagonal
!     matrix t. beta(1)=0.  and a typical row is
!         beta(i)     alfa(i)     beta(i+1)
      integer :: nfix = 0,nxtbnd,i,nlan,jlan1,jlan,mlan,l,kx,maxrz,         &
                 nuk,k=0,idummy,lfp,m=0,m1,j,ne=0,nxtr=0,match,nur,ncand,   &
                 nritz,maxrzo
      real    ::       drelpr,vv,uu,alf,bet,vi,tol,be2,anorm,oldel, &
                       older,dr,xmid,xa,xc,fa,fc,tem,den,pi,sig,            &
                       tau,disc,xav,pole,ritz,xl,xr,diff,                   &
                       err,bnd,t,tolc,w,xli,xri,en,g,fa01as
      common/ea15bd/anorm,oldel,older,tolc,kx,maxrz,maxrzo,nlan,mlan,nxtbnd
! common block ea15bd contains the following quantities,
!     which are required to be preserved between entries:
!   anorm is an estimate of the norm of the matrix a, based
!      on a gregorshin bound applied to t.
!   oldel,older are the values of el,er on the previous entry.
!   tolc is the tolerance for agreement between sucessive ritz
!      values that decides whether to call the error bounding routine
!      ea15cd (but see also nxtbnd)
!   kx is the number of points currently stored in x, with associated
!      information in del and nu.
!   maxrz is the maximal number of ritz values in a gap between
!      fixed intervals.
!   maxrzo is the previous value of maxrz
!   nlan is the number of lanczos steps performed. alfa(i)<
!      beta(i),i=1,2,...,nlan have been set.
!   mlan is the order of tridiagonal matrix so far used for this
!      interval of the spectrum.
!   nxtbnd is used to delay error bounding when tolc is at roundoff
!      error level. ea15cd is not called until we are looking
!      at the tridiagonal matrix of order nxtbnd.
      real  ::     zero,one,two,half,three
! drelpr is the relative precision
      data drelpr/2.2e-16/, zero/0.0e0/, one/1.0e0/, two/2.0e0/
!cray data drelpr/1.4e-14/, zero/0.0e0/, one/1.0e0/, two/2.0e0/
!ibm  data drelpr/2.2e-16/, zero/0.0e0/, one/1.0e0/, two/2.0e0/
      data half/0.5e0/, three/3.0e0/
      en=float(n)
      if(n.le.0 .or. lx.le.5 .or. lalfa.le.0)go to 950
      if(iflag.eq.1)go to 35
      vv=zero
      if(iflag.eq.-2)go to 5
      if(iflag.eq.-1)go to 15
      if(iflag.eq.0)go to 70
      go to 970
!
! find the norm of the user's start vector.
    5 do 10 i=1,n
      vv=vv+v(i)**2
   10 continue
      if(vv.gt.zero)go to 25
!
! generate pseudo-random start vector
   15 g=1431655765.e0
      do 20 i=1,n
      g=mod(g*9228907.e0,4294967296.e0)
      if(i.ge.0)fa01as=g/4294967296.e0
      if(i.lt.0)fa01as=2.e0*g/4294967296.e0-1.e0
      v(i)=fa01as
      vv=vv+v(i)**2
   20 continue  
!
! normalize start vector and perform initializations
   25 vv=one/ sqrt(vv)
      do 30 i=1,n
      v(i)=v(i)*vv
      u(i)=zero
   30 continue
      nlan=0
      anorm=zero
      beta(1)=zero
      neig=0
      go to 912
!
! perform a lanczos step
   35 nlan=nlan+1
      call ea15fd(n,itape,nlan,1,v)
      if(nlan.ge.lalfa)go to 940
      alf=zero
      do 40 i=1,n
      alf=alf+u(i)*v(i)
   40 continue
      alfa(nlan)=alf
      uu=zero
      do 50 i=1,n
      u(i)=u(i)-alf*v(i)
      uu=uu+u(i)**2
   50 continue
      uu= sqrt(uu)
      bet=beta(nlan)
      anorm=max(anorm,bet+ abs(alf)+uu)
      bet=max(uu,anorm*drelpr)
      uu=one/bet
      beta(nlan+1)=bet
      do 60 i=1,n
      vi=u(i)*uu
      u(i)=-bet*v(i)
      v(i)=vi
   60 continue
!
! this is the beginning of the loop that analyses t
! normally we advance the analysis of t by one row, but the
!     user may have requested a restart
   70 jlan1=nlan
      if(oldel.eq.el .and. older.eq.er .and. iflag.eq.1)go to 71
      jlan1=1
      neig=0
   71 do 910 jlan=jlan1,nlan
      mlan=jlan
      if(jlan-2<0)go to 73 ; if(jlan-2==0)go to 76; if(jlan-2>0)go to 85
!
! jlan=1. check for the trivial case.
   73 eig(1)=alfa(1)
      jeig(1,1)=jlan
      jeig(2,1)=jlan
      if(beta(2).gt.drelpr*anorm)go to 910
      if(el.ge.er .or. (el.le.eig(1) .and. eig(1).le.er))neig=1
      go to 960
!
! jlan=2. set up four points that span and separate  the
!     three ritz values at levels 1 and 2
   76 w=half* abs(alfa(1)-alfa(2))
      t= sqrt(w*w+beta(2)**2)
      x(2)=half*(alfa(1)+alfa(2)-t-w)
      x(3)=x(2)+t+w
      x(1)=x(2)+w-1.1*t
      x(4)=x(3)-w+1.1*t
      do 80 l=1,4
      del(l)=alfa(2)-x(l)-beta(2)**2/(alfa(1)-x(l))
      nu(l)=l-(l-1)/2
   80 continue
      kx=4
      maxrz=2
      maxrzo=1
      if(x(1).lt.x(2) .and. x(3).lt.x(4))go to 910
      eig(1)=(alfa(1)+alfa(2))*half
      jeig(1,1)=jlan
      jeig(2,1)=jlan
      neig=1
      go to 960
!
! jlan.gt.2
!
! add or remove points to the right, if appropriate
   85 nuk=jlan-1
      if(el.ge.er)go to 100
! find interval (x(k-1),x(k)) containing er
      do 92 i=2,kx
      k=2+kx-i
      if(x(k-1).le.er)go to 95
   92 continue
      k=1
   95 nuk=min0(iabs(nu(k)),jlan-1)
! look for point with nu value at least 2 greater than any in interval
!     containing er
      do 97 i=k,kx
      if(iabs(nu(i)).gt.nuk+1)go to 98
   97 continue
      go to 100
!  reduce kx to get rid of unnecessary points
   98 kx=i
      go to 110
! if necessary add extra points to right
  100 do 105 idummy=1,lx
      if(iabs(nu(kx-1)).gt.nuk)go to 110
      if(kx.ge.lx)go to 930
      kx=kx+1
      x(kx)=x(kx-1)*three-x(kx-2)*two
      call ea15dd(jlan-1,alfa,beta,x(kx),del(kx),nu(kx),dr,nur)
  105  continue
!
! copy information to ends of arrays x,del,nu. if any
!     gap between fixed intervals contains no ritz values,
!     remove all points in that gap
  110 lfp=lx
! lfp is the last left-end of a fixed interval
      m=lx
      do 140 idummy=1,lx
      if(nu(kx-1).ge.0)go to 130
      if(nu(kx).eq.iabs(nu(lfp)))m=lfp-1
      lfp=m-1
  130 x(m)=x(kx)
      del(m)=del(kx)
      nu(m)=nu(kx)
      kx=kx-1
      m=m-1
      if(kx.eq.1)go to 150
  140 continue
  150 x(m)=x(1)
      del(m)=del(1)
      nu(m)=nu(1)
!
! add or remove points to the left, if appropriate
      nuk=2
      if(el.ge.er)go to 175
! find interval (x(k-1),x(k)) containing el
      m1=m+1
      do 155 k=m1,lx
      if(x(k).ge.el)go to 160
  155 continue
      k=lx+1
  160 nuk=max0(iabs(nu(k-1)),2)
! look for point with nu value at least 2 less than any in interval
!     containing el
      do 165 j=m1,k
      i=m+k-j
      if(iabs(nu(i)).lt.nuk-1)go to 170
  165 continue
      go to 175
! increase m to get rid of unnecessary points
  170 m=i
      go to 190
! if necessary, add extra points to left
  175 do 180 idummy=1,lx
      if(iabs(nu(m+1)).lt.nuk)go to 190
      if(m.eq.3)go to 930
      m=m-1
      x(m)=three*x(m+1)-two*x(m+2)
      call ea15dd(jlan-1,alfa,beta,x(m),del(m),nu(m),dr,nur)
  180 continue
  190 k=3
      if(m.le.5)go to 930
      xri=x(m)
      alf=alfa(jlan)
      bet=beta(jlan)
      be2=bet**2
! tol holds the radius of fixed intervals set up around converged
!     eigenvalues.
      tol=max(acc,en*two*drelpr)*anorm
! tolc is the tolerance used for accepting eigenvalues.
      if(jlan.lt.10)tolc=max(tol**2*en/anorm,anorm*drelpr*5.)
! move leading three points to beginning of store and update
!     associated data.
      m=m-1
      do  210 i=1,3
        m=m+1
        x(i)=x(m)
        nu(i)=nu(m)