Annotation of rpl/lapack/lapack/dla_syrpvgrw.f, revision 1.1
1.1 ! bertrand 1: DOUBLE PRECISION FUNCTION DLA_SYRPVGRW( UPLO, N, INFO, A, LDA, AF,
! 2: $ LDAF, IPIV, WORK )
! 3: *
! 4: * -- LAPACK routine (version 3.2.2) --
! 5: * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
! 6: * -- Jason Riedy of Univ. of California Berkeley. --
! 7: * -- June 2010 --
! 8: *
! 9: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 10: * -- Univ. of California Berkeley and NAG Ltd. --
! 11: *
! 12: IMPLICIT NONE
! 13: * ..
! 14: * .. Scalar Arguments ..
! 15: CHARACTER*1 UPLO
! 16: INTEGER N, INFO, LDA, LDAF
! 17: * ..
! 18: * .. Array Arguments ..
! 19: INTEGER IPIV( * )
! 20: DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), WORK( * )
! 21: * ..
! 22: *
! 23: * Purpose
! 24: * =======
! 25: *
! 26: * DLA_SYRPVGRW computes the reciprocal pivot growth factor
! 27: * norm(A)/norm(U). The "max absolute element" norm is used. If this is
! 28: * much less than 1, the stability of the LU factorization of the
! 29: * (equilibrated) matrix A could be poor. This also means that the
! 30: * solution X, estimated condition numbers, and error bounds could be
! 31: * unreliable.
! 32: *
! 33: * Arguments
! 34: * =========
! 35: *
! 36: * UPLO (input) CHARACTER*1
! 37: * = 'U': Upper triangle of A is stored;
! 38: * = 'L': Lower triangle of A is stored.
! 39: *
! 40: * N (input) INTEGER
! 41: * The number of linear equations, i.e., the order of the
! 42: * matrix A. N >= 0.
! 43: *
! 44: * INFO (input) INTEGER
! 45: * The value of INFO returned from DSYTRF, .i.e., the pivot in
! 46: * column INFO is exactly 0.
! 47: *
! 48: * NCOLS (input) INTEGER
! 49: * The number of columns of the matrix A. NCOLS >= 0.
! 50: *
! 51: * A (input) DOUBLE PRECISION array, dimension (LDA,N)
! 52: * On entry, the N-by-N matrix A.
! 53: *
! 54: * LDA (input) INTEGER
! 55: * The leading dimension of the array A. LDA >= max(1,N).
! 56: *
! 57: * AF (input) DOUBLE PRECISION array, dimension (LDAF,N)
! 58: * The block diagonal matrix D and the multipliers used to
! 59: * obtain the factor U or L as computed by DSYTRF.
! 60: *
! 61: * LDAF (input) INTEGER
! 62: * The leading dimension of the array AF. LDAF >= max(1,N).
! 63: *
! 64: * IPIV (input) INTEGER array, dimension (N)
! 65: * Details of the interchanges and the block structure of D
! 66: * as determined by DSYTRF.
! 67: *
! 68: * WORK (input) DOUBLE PRECISION array, dimension (2*N)
! 69: *
! 70: * =====================================================================
! 71: *
! 72: * .. Local Scalars ..
! 73: INTEGER NCOLS, I, J, K, KP
! 74: DOUBLE PRECISION AMAX, UMAX, RPVGRW, TMP
! 75: LOGICAL UPPER
! 76: * ..
! 77: * .. Intrinsic Functions ..
! 78: INTRINSIC ABS, MAX, MIN
! 79: * ..
! 80: * .. External Functions ..
! 81: EXTERNAL LSAME, DLASET
! 82: LOGICAL LSAME
! 83: * ..
! 84: * .. Executable Statements ..
! 85: *
! 86: UPPER = LSAME( 'Upper', UPLO )
! 87: IF ( INFO.EQ.0 ) THEN
! 88: IF ( UPPER ) THEN
! 89: NCOLS = 1
! 90: ELSE
! 91: NCOLS = N
! 92: END IF
! 93: ELSE
! 94: NCOLS = INFO
! 95: END IF
! 96:
! 97: RPVGRW = 1.0D+0
! 98: DO I = 1, 2*N
! 99: WORK( I ) = 0.0D+0
! 100: END DO
! 101: *
! 102: * Find the max magnitude entry of each column of A. Compute the max
! 103: * for all N columns so we can apply the pivot permutation while
! 104: * looping below. Assume a full factorization is the common case.
! 105: *
! 106: IF ( UPPER ) THEN
! 107: DO J = 1, N
! 108: DO I = 1, J
! 109: WORK( N+I ) = MAX( ABS( A( I, J ) ), WORK( N+I ) )
! 110: WORK( N+J ) = MAX( ABS( A( I, J ) ), WORK( N+J ) )
! 111: END DO
! 112: END DO
! 113: ELSE
! 114: DO J = 1, N
! 115: DO I = J, N
! 116: WORK( N+I ) = MAX( ABS( A( I, J ) ), WORK( N+I ) )
! 117: WORK( N+J ) = MAX( ABS( A( I, J ) ), WORK( N+J ) )
! 118: END DO
! 119: END DO
! 120: END IF
! 121: *
! 122: * Now find the max magnitude entry of each column of U or L. Also
! 123: * permute the magnitudes of A above so they're in the same order as
! 124: * the factor.
! 125: *
! 126: * The iteration orders and permutations were copied from dsytrs.
! 127: * Calls to SSWAP would be severe overkill.
! 128: *
! 129: IF ( UPPER ) THEN
! 130: K = N
! 131: DO WHILE ( K .LT. NCOLS .AND. K.GT.0 )
! 132: IF ( IPIV( K ).GT.0 ) THEN
! 133: ! 1x1 pivot
! 134: KP = IPIV( K )
! 135: IF ( KP .NE. K ) THEN
! 136: TMP = WORK( N+K )
! 137: WORK( N+K ) = WORK( N+KP )
! 138: WORK( N+KP ) = TMP
! 139: END IF
! 140: DO I = 1, K
! 141: WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
! 142: END DO
! 143: K = K - 1
! 144: ELSE
! 145: ! 2x2 pivot
! 146: KP = -IPIV( K )
! 147: TMP = WORK( N+K-1 )
! 148: WORK( N+K-1 ) = WORK( N+KP )
! 149: WORK( N+KP ) = TMP
! 150: DO I = 1, K-1
! 151: WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
! 152: WORK( K-1 ) = MAX( ABS( AF( I, K-1 ) ), WORK( K-1 ) )
! 153: END DO
! 154: WORK( K ) = MAX( ABS( AF( K, K ) ), WORK( K ) )
! 155: K = K - 2
! 156: END IF
! 157: END DO
! 158: K = NCOLS
! 159: DO WHILE ( K .LE. N )
! 160: IF ( IPIV( K ).GT.0 ) THEN
! 161: KP = IPIV( K )
! 162: IF ( KP .NE. K ) THEN
! 163: TMP = WORK( N+K )
! 164: WORK( N+K ) = WORK( N+KP )
! 165: WORK( N+KP ) = TMP
! 166: END IF
! 167: K = K + 1
! 168: ELSE
! 169: KP = -IPIV( K )
! 170: TMP = WORK( N+K )
! 171: WORK( N+K ) = WORK( N+KP )
! 172: WORK( N+KP ) = TMP
! 173: K = K + 2
! 174: END IF
! 175: END DO
! 176: ELSE
! 177: K = 1
! 178: DO WHILE ( K .LE. NCOLS )
! 179: IF ( IPIV( K ).GT.0 ) THEN
! 180: ! 1x1 pivot
! 181: KP = IPIV( K )
! 182: IF ( KP .NE. K ) THEN
! 183: TMP = WORK( N+K )
! 184: WORK( N+K ) = WORK( N+KP )
! 185: WORK( N+KP ) = TMP
! 186: END IF
! 187: DO I = K, N
! 188: WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
! 189: END DO
! 190: K = K + 1
! 191: ELSE
! 192: ! 2x2 pivot
! 193: KP = -IPIV( K )
! 194: TMP = WORK( N+K+1 )
! 195: WORK( N+K+1 ) = WORK( N+KP )
! 196: WORK( N+KP ) = TMP
! 197: DO I = K+1, N
! 198: WORK( K ) = MAX( ABS( AF( I, K ) ), WORK( K ) )
! 199: WORK( K+1 ) = MAX( ABS( AF(I, K+1 ) ), WORK( K+1 ) )
! 200: END DO
! 201: WORK( K ) = MAX( ABS( AF( K, K ) ), WORK( K ) )
! 202: K = K + 2
! 203: END IF
! 204: END DO
! 205: K = NCOLS
! 206: DO WHILE ( K .GE. 1 )
! 207: IF ( IPIV( K ).GT.0 ) THEN
! 208: KP = IPIV( K )
! 209: IF ( KP .NE. K ) THEN
! 210: TMP = WORK( N+K )
! 211: WORK( N+K ) = WORK( N+KP )
! 212: WORK( N+KP ) = TMP
! 213: END IF
! 214: K = K - 1
! 215: ELSE
! 216: KP = -IPIV( K )
! 217: TMP = WORK( N+K )
! 218: WORK( N+K ) = WORK( N+KP )
! 219: WORK( N+KP ) = TMP
! 220: K = K - 2
! 221: ENDIF
! 222: END DO
! 223: END IF
! 224: *
! 225: * Compute the *inverse* of the max element growth factor. Dividing
! 226: * by zero would imply the largest entry of the factor's column is
! 227: * zero. Than can happen when either the column of A is zero or
! 228: * massive pivots made the factor underflow to zero. Neither counts
! 229: * as growth in itself, so simply ignore terms with zero
! 230: * denominators.
! 231: *
! 232: IF ( UPPER ) THEN
! 233: DO I = NCOLS, N
! 234: UMAX = WORK( I )
! 235: AMAX = WORK( N+I )
! 236: IF ( UMAX /= 0.0D+0 ) THEN
! 237: RPVGRW = MIN( AMAX / UMAX, RPVGRW )
! 238: END IF
! 239: END DO
! 240: ELSE
! 241: DO I = 1, NCOLS
! 242: UMAX = WORK( I )
! 243: AMAX = WORK( N+I )
! 244: IF ( UMAX /= 0.0D+0 ) THEN
! 245: RPVGRW = MIN( AMAX / UMAX, RPVGRW )
! 246: END IF
! 247: END DO
! 248: END IF
! 249:
! 250: DLA_SYRPVGRW = RPVGRW
! 251: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dla_syrpvgrw.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack