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