symmlq.f

一个画有向图的程序。里面含有力导引画图算法等多个经典算法。

*     y and beta will be zero or very small if b is an eigenvector.

      call daxpy ( n, (- shift), v, 1, y, 1 )
      alfa   = ddot  ( n, v, 1, y, 1 )
      call daxpy ( n, (- alfa / beta1), r1, 1, y, 1 )

*     Make sure  r2  will be orthogonal to the first  v.

      z      = ddot  ( n, v, 1, y, 1 )
      s      = ddot  ( n, v, 1, v, 1 )
      call daxpy ( n, (- z / s), v, 1, y, 1 )

      call dcopy ( n, y, 1, r2, 1 )
      if ( precon ) call msolve( n, r2, y, A, vwsqrt, work )
      oldb   = beta1
      beta   = ddot  ( n, r2, 1, y, 1 )

      if (beta .lt. zero) then
         istop = 8
         go to 900
      end if

*     Cause termination (later) if beta is essentially zero.

      beta   = sqrt( beta )
      if (beta .le. eps) then
         istop = -1
      end if

*     See if the local reorthogonalization achieved anything.

      denom  = sqrt( s ) * dnrm2( n, r2, 1 )  +  eps
      s      = z / denom
      t      = ddot  ( n, v, 1, r2, 1 ) / denom

c     if (nout .gt. 0  .and.  goodb) then
c        write(nout, 1100) beta1, alfa, s, t
c     end if

*     Initialize other quantities.

      cgnorm = beta1
      gbar   = alfa
      dbar   = beta
      rhs1   = beta1
      rhs2   = zero
      bstep  = zero
      snprod = one
      tnorm  = alfa**2
      ynorm2 = zero
      gmax   = abs( alfa )
      gmin   = gmax

      if ( goodb ) then
         do 200 i = 1, n
            w(i)  = zero
  200    continue
      else
         call dcopy ( n, v, 1, w, 1 )
      end if

*     ------------------------------------------------------------------
*     Main iteration loop.
*     ------------------------------------------------------------------

*     Estimate various norms and test for convergence.

  300 anorm  = sqrt( tnorm  )
      ynorm  = sqrt( ynorm2 )
      epsa   = anorm * eps
      epsx   = anorm * ynorm * eps
      epsr   = anorm * ynorm * rtol
      diag   = gbar
      if (diag .eq. zero) diag = epsa

      lqnorm = sqrt( rhs1**2 + rhs2**2 )
      qrnorm = snprod * beta1
      cgnorm = qrnorm * beta / abs( diag )

*     Estimate  cond(A).
*     In this version we look at the diagonals of  L  in the
*     factorization of the tridiagonal matrix,  T = L*Q.
*     Sometimes, T(k) can be misleadingly ill-conditioned when
*     T(k+1) is not, so we must be careful not to overestimate acond.

      if (lqnorm .le. cgnorm) then
         acond  = gmax / gmin
      else
         denom  = min( gmin, abs( diag ) )
         acond  = gmax / denom
      end if

*     See if any of the stopping criteria are satisfied.
*     In rare cases, istop is already -1 from above (Abar = const * I).

      if (istop .eq. 0) then
         if (itn    .ge. itnlim ) istop = 5
c         if (acond  .ge. 0.1/eps) istop = 4
         if (epsx   .ge. beta1  ) istop = 3
         if (cgnorm .le. epsx   ) istop = 2
         if (cgnorm .le. epsr   ) istop = 1
	 if ((istop .ne. 4).and.(ynorm .ge. normxlim)
     1	      .and. (normxlim .ne. 0.0)) istop = 9
      end if
      if (itn .lt. itnmin) istop = 0
*     ==================================================================

*     See if it is time to print something.

      if (nout .le.  0)          go to 600
      if (n    .le. 40)          go to 400
      if (itn  .le. 10)          go to 400
      if (itn  .ge. itnlim - 10) go to 400
      if (mod(itn,10)  .eq.   0) go to 400
      if (cgnorm .le. 10.0*epsx) go to 400
      if (cgnorm .le. 10.0*epsr) go to 400
      if (acond  .ge. 0.01/eps ) go to 400
      if (istop  .ne. 0)         go to 400
      go to 600

*     Print a line for this iteration.

  400 zbar   = rhs1 / diag
      z      = (snprod * zbar  +  bstep) / beta1
c     x1lq   = x(1)  +  b1 * bstep / beta1
c     x1cg   = x(1)  +  w(1) * zbar  +  b1 * z

c     if (    itn     .eq. 0) write(nout, 1200)
c     write(nout, 1300) itn, x1cg, cgnorm, bstep/beta1, anorm, acond
c     if (mod(itn,10) .eq. 0) write(nout, 1500)
*     ==================================================================


