![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DLACON( N, V, X, ISGN, EST, KASE ) 2: * 3: * -- LAPACK auxiliary routine (version 3.2) -- 4: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 6: * November 2006 7: * 8: * .. Scalar Arguments .. 9: INTEGER KASE, N 10: DOUBLE PRECISION EST 11: * .. 12: * .. Array Arguments .. 13: INTEGER ISGN( * ) 14: DOUBLE PRECISION V( * ), X( * ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * DLACON estimates the 1-norm of a square, real matrix A. 21: * Reverse communication is used for evaluating matrix-vector products. 22: * 23: * Arguments 24: * ========= 25: * 26: * N (input) INTEGER 27: * The order of the matrix. N >= 1. 28: * 29: * V (workspace) DOUBLE PRECISION array, dimension (N) 30: * On the final return, V = A*W, where EST = norm(V)/norm(W) 31: * (W is not returned). 32: * 33: * X (input/output) DOUBLE PRECISION array, dimension (N) 34: * On an intermediate return, X should be overwritten by 35: * A * X, if KASE=1, 36: * A' * X, if KASE=2, 37: * and DLACON must be re-called with all the other parameters 38: * unchanged. 39: * 40: * ISGN (workspace) INTEGER array, dimension (N) 41: * 42: * EST (input/output) DOUBLE PRECISION 43: * On entry with KASE = 1 or 2 and JUMP = 3, EST should be 44: * unchanged from the previous call to DLACON. 45: * On exit, EST is an estimate (a lower bound) for norm(A). 46: * 47: * KASE (input/output) INTEGER 48: * On the initial call to DLACON, KASE should be 0. 49: * On an intermediate return, KASE will be 1 or 2, indicating 50: * whether X should be overwritten by A * X or A' * X. 51: * On the final return from DLACON, KASE will again be 0. 52: * 53: * Further Details 54: * ======= ======= 55: * 56: * Contributed by Nick Higham, University of Manchester. 57: * Originally named SONEST, dated March 16, 1988. 58: * 59: * Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of 60: * a real or complex matrix, with applications to condition estimation", 61: * ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988. 62: * 63: * ===================================================================== 64: * 65: * .. Parameters .. 66: INTEGER ITMAX 67: PARAMETER ( ITMAX = 5 ) 68: DOUBLE PRECISION ZERO, ONE, TWO 69: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 ) 70: * .. 71: * .. Local Scalars .. 72: INTEGER I, ITER, J, JLAST, JUMP 73: DOUBLE PRECISION ALTSGN, ESTOLD, TEMP 74: * .. 75: * .. External Functions .. 76: INTEGER IDAMAX 77: DOUBLE PRECISION DASUM 78: EXTERNAL IDAMAX, DASUM 79: * .. 80: * .. External Subroutines .. 81: EXTERNAL DCOPY 82: * .. 83: * .. Intrinsic Functions .. 84: INTRINSIC ABS, DBLE, NINT, SIGN 85: * .. 86: * .. Save statement .. 87: SAVE 88: * .. 89: * .. Executable Statements .. 90: * 91: IF( KASE.EQ.0 ) THEN 92: DO 10 I = 1, N 93: X( I ) = ONE / DBLE( N ) 94: 10 CONTINUE 95: KASE = 1 96: JUMP = 1 97: RETURN 98: END IF 99: * 100: GO TO ( 20, 40, 70, 110, 140 )JUMP 101: * 102: * ................ ENTRY (JUMP = 1) 103: * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X. 104: * 105: 20 CONTINUE 106: IF( N.EQ.1 ) THEN 107: V( 1 ) = X( 1 ) 108: EST = ABS( V( 1 ) ) 109: * ... QUIT 110: GO TO 150 111: END IF 112: EST = DASUM( N, X, 1 ) 113: * 114: DO 30 I = 1, N 115: X( I ) = SIGN( ONE, X( I ) ) 116: ISGN( I ) = NINT( X( I ) ) 117: 30 CONTINUE 118: KASE = 2 119: JUMP = 2 120: RETURN 121: * 122: * ................ ENTRY (JUMP = 2) 123: * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. 124: * 125: 40 CONTINUE 126: J = IDAMAX( N, X, 1 ) 127: ITER = 2 128: * 129: * MAIN LOOP - ITERATIONS 2,3,...,ITMAX. 130: * 131: 50 CONTINUE 132: DO 60 I = 1, N 133: X( I ) = ZERO 134: 60 CONTINUE 135: X( J ) = ONE 136: KASE = 1 137: JUMP = 3 138: RETURN 139: * 140: * ................ ENTRY (JUMP = 3) 141: * X HAS BEEN OVERWRITTEN BY A*X. 142: * 143: 70 CONTINUE 144: CALL DCOPY( N, X, 1, V, 1 ) 145: ESTOLD = EST 146: EST = DASUM( N, V, 1 ) 147: DO 80 I = 1, N 148: IF( NINT( SIGN( ONE, X( I ) ) ).NE.ISGN( I ) ) 149: $ GO TO 90 150: 80 CONTINUE 151: * REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED. 152: GO TO 120 153: * 154: 90 CONTINUE 155: * TEST FOR CYCLING. 156: IF( EST.LE.ESTOLD ) 157: $ GO TO 120 158: * 159: DO 100 I = 1, N 160: X( I ) = SIGN( ONE, X( I ) ) 161: ISGN( I ) = NINT( X( I ) ) 162: 100 CONTINUE 163: KASE = 2 164: JUMP = 4 165: RETURN 166: * 167: * ................ ENTRY (JUMP = 4) 168: * X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. 169: * 170: 110 CONTINUE 171: JLAST = J 172: J = IDAMAX( N, X, 1 ) 173: IF( ( X( JLAST ).NE.ABS( X( J ) ) ) .AND. ( ITER.LT.ITMAX ) ) THEN 174: ITER = ITER + 1 175: GO TO 50 176: END IF 177: * 178: * ITERATION COMPLETE. FINAL STAGE. 179: * 180: 120 CONTINUE 181: ALTSGN = ONE 182: DO 130 I = 1, N 183: X( I ) = ALTSGN*( ONE+DBLE( I-1 ) / DBLE( N-1 ) ) 184: ALTSGN = -ALTSGN 185: 130 CONTINUE 186: KASE = 1 187: JUMP = 5 188: RETURN 189: * 190: * ................ ENTRY (JUMP = 5) 191: * X HAS BEEN OVERWRITTEN BY A*X. 192: * 193: 140 CONTINUE 194: TEMP = TWO*( DASUM( N, X, 1 ) / DBLE( 3*N ) ) 195: IF( TEMP.GT.EST ) THEN 196: CALL DCOPY( N, X, 1, V, 1 ) 197: EST = TEMP 198: END IF 199: * 200: 150 CONTINUE 201: KASE = 0 202: RETURN 203: * 204: * End of DLACON 205: * 206: END