libks.f90

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

        del(i)=alf-x(m)-be2/del(m)
        if(del(i))195,200,210
   195  nu(i)=nu(i)+isign(1,nu(i))
        go to 210
   200  del(i)=drelpr*anorm
   210 continue
!
! process the points from left to right. the current interval
!     is (x(k-1),x(k))=(x(m-1),x(m)), k.lt.m-2.
! nu(i),del(i),i=1,2,...,k hold new values of nu,delta.
! nu(i),del(i),i=m-2,m-1,...,lx hold old values.
      ne=1
! ne normally holds the number of empty intervals adjacent on the
!     left of the current interval.
      lfp=1
! lfp points to the last fixed point encountered.
!      do 530 until m.eq.lx
  220 if(m.ge.lx)go to 600
      if(nu(k-1).ge.0)go to 240
! we have a fixed interval
      if(iabs(nu(k-1)).eq.iabs(nu(k)))go to 230
! accept fixed interval
      lfp=k
      go to 500
! fixed interval no longer contains ritz value. free it.
  230 nu(k-1)=iabs(nu(k-1))
      do 235 i=1,neig
      if(eig(i).le.x(k))go to 235
      eig(i-1)=eig(i)
      jeig(1,i-1)=jeig(1,i)
      jeig(2,i-1)=jeig(2,i)
  235 continue
      neig=neig-1
  240 if(nu(k-1).lt.iabs(nu(k)))go to 250
! current interval contains no ritz value.
      ne=ne+1
      if(ne.le.3)go to 502
! there are four adjacent empty intervals. remove middle point.
      x(k-2)=x(k-1)
      nu(k-2)=nu(k-1)
      del(k-2)=del(k-1)
      x(k-1)=x(k)
      nu(k-1)=nu(k)
      del(k-1)=del(k)
      ne=3
      k=k-1
      go to 502
  250 if(iabs(nu(m-1)).eq.iabs(nu(m)))go to 500
!
! interval is 'interesting'. it contains at least one new and
!     at least one old ritz value.
! jump if interval is a gap between fixed intervals,
!     and contains just one ritz value, unless this is the
!     first time that all such gaps have less than two
!     ritz values.
      if(maxrz.le.1 .and. maxrzo.gt.1)go to 252
      if(nu(k-1)+nu(k).eq.-1)go to 500
  252 xmid=(x(m-1)+x(m))*half
! bisect if interval contains more than one ritz value.
      ritz=xmid
      if(iabs(nu(k)).gt.nu(k-1)+1)go to 275
! check whether progress is so slow that bisection is needed
!     the criterion is that three steps have been taken without
!     reducing the interval length by at least the factor 0.6.
!     (xli,xri) is the original interval of interest or the
!     interval of interest when it was last reduced in length
!      by at least the factor 0.6
      if(x(k-1).ge.xri)go to 255
      if(x(k)-x(k-1).ge.0.6*(xri-xli))go to 257
  255 nxtr=0
      xri=x(k-1)
      xli=x(k)
! bisect if progress on this interval is slow
  257 if(nxtr.ge.2)go to 275
      nxtr=nxtr+1
! estimate position (pole) of old ritz value by 2-1 rational
!     interpolation at x(l-1), x(l), x(l+1)
! choose l so that extra point is near
      l=m
      if(xmid-x(m-2).lt.x(m+1)-xmid)l=m-1
      xa=x(l-1)-x(l)
      xc=x(l+1)-x(l)
      fa=xa*(xa+del(l-1))
      fc=xc*(xc+del(l+1))
      tem=xa/(xa-xc)
      den=del(l)-del(l-1)-tem*(del(l+1)-del(l-1))
      if(den.eq.zero)go to 275
      pi=(tem*(fa-fc)-fa)/den
      sig=half*((fc-fa)-pi*(del(l+1)-del(l-1)))/(xc-xa)
      tau=fa-two*xa*sig-del(l-1)*pi
      disc=sig**2+tau
      if(disc.lt.zero)go to 275
      xav=sig+ sign( sqrt(disc),sig)
      pole=-tau/xav
! jump if pole is in required interval
      if( abs(pole+x(l)-xmid).le.x(m)-xmid)go to 260
