libks.f90

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

!     lanczos tridiagonal matrix t. they are not altered.  note that a
!     typical row of t is
!          beta(i)  alfa(i)  beta(i+1)
! e must be set to the approximate eigenvalue of t. it is not altered.
! enorm must be set to the relative precision times the norm of a times
!     the order of a. it is not altered.
! err is set to a bound for the error of e.
! bnd is set to an estimate of the amount that err might be reduced by
!     using a more accurate e and/or a more accurate eigenvector
!     calculation.
! match is set to the matching point between use of forward and backward
!     recurrences.
   implicit none
      integer ::       lalfa,j,match
      real    ::       alfa(lalfa),beta(lalfa),e,s,psi,bnd,err,enorm
      real    ::       w,w0,w1,w2,sw,z1,z2,bndk,errk,wmax10  ,zero,one
      integer ::       j1,k,ktry,k1,kk,i
      data zero/0.0e0/, one/1.0e0/
      j1=j-1
!
! start the forward recurrence for solving (t-e)*z=psi*ek
      w1=zero
      w2=one
      sw=zero
      k=0
      w=alfa(1)-e
      wmax10=zero
!
      do 100 ktry=1,10
      if(k-j1<0)goto 5; if(k-j1==0)goto 15; if(k-j1>0)goto 110
    5 k1=k+1
!     continue the forward recurrence until the first maximum
!     or a maximum at least ten times bigger than the previous
!     maximum is found.
      do 10 k=k1,j1
      w0=w1
      w1=w2
      w2=-w/beta(k+1)
      sw=sw+w1**2
      w=(alfa(k+1)-e)*w2+beta(k+1)*w1
      if( abs(w1).lt.wmax10)go to 10
      if( abs(w1).ge.max( abs(w0), abs(w2)))go to 20
   10 continue
   15 k=j
      w0=w1
      w1=w2
      sw=sw+w1**2
      if( abs(w1).lt.wmax10)go to 100
      if( abs(w1).lt. abs(w0))go to 100
!
!     forward recurrence shows a maximum at component k. continue
!     solving (t-e)*z=psi*ek by backward recurrence
   20 s=zero
      wmax10= abs(w1)*10.
      z1=one
      psi=alfa(j)-e
      if(k.gt.j1)go to 40
      do 30 kk=k,j1
      i=k+j1-kk
      z2=z1
      z1=-psi/beta(i+1)
      s=s+z2**2
      psi=(alfa(i)-e)*z1+beta(i+1)*z2
   30 continue
   40 if(w1.eq.zero)go to 100
!
! form present error bound and store the best so far found
      z1=z1/w1
      s=s+sw*z1**2
      s= sqrt(s)
      if(k.gt.1)psi=psi+beta(k)*z1*w0
      bndk=1.25* abs(psi)/s
      errk=1.25*(enorm+beta(j+1)/s)+bndk
      if(ktry.eq.1)go to 50
      if(err.lt.errk)go to 90
   50 err=errk
      bnd=bndk
      match=k
!
! test result for acceptability
   90 if(bndk.lt.errk/5.)go to 110
  100 continue
  110 return
      end subroutine ea15cd
      subroutine ea15dd(n,alfa,beta,x,del,nu,dr,nur)
! computes the triangular factorization of t-x*i, where t is a
!     tridiagonal matrix, x is a scalar and i is the identity matrix. it
!     records the number of negative pivots and the last two pivots.
  implicit none
      integer :: n,nu,nur
      real    :: alfa(n),beta(n),x,del,dr
! n must be set to the order of the matrix. it is not altered.
! alfa must be set ot contain the diagonal elements of the tridiagonal
!     matrix. it is not altered.
! beta(1) must have the value zero and beta(2),...,beta(nlan) must be
!     set to contain the off-diagonal elements. beta is not altered.
! x   must be set to the scalar. it is not altered.
! del is set to the last pivot.
! nu is set to 1+(the number of negative pivots).
! dr is set to the last-but-one pivot.
! nur is set to 1+(the number of negative pivots when the last is
!     excluded).
      real ::      zero,one        ; integer :: k
! drelpr is the relative precision.
      real ::      drelpr
!ibm  data drelpr/2.2e-16/
!ray  data drelpr/1.4e-14/
      data drelpr/2.2e-16/
      data zero/0.0e0/, one/1.0e0/
      nu=1
      del=one
      do 10 k=1,n
      dr=del
      del=alfa(k)-x-beta(k)**2/del
      if(del)6,3,10
    3 del=drelpr*( abs(alfa(k))+ abs(x))
      go to 10
    6 nu=nu+1
   10 continue
   20 nur=nu
      if(del.lt.zero)nur=nu-1
      return
      end subroutine ea15dd
      subroutine lancz2(n,lalfa,lp,itape,eig,jeig,neig,  &
                       alfa,beta,lz,jflag,y,w,z)
