LAPACK: SRC/zlacon.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZLACON( N, V, X, EST, KASE )
00002 *
00003 *  -- LAPACK auxiliary routine (version 3.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       INTEGER            KASE, N
00010       DOUBLE PRECISION   EST
00011 *     ..
00012 *     .. Array Arguments ..
00013       COMPLEX*16         V( N ), X( N )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  ZLACON estimates the 1-norm of a square, complex matrix A.
00020 *  Reverse communication is used for evaluating matrix-vector products.
00021 *
00022 *  Arguments
00023 *  =========
00024 *
00025 *  N      (input) INTEGER
00026 *         The order of the matrix.  N >= 1.
00027 *
00028 *  V      (workspace) COMPLEX*16 array, dimension (N)
00029 *         On the final return, V = A*W,  where  EST = norm(V)/norm(W)
00030 *         (W is not returned).
00031 *
00032 *  X      (input/output) COMPLEX*16 array, dimension (N)
00033 *         On an intermediate return, X should be overwritten by
00034 *               A * X,   if KASE=1,
00035 *               A' * X,  if KASE=2,
00036 *         where A' is the conjugate transpose of A, and ZLACON must be
00037 *         re-called with all the other parameters unchanged.
00038 *
00039 *  EST    (input/output) DOUBLE PRECISION
00040 *         On entry with KASE = 1 or 2 and JUMP = 3, EST should be
00041 *         unchanged from the previous call to ZLACON.
00042 *         On exit, EST is an estimate (a lower bound) for norm(A). 
00043 *
00044 *  KASE   (input/output) INTEGER
00045 *         On the initial call to ZLACON, KASE should be 0.
00046 *         On an intermediate return, KASE will be 1 or 2, indicating
00047 *         whether X should be overwritten by A * X  or A' * X.
00048 *         On the final return from ZLACON, KASE will again be 0.
00049 *
00050 *  Further Details
00051 *  ======= =======
00052 *
00053 *  Contributed by Nick Higham, University of Manchester.
00054 *  Originally named CONEST, dated March 16, 1988.
00055 *
00056 *  Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of
00057 *  a real or complex matrix, with applications to condition estimation",
00058 *  ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
00059 *
00060 *  Last modified:  April, 1999
00061 *
00062 *  =====================================================================
00063 *
00064 *     .. Parameters ..
00065       INTEGER            ITMAX
00066       PARAMETER          ( ITMAX = 5 )
00067       DOUBLE PRECISION   ONE, TWO
00068       PARAMETER          ( ONE = 1.0D0, TWO = 2.0D0 )
00069       COMPLEX*16         CZERO, CONE
00070       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
00071      $                   CONE = ( 1.0D0, 0.0D0 ) )
00072 *     ..
00073 *     .. Local Scalars ..
00074       INTEGER            I, ITER, J, JLAST, JUMP
00075       DOUBLE PRECISION   ABSXI, ALTSGN, ESTOLD, SAFMIN, TEMP
00076 *     ..
00077 *     .. External Functions ..
00078       INTEGER            IZMAX1
00079       DOUBLE PRECISION   DLAMCH, DZSUM1
00080       EXTERNAL           IZMAX1, DLAMCH, DZSUM1
00081 *     ..
00082 *     .. External Subroutines ..
00083       EXTERNAL           ZCOPY
00084 *     ..
00085 *     .. Intrinsic Functions ..
00086       INTRINSIC          ABS, DBLE, DCMPLX, DIMAG
00087 *     ..
00088 *     .. Save statement ..
00089       SAVE
00090 *     ..
00091 *     .. Executable Statements ..
00092 *
00093       SAFMIN = DLAMCH( 'Safe minimum' )
00094       IF( KASE.EQ.0 ) THEN
00095          DO 10 I = 1, N
00096             X( I ) = DCMPLX( ONE / DBLE( N ) )
00097    10    CONTINUE
00098          KASE = 1
00099          JUMP = 1
00100          RETURN
00101       END IF
00102 *
00103       GO TO ( 20, 40, 70, 90, 120 )JUMP
00104 *
00105 *     ................ ENTRY   (JUMP = 1)
00106 *     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY A*X.
00107 *
00108    20 CONTINUE
00109       IF( N.EQ.1 ) THEN
00110          V( 1 ) = X( 1 )
00111          EST = ABS( V( 1 ) )
00112 *        ... QUIT
00113          GO TO 130
00114       END IF
00115       EST = DZSUM1( N, X, 1 )
00116 *
00117       DO 30 I = 1, N
00118          ABSXI = ABS( X( I ) )
00119          IF( ABSXI.GT.SAFMIN ) THEN
00120             X( I ) = DCMPLX( DBLE( X( I ) ) / ABSXI,
00121      $               DIMAG( X( I ) ) / ABSXI )
00122          ELSE
00123             X( I ) = CONE
00124          END IF
00125    30 CONTINUE
00126       KASE = 2
00127       JUMP = 2
00128       RETURN
00129 *
00130 *     ................ ENTRY   (JUMP = 2)
00131 *     FIRST ITERATION.  X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
00132 *
00133    40 CONTINUE
00134       J = IZMAX1( N, X, 1 )
00135       ITER = 2
00136 *
00137 *     MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
00138 *
00139    50 CONTINUE
00140       DO 60 I = 1, N
00141          X( I ) = CZERO
00142    60 CONTINUE
00143       X( J ) = CONE
00144       KASE = 1
00145       JUMP = 3
00146       RETURN
00147 *
00148 *     ................ ENTRY   (JUMP = 3)
00149 *     X HAS BEEN OVERWRITTEN BY A*X.
00150 *
00151    70 CONTINUE
00152       CALL ZCOPY( N, X, 1, V, 1 )
00153       ESTOLD = EST
00154       EST = DZSUM1( N, V, 1 )
00155 *
00156 *     TEST FOR CYCLING.
00157       IF( EST.LE.ESTOLD )
00158      $   GO TO 100
00159 *
00160       DO 80 I = 1, N
00161          ABSXI = ABS( X( I ) )
00162          IF( ABSXI.GT.SAFMIN ) THEN
00163             X( I ) = DCMPLX( DBLE( X( I ) ) / ABSXI,
00164      $               DIMAG( X( I ) ) / ABSXI )
00165          ELSE
00166             X( I ) = CONE
00167          END IF
00168    80 CONTINUE
00169       KASE = 2
00170       JUMP = 4
00171       RETURN
00172 *
00173 *     ................ ENTRY   (JUMP = 4)
00174 *     X HAS BEEN OVERWRITTEN BY CTRANS(A)*X.
00175 *
00176    90 CONTINUE
00177       JLAST = J
00178       J = IZMAX1( N, X, 1 )
00179       IF( ( ABS( X( JLAST ) ).NE.ABS( X( J ) ) ) .AND.
00180      $    ( ITER.LT.ITMAX ) ) THEN
00181          ITER = ITER + 1
00182          GO TO 50
00183       END IF
00184 *
00185 *     ITERATION COMPLETE.  FINAL STAGE.
00186 *
00187   100 CONTINUE
00188       ALTSGN = ONE
00189       DO 110 I = 1, N
00190          X( I ) = DCMPLX( ALTSGN*( ONE+DBLE( I-1 ) / DBLE( N-1 ) ) )
00191          ALTSGN = -ALTSGN
00192   110 CONTINUE
00193       KASE = 1
00194       JUMP = 5
00195       RETURN
00196 *
00197 *     ................ ENTRY   (JUMP = 5)
00198 *     X HAS BEEN OVERWRITTEN BY A*X.
00199 *
00200   120 CONTINUE
00201       TEMP = TWO*( DZSUM1( N, X, 1 ) / DBLE( 3*N ) )
00202       IF( TEMP.GT.EST ) THEN
00203          CALL ZCOPY( N, X, 1, V, 1 )
00204          EST = TEMP
00205       END IF
00206 *
00207   130 CONTINUE
00208       KASE = 0
00209       RETURN
00210 *
00211 *     End of ZLACON
00212 *
00213       END