! try other root.
      pole=xav
      if( abs(pole+x(l)-xmid).gt.x(m)-xmid)go to 275
! estimate position (ritz) of ritz value by 2-1 rational
!     interpolation at x(k-1),x(k) using (x-pole) for
!     denominator.
  260 if(l.ne.m)go to 265
      tau=-pole*del(k)
      sig=xa+((xa-pole)*del(k-1)-tau)/xa
      go to 270
  265 tau=-pole*del(k-1)
      sig=xc+((xc-pole)*del(k)-tau)/xc
  270 sig=sig*half
      disc=sig**2+tau
      if(disc.lt.zero)go to 275
      xav=sig+ sign( sqrt(disc),sig)
      tau=-tau/xav
      if( abs(tau+x(l)-xmid).le.x(m)-xmid)go to 280
      tau=xav
      if( abs(tau+x(l)-xmid).le.x(m)-xmid)go to 280
! calculation has failed
! if third point is just outside current interval then take
!     point twice as far inside interval. otherwise bisect.
  275 ritz=xmid
      pole=xmid
      if(x(m)-x(m-1).le.tolc)go to 300
      go to 490
  280 ritz=x(l)+tau
      diff= abs(tau-pole)
      pole=x(l)+pole
      if( abs(tau).ge.50.*max(diff,tolc))go to 480
      if(diff.gt.tolc)go to 500
! we may have a converged eigenvalue
  300 if(jlan.lt.nxtbnd)go to 500
      call ea15cd(jlan,lalfa,alfa,beta,ritz,     &
                   anorm*drelpr*en,err,bnd,match)
      if(bnd.lt.err*0.1)tolc=diff*(tol/err)**2
      if(tolc.gt.anorm*drelpr*5.0)go to 305
      tolc=anorm*drelpr*5.0
! tolc has become too small to give good criterion for calling ea15cd.
!    do not call ea15cd again until 1% more steps performed.
      nxtbnd=jlan+jlan/100
  305 if(err.le.tol)go to 310
      if(bnd.gt.err*0.1 .and. x(m-1).lt.ritz .and. ritz.lt.x(m)) go to 480
      go to 500
!
! we have an accepted point
! test whether (ritz-err,ritz+err) overlaps a fixed interval
  310 if(ritz-err.lt.x(lfp) .and. lfp.gt.1)go to 500
      do 410 i=m,lx
      if(x(i).gt.ritz+err)go to 420
      if(nu(i).lt.0)go to 500
  410 continue
! set up new fixed interval
  420 xl=ritz-tol
      if(xl.le.x(lfp))go to 450
! remove points in interval (xl,ritz)
      do 430 idummy=1,lx
      if(x(k-1).lt.xl)go to 440
      k=k-1
  430 continue
  440 x(k)=xl
      call ea15dd(jlan,alfa,beta,xl,del(k),nu(k),dr,nur)
      go to 455
  450 k=lfp
  455 nu(k)=-nu(k)
      xr=ritz+tol
! remove points in interval (ritz,xr)
      i=m
      do 460 m=i,lx
      if(x(m).gt.xr)go to 470
      if(nu(m).lt.0)go to 473
  460 continue
      m=lx+1
  470 m=m-1
      x(m)=xr
      call ea15dd(jlan-1,alfa,beta,x(m),del(m),nu(m),dr,nur)
  473 nfix=nfix+1
      if(neig.ge.leig)go to 920
      if(el.ge.er)go to 474
      if(ritz.lt.el .or. ritz.gt.er)go to 478
  474 i=neig
      if(neig.le.0)go to 477
      do 475 j=1,neig
      if(eig(i).lt.ritz)go to 477
      eig(i+1)=eig(i)
      jeig(1,i+1)=jeig(1,i)
      jeig(2,i+1)=jeig(2,i)
      i=i-1
  475 continue
  477 eig(i+1)=ritz
      jeig(1,i+1)=jlan
      jeig(2,i+1)=match
      neig=neig+1
  478 if(m.gt.k+2)go to 505
      go to 930