! find approximate eigenvectors corresponding to eigenvalues found
!     by lancz1.
! n must be set to the matrix order. it is not altered.
! lalfa must be set to the length of arrays alfa and beta. 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 the sequential file used
!     by lancz1. it is not altered.
! eig must contain eigenvalues as found by lancz1. it is not altered.
! jeig must hold the information passed down by lancz1, namely jeig(1,i)
!     must contain  the lanczos step at which eig(i) was accepted and
!     jeig(2,i) must contain the matching point for the eigenvector
!     corresponding to eig(i). jeig is not altered.
! neig must hold the number of eigenvectors wanted. it is not altered.
! alfa,beta must be as left by lancz1. they are not altered.
! lz must hold the first dimension of z. it must be at least
!     max(jeig(1,i),i=1,neig). it is not altered.
! iflag need not be set. on return it has one of the values
!      0    successful completion
!      1    n.le.0 .or. lalfa.le.0 .or. itape.le.0 .or. neig.le.0
!      2    ly too small
!      3    lz too small
! y is set to the required eigenvectors.
! w is used as a workvector.
! z is used as a workvector.
! lz must be set to the first dimension of z. it must be at least
!     max(jeig(1,i),i=1,2,...,neig).
   implicit none 
      integer ::       neig,lalfa,n,lz,jflag,itape,lp  
      real    ::       eig(neig),alfa(lalfa),beta(lalfa),   &
                       y(0:n,neig),w(0:n),z(lz,neig)
      integer ::       jeig(2,neig)
      real    ::       s,zero,one  ; integer :: i,j,l,m
      data zero/0.0e0/, one/1.0e0/
      jflag=0
      if(n.le.0 .or. lalfa.le.0 .or. itape.le.0 .or.neig.le.0)go to 100
      if(n.lt.n)go to 120
      m=0
      do 20 l=1,neig
      m=max0(m,jeig(1,l))
      if(lz.lt.m)go to 140
      do 10 i=1,n
      y(i,l)=zero
   10 continue
      call ea15gd(jeig(1,l),lalfa,alfa,beta,eig(l),jeig(2,l),z(1,l))
   20 continue
      do 50 j=1,m
      call ea15fd(n,itape,j,2,w)
      do 40 l=1,neig
      if(j.gt.jeig(1,l))go to 40
      do 30 i=1,n
      y(i,l)=y(i,l)+z(j,l)*w(i)
   30 continue
   40 continue
   50 continue
!
! normalize the vectors
      do 90 l=1,neig
      s=zero
      do 70 i=1,n
      s=s+y(i,l)**2
   70 continue
      s=one/ sqrt(s)
      do 80 i=1,n
      y(i,l)=y(i,l)*s
   80 continue
   90 continue
      go to 150
  100 jflag=1
      if(lp.gt.0)write(lp,110)'error return from lancz2. &
          & n=',n,'lalfa=',lalfa,'itape=',itape,'neig=',neig
  110 format( a,i6,a,i6,a,i6,a,i6 )
      go to 150
  120 jflag=2
      if(lp.gt.0)write(lp,130)'error return 2 from lancz2. ly =',n,&
         & 'and should be at least',n
  130 format( a , i6 , a , i6)
      go to 150
  140 jflag=3
      if(lp.gt.0)write(lp,145)'error return 3 from lancz2. lz =',lz,&
         & 'and should be at least',m
  145 format( a , i6, a , i6)
  150 return
      end  subroutine lancz2
      subroutine ea15fd(n,itape,nlan,io,v)
! store or recover a lanczos vector.
   implicit none
      integer :: n,itape,nlan,io
      real    :: v(0:n)
      integer :: last,ngap,i
! n is the order of the matrix (also the order of the vector being
!     stored or recovered). it must not be altered.
! itape is the unit number of the tape unit. it is not altered.
! nlan is the lanczos iteration number. the vectors are always sent for
!     storage in order (nlan=1,2,...) and are recovered in order
!     (nlan=1,2,...). if the first vector (nlan=1) is sent for storage
!     then it can be assumed that a new problem is in hand and the old
!     vectors may be overwritten. nlan must not be altered.
! io is 1 if a vector is to be stored and to 2 if one is to be
!     recovered. it must not be altered.
! v holds the vector to be stored or is set to the vector recovered.
      common/ea15hd/last
! last is the value of nlan on the last call.
      if(itape.le.0)go to 50
      if(nlan.ne.1)go to 10
      rewind itape
      last=0
! if necessary, skip over the vectors written earlier
   10 ngap=nlan-last-1
      if(ngap.le.0)go to 30
      do 20 i=1,ngap
      read(itape)
   20 continue
! perform actual read or write
   30 if(io.eq.1)write(itape)v
      if(io.eq.2)read(itape)v
   40 last=nlan
   50 return
      end  subroutine ea15fd
      subroutine ea15gd(j,lalfa,alfa,beta,e,match,z)
! find an approximate eigenvector of a lanczos tridiagonal matrix by
!     forward and backward recurrence.
! j is the order of the matrix and is not altered.
! alfa,beta must be exactly as left by lancz1. they are not altered.
! e is the eigenvalue corresponding to which an eigenvector is wanted.
!     it is not altered.
! match is the matching point between the recurrences. it is not
!     altered.
! z is set to the required eigenvector.
  implicit none 
      integer ::   j,lalfa,match
      real    ::   alfa(lalfa),beta(lalfa),e,z(j)
      real    ::   psi,w,w0,w1,w2,z1,z2  ,zero,one
      integer ::   j1,k,kk,i 
      data zero/0.0e0/, one/1.0e0/
      j1=j-1
!
! forward recurrence.
      w1=zero
      w2=one
      w=alfa(1)-e
      if(j1.le.0)go to 15
      do 10 k=1,j1
      w0=w1
      w1=w2
      w2=-w/beta(k+1)
      w=(alfa(k+1)-e)*w2+beta(k+1)*w1
      z(k)=w1
      if(k.eq.match)go to 20
   10 continue
   15 k=j
      w0=w1
      w1=w2
      z(k)=w1
! backward recurrence
   20 z1=one
      psi=alfa(j)-e
      if(k.gt.j1)go to 40
      do 30 kk=k,j1
      i=k+j1-kk
      z2=z1
      z1=-psi/beta(i+1)
      z(i+1)=z2
      psi=(alfa(i)-e)*z1+beta(i+1)*z2
   30 continue
! rescale the last set of components
   40 w1=w1/z1
      if(k.gt.j1)go to 60
      do 50 i=k,j1
      z(i+1)=z(i+1)*w1
   50 continue
   60 return
      end subroutine ea15gd 
end module libks