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