! extrapolate to estimate position of eigenvalue
  480 ritz=ritz+ritz-pole
      ritz=max(ritz,x(m-1)+(x(m)-x(m-1))*0.01)
      ritz=min(ritz,x(m)-(x(m)-x(m-1))*0.01)
! insert new point
  490 if(k.ge.m-3)go to 930
      m=m-1
      x(m-2)=x(m-1)
      del(m-2)=del(m-1)
      nu(m-2)=nu(m-1)
      x(m-1)=x(m)
      del(m-1)=del(m)
      nu(m-1)=nu(m)
      x(m)=ritz
      x(k)=ritz
      call ea15dd(jlan,alfa,beta,ritz,del(k),nu(k),dr,nur)
      del(m)=dr
      nu(m)=nur
      go to 530
! interval is acceptable. advance by one point
  500 ne=0
  502 m=m+1
  505 k=k+1
      x(k)=x(m)
      del(k)=alf-x(k)-be2/del(m)
      nu(k)=nu(m)
      if(del(k))510,520,530
  510 nu(k)=nu(k)+isign(1,nu(k))
      go to 530
  520 del(k)=drelpr*anorm
  530 go to 220
!
! scan to find the maximum number of ritz values (maxrz) in a gap
!     between fixed intervals, the number of gaps containing candidates
!     that are being watched (ncand) and the number of fixed intervals
!     (nfix).
  600 kx=k
      ncand=0
      maxrzo=maxrz
      maxrz=0
      nfix=0
      lfp=2
      k=k-1
      do 710 i=1,k
      if(nu(i).gt.0)go to 710
      nfix=nfix+1
      nritz=iabs(nu(i))-iabs(nu(lfp))
      maxrz=max0(maxrz,nritz)
      if(nritz.gt.0 .and. i.gt.lfp+1)ncand=ncand+1
      do 680 lfp=i,kx
      if(nu(lfp).gt.0)go to 710
  680 continue
      lfp=kx
  710 continue
      nritz=iabs(nu(k))-iabs(nu(lfp))
      maxrz=max0(maxrz,nritz)
      if(nritz.gt.0 .and. k.gt.lfp+1)ncand=ncand+1
      if(bet.le.en*anorm*drelpr)go to 960
      if(maxrz.gt.1 .or. ncand.gt.0)go to 910
      if(nfix.gt.0)go to 915
  910 continue
!
! normal returns
  912 iflag=1
      go to 1000
  915 iflag=0
      go to 1000
!
! error returns
  920 iflag=2
      if(lp.gt.0)write(lp,925)'error return 2 from lancz1. leig=',leig
  925 format( a , i6)
      go to 1000
  930 iflag=3
      if(lp.gt.0)write(lp,935)'error return 3 from lancz1. lx=',lx
  935 format( a , i6)
      go to 1000
  940 iflag=4
      if(lp.gt.0)write(lp,945)'error return 4 from lancz1. lalfa=',lalfa
  945 format( a , i6)
      go to 1000
  950 iflag=5
      if(lp.gt.0)write(lp,955)'error return 5 from lancz1. n,lx,lalfa=',&
                               n,lx,lalfa
  955 format( a , 3i6)
      go to 1000
  960 iflag=6
      if(lp.gt.0)write(lp,965)'error return 6 from lancz1. premature &
         &termination.',neig,'eigenvalues found'
  965 format( a , i6 , a)
      go to 1000
  970 if(lp.gt.0)write(lp,975)'error return 7 from lancz1. &
                             & on entry iflag is ',iflag
  975 format( a , i4)
      iflag=7
 1000 oldel=el
      older=er
      return
      end subroutine lancz1
      subroutine ea15cd(j,lalfa,alfa,beta,e,enorm,err,bnd,match)
! given a tridiagonal matrix t generated by the lanczos process applied
!     to a given symmetric matrix a and an approximate eigenvalue of t,
!     find a bound for its error as an eigenvalue of a.
! j is the index of the lanczos step. it is not altered.
! lalfa must be set to the lengths of arrays alfa and beta. it is not
!     altered.
! alfa(i),i=1,j must be set to the diagonal elements of the lanczos
!     tridiagonal matrix t. they are not altered.
! beta(i),i=2,j+1 must be set to the off-diagonal elements of the