|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
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
1.7.3 