*     Obtain the current Lanczos vector  v = (1 / beta)*y
*     and set up  y  for the next iteration.

  600 if (istop .ne. 0) go to 800
      s      = one / beta

      do 620 i = 1, n
         v(i)  = s * y(i)
  620 continue

      call aprod ( n, v, y, A, vwsqrt, work, orthlist )
      call daxpy ( n, (- shift), v, 1, y, 1 )
      call daxpy ( n, (- beta / oldb), r1, 1, y, 1 )
      alfa   = ddot( n, v, 1, y, 1 )
      tnorm  = tnorm  +  alfa**2  +  two * beta**2
      call daxpy ( n, (- alfa / beta), r2, 1, y, 1 )
      call dcopy ( n, r2, 1, r1, 1 )
      call dcopy ( n, y, 1, r2, 1 )
      if ( precon ) call msolve( n, r2, y, A, vwsqrt, work )
      oldb   = beta
      beta   = ddot  ( n, r2, 1, y, 1 )
      if (beta .lt. zero) then
         istop = 6
         go to 800
      end if
      beta   = sqrt( beta )

*     Compute the next plane rotation for  Q.

      gamma  = sqrt( gbar**2 + oldb**2 )
      cs     = gbar / gamma
      sn     = oldb / gamma
      delta  = cs * dbar  +  sn * alfa
      gbar   = sn * dbar  -  cs * alfa
      epsln  = sn * beta
      dbar   =            -  cs * beta

*     Update  x.

      z      = rhs1 / gamma
      s      = z * cs
      t      = z * sn

      do 700 i = 1, n
         x(i)  = (w(i) * s   +   v(i) * t)  +  x(i)
         w(i)  =  w(i) * sn  -   v(i) * cs
  700 continue

*     Accumulate the step along the direction  b,
*     and go round again.

      bstep  = snprod * cs * z  +  bstep
      snprod = snprod * sn
      gmax   = max( gmax, gamma )
      gmin   = min( gmin, gamma )
      ynorm2 = z**2  +  ynorm2
      rhs1   = rhs2  -  delta * z
      rhs2   =       -  epsln * z
      itn    = itn   +  1
      go to 300

*     ------------------------------------------------------------------
*     End of main iteration loop.
*     ------------------------------------------------------------------

*     Move to the CG point if it seems better.
*     In this version of SYMMLQ, the convergence tests involve
*     only cgnorm, so we're unlikely to stop at an LQ point,
*     EXCEPT if the iteration limit interferes.

  800 if (cgnorm .le. lqnorm) then
         zbar   = rhs1 / diag
         bstep  = snprod * zbar  +  bstep
         ynorm  = sqrt( ynorm2  +  zbar**2 )
         rnorm  = cgnorm
         call daxpy ( n, zbar, w, 1, x, 1 )
      else
         rnorm  = lqnorm
      end if

      if ( goodb ) then

*        Add the step along  b.

         bstep  = bstep / beta1
         call dcopy ( n, b, 1, y, 1 )
         if ( precon ) call msolve( n, b, y, A, vwsqrt, work )
         call daxpy ( n, bstep, y, 1, x, 1 )
      end if

*     ==================================================================
*     Display final status.
*     ==================================================================
  900 continue
c     if (nout  .gt. 0) then
c        write(nout, 2000) exit, istop, itn,
c    $                     exit, anorm, acond,
c    $                     exit, rnorm, ynorm
c        write(nout, 3000) exit, msg(istop)
c     end if

      return

*     ------------------------------------------------------------------
c 1000 format(// 1p,    a, 5x, 'Solution of symmetric   Ax = b'
c     $       / ' n      =', i7, 5x, 'checka =', l4, 12x,
c     $          'goodb  =', l4, 7x, 'precon =', l4
c     $       / ' itnlim =', i7, 5x, 'rtol   =', e11.2, 5x,
c     $          'shift  =', e23.14)
c 1100 format(/ 1p, ' beta1  =', e10.2, 3x, 'alpha1 =', e10.2
c     $       / ' (v1,v2) before and after ', e14.2
c     $       / ' local reorthogonalization', e14.2)
c 1200 format(// 5x, 'itn', 7x, 'x1(cg)', 10x,
c     $         'norm(r)', 5x, 'bstep', 7x, 'norm(A)', 3X, 'cond(A)')
c 1300 format(1p, i8, e19.10, e11.2, e14.5, 2e10.2)
c 1500 format(1x)
c 2000 format(/ 1p, a, 6x, 'istop =', i3,   15x, 'itn   =', i8
c     $       /     a, 6x, 'anorm =', e12.4, 6x, 'acond =', e12.4
c     $       /     a, 6x, 'rnorm =', e12.4, 6x, 'ynorm =', e12.4)
c 3000 format(      a, 6x, a )
*     ------------------------------------------------------------------
*     end of SYMMLQ
      end