![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DLACN2( N, V, X, ISGN, EST, KASE, ISAVE ) 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( * ), ISAVE( 3 ) 14: DOUBLE PRECISION V( * ), X( * ) 15: * .. 16: * 17: * Purpose 18: * ======= 19: * 20: * DLACN2 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 DLACN2 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 ISAVE(1) = 3, EST should be 44: * unchanged from the previous call to DLACN2. 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 DLACN2, 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 DLACN2, KASE will again be 0. 52: * 53: * ISAVE (input/output) INTEGER array, dimension (3) 54: * ISAVE is used to save variables between calls to DLACN2 55: * 56: * Further Details 57: * ======= ======= 58: * 59: * Contributed by Nick Higham, University of Manchester. 60: * Originally named SONEST, dated March 16, 1988. 61: * 62: * Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of 63: * a real or complex matrix, with applications to condition estimation", 64: * ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988. 65: * 66: * This is a thread safe version of DLACON, which uses the array ISAVE 67: * in place of a SAVE statement, as follows: 68: * 69: * DLACON DLACN2 70: * JUMP ISAVE(1) 71: * J ISAVE(2) 72: * ITER ISAVE(3) 73: * 74: * ===================================================================== 75: * 76: * .. Parameters .. 77: INTEGER ITMAX 78: PARAMETER ( ITMAX = 5 ) 79: DOUBLE PRECISION ZERO, ONE, TWO 80: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 ) 81: * .. 82: * .. Local Scalars .. 83: INTEGER I, JLAST 84: DOUBLE PRECISION ALTSGN, ESTOLD, TEMP 85: * .. 86: * .. External Functions .. 87: INTEGER IDAMAX 88: DOUBLE PRECISION DASUM 89: EXTERNAL IDAMAX, DASUM 90: * .. 91: * .. External Subroutines .. 92: EXTERNAL DCOPY 93: * .. 94: * .. Intrinsic Functions .. 95: INTRINSIC ABS, DBLE, NINT, SIGN 96: * .. 97: * .. Executable Statements .. 98: * 99: IF( KASE.EQ.0 ) THEN 100: DO 10 I = 1, N 101: X( I ) = ONE / DBLE( N ) 102: 10 CONTINUE 103: KASE = 1 104: ISAVE( 1 ) = 1 105: RETURN 106: END IF 107: * 108: GO TO ( 20, 40, 70, 110, 140 )ISAVE( 1 ) 109: * 110: * ................ ENTRY (ISAVE( 1 ) = 1) 111: * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X. 112: * 113: 20 CONTINUE 114: IF( N.EQ.1 ) THEN 115: V( 1 ) = X( 1 ) 116: EST = ABS( V( 1 ) ) 117: * ... QUIT 118: GO TO 150 119: END IF 120: EST = DASUM( N, X, 1 ) 121: * 122: DO 30 I = 1, N 123: X( I ) = SIGN( ONE, X( I ) ) 124: ISGN( I ) = NINT( X( I ) ) 125: 30 CONTINUE 126: KASE = 2 127: ISAVE( 1 ) = 2 128: RETURN 129: * 130: * ................ ENTRY (ISAVE( 1 ) = 2) 131: * FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. 132: * 133: 40 CONTINUE 134: ISAVE( 2 ) = IDAMAX( N, X, 1 ) 135: ISAVE( 3 ) = 2 136: * 137: * MAIN LOOP - ITERATIONS 2,3,...,ITMAX. 138: * 139: 50 CONTINUE 140: DO 60 I = 1, N 141: X( I ) = ZERO 142: 60 CONTINUE 143: X( ISAVE( 2 ) ) = ONE 144: KASE = 1 145: ISAVE( 1 ) = 3 146: RETURN 147: * 148: * ................ ENTRY (ISAVE( 1 ) = 3) 149: * X HAS BEEN OVERWRITTEN BY A*X. 150: * 151: 70 CONTINUE 152: CALL DCOPY( N, X, 1, V, 1 ) 153: ESTOLD = EST 154: EST = DASUM( N, V, 1 ) 155: DO 80 I = 1, N 156: IF( NINT( SIGN( ONE, X( I ) ) ).NE.ISGN( I ) ) 157: $ GO TO 90 158: 80 CONTINUE 159: * REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED. 160: GO TO 120 161: * 162: 90 CONTINUE 163: * TEST FOR CYCLING. 164: IF( EST.LE.ESTOLD ) 165: $ GO TO 120 166: * 167: DO 100 I = 1, N 168: X( I ) = SIGN( ONE, X( I ) ) 169: ISGN( I ) = NINT( X( I ) ) 170: 100 CONTINUE 171: KASE = 2 172: ISAVE( 1 ) = 4 173: RETURN 174: * 175: * ................ ENTRY (ISAVE( 1 ) = 4) 176: * X HAS BEEN OVERWRITTEN BY TRANSPOSE(A)*X. 177: * 178: 110 CONTINUE 179: JLAST = ISAVE( 2 ) 180: ISAVE( 2 ) = IDAMAX( N, X, 1 ) 181: IF( ( X( JLAST ).NE.ABS( X( ISAVE( 2 ) ) ) ) .AND. 182: $ ( ISAVE( 3 ).LT.ITMAX ) ) THEN 183: ISAVE( 3 ) = ISAVE( 3 ) + 1 184: GO TO 50 185: END IF 186: * 187: * ITERATION COMPLETE. FINAL STAGE. 188: * 189: 120 CONTINUE 190: ALTSGN = ONE 191: DO 130 I = 1, N 192: X( I ) = ALTSGN*( ONE+DBLE( I-1 ) / DBLE( N-1 ) ) 193: ALTSGN = -ALTSGN 194: 130 CONTINUE 195: KASE = 1 196: ISAVE( 1 ) = 5 197: RETURN 198: * 199: * ................ ENTRY (ISAVE( 1 ) = 5) 200: * X HAS BEEN OVERWRITTEN BY A*X. 201: * 202: 140 CONTINUE 203: TEMP = TWO*( DASUM( N, X, 1 ) / DBLE( 3*N ) ) 204: IF( TEMP.GT.EST ) THEN 205: CALL DCOPY( N, X, 1, V, 1 ) 206: EST = TEMP 207: END IF 208: * 209: 150 CONTINUE 210: KASE = 0 211: RETURN 212: * 213: * End of DLACN2 214: * 215: END