zlacon.f

Go to the documentation of this file.

      SUBROUTINE zlacon( N, V, X, EST, KASE )
*
*  -- LAPACK auxiliary routine (version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     June 30, 1999
*
*     .. Scalar Arguments ..
      INTEGER            KASE, N
      DOUBLE PRECISION   EST
*     ..
*     .. Array Arguments ..
      COMPLEX*16         V( N ), X( N )
*     ..
*
*  Purpose
*  =======
*
*  ZLACON estimates the 1-norm of a square, complex matrix A.
*  Reverse communication is used for evaluating matrix-vector products.
*
*  Arguments
*  =========
*
*  N      (input) INTEGER
*         The order of the matrix.  N >= 1.
*
*  V      (workspace) COMPLEX*16 array, dimension (N)
*         On the final return, V = A*W,  where  EST = norm(V)/norm(W)
*         (W is not returned).
*
*  X      (input/output) COMPLEX*16 array, dimension (N)
*         On an intermediate return, X should be overwritten by
*               A * X,   if KASE=1,
*               A' * X,  if KASE=2,
*         where A' is the conjugate transpose of A, and ZLACON must be
*         re-called with all the other parameters unchanged.
*
*  EST    (output) DOUBLE PRECISION
*         An estimate (a lower bound) for norm(A).
*
*  KASE   (input/output) INTEGER
*         On the initial call to ZLACON, KASE should be 0.
*         On an intermediate return, KASE will be 1 or 2, indicating
*         whether X should be overwritten by A * X  or A' * X.
*         On the final return from ZLACON, KASE will again be 0.
*
*  Further Details
*  ======= =======
*
*  Contributed by Nick Higham, University of Manchester.
*  Originally named CONEST, dated March 16, 1988.
*
*  Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of
*  a real or complex matrix, with applications to condition estimation",
*  ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
*
*  Last modified:  April, 1999
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            ITMAX
      parameter( itmax = 5 )
      DOUBLE PRECISION   ONE, TWO
      parameter( one = 1.0d0, two = 2.0d0 )
      COMPLEX*16         CZERO, CONE
      parameter( czero = ( 0.0d0, 0.0d0 ),
     $                   cone = ( 1.0d0, 0.0d0 ) )
*     ..
*     .. Local Scalars ..
      INTEGER            I, ITER, J, JLAST, JUMP
      DOUBLE PRECISION   ABSXI, ALTSGN, ESTOLD, SAFMIN, TEMP
*     ..
*     .. External Functions ..
      INTEGER            IZMAX1
      DOUBLE PRECISION   DLAMCH, DZSUM1
      EXTERNAL           izmax1, dlamch, dzsum1
*     ..
*     .. External Subroutines ..
      EXTERNAL           zcopy
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, dcmplx, dimag
*     ..
*     .. Save statement ..
      SAVE
*     ..
*     .. Executable Statements ..
*
      safmin = dlamch( 'Safe minimum' )
      IF( kase.EQ.0 ) THEN
         DO 10 i = 1, n
            x( i ) = dcmplx( one / dble( n ) )
   10    CONTINUE
         kase = 1
         jump = 1
         RETURN
      END IF
*
      GO TO ( 20, 40, 70, 90, 120 )jump
*
*     ................ ENTRY   (JUMP = 1)
*     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X.
*
   20 CONTINUE
      IF( n.EQ.1 ) THEN
         v( 1 ) = x( 1 )
         est = abs( v( 1 ) )
*        ... QUIT
         GO TO 130
      END IF
      est = dzsum1( n, x, 1 )
*
      DO 30 i = 1, n
         absxi = abs( x( i ) )
         IF( absxi.GT.safmin ) THEN
            x( i ) = dcmplx( dble( x( i ) ) / absxi,
     $               dimag( x( i ) ) / absxi )
         ELSE
            x( i ) = cone
         END IF
   30 CONTINUE
      kase = 2
      jump = 2
      RETURN
*
*     ................ ENTRY   (JUMP = 2)
*     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY ZTRANS(A)*X.
*
   40 CONTINUE
      j = izmax1( n, x, 1 )
      iter = 2
*
*     MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
*
   50 CONTINUE
      DO 60 i = 1, n
         x( i ) = czero
   60 CONTINUE
      x( j ) = cone
      kase = 1
      jump = 3
      RETURN
*
*     ................ ENTRY   (JUMP = 3)
*     X HAS BEEN OVERWRITTEN BY A*X.
*
   70 CONTINUE
      CALL zcopy( n, x, 1, v, 1 )
      estold = est
      est = dzsum1( n, v, 1 )
*
*     TEST FOR CYCLING.
      IF( est.LE.estold )
     $   GO TO 100
*
      DO 80 i = 1, n
         absxi = abs( x( i ) )
         IF( absxi.GT.safmin ) THEN
            x( i ) = dcmplx( dble( x( i ) ) / absxi,
     $               dimag( x( i ) ) / absxi )
         ELSE
            x( i ) = cone
         END IF
   80 CONTINUE
      kase = 2
      jump = 4
      RETURN
*
*     ................ ENTRY   (JUMP = 4)
*     X HAS BEEN OVERWRITTEN BY ZTRANS(A)*X.
*
   90 CONTINUE
      jlast = j
      j = izmax1( n, x, 1 )
      IF( ( abs( x( jlast ) ).NE.abs( x( j ) ) ) .AND.
     $    ( iter.LT.itmax ) ) THEN
         iter = iter + 1
         GO TO 50
      END IF
*
*     ITERATION COMPLETE.  FINAL STAGE.
*
  100 CONTINUE
      altsgn = one
      DO 110 i = 1, n
         x( i ) = dcmplx( altsgn*( one+dble( i-1 ) / dble( n-1 ) ) )
         altsgn = -altsgn
  110 CONTINUE
      kase = 1
      jump = 5
      RETURN
*
*     ................ ENTRY   (JUMP = 5)
*     X HAS BEEN OVERWRITTEN BY A*X.
*
  120 CONTINUE
      temp = two*( dzsum1( n, x, 1 ) / dble( 3*n ) )
      IF( temp.GT.est ) THEN
         CALL zcopy( n, x, 1, v, 1 )
         est = temp
      END IF
*
  130 CONTINUE
      kase = 0
      RETURN
*
*     End of ZLACON
*
      END