Annotation of rpl/lapack/lapack/zla_syrpvgrw.f, revision 1.1
1.1 ! bertrand 1: DOUBLE PRECISION FUNCTION ZLA_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: COMPLEX*16 A( LDA, * ), AF( LDAF, * )
! 20: DOUBLE PRECISION WORK( * )
! 21: INTEGER IPIV( * )
! 22: * ..
! 23: *
! 24: * Purpose
! 25: * =======
! 26: *
! 27: * ZLA_SYRPVGRW 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 ZSYTRF, .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 ZSYTRF.
! 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 ZSYTRF.
! 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
! 77: COMPLEX*16 ZDUM
! 78: * ..
! 79: * .. Intrinsic Functions ..
! 80: INTRINSIC ABS, REAL, DIMAG, MAX, MIN
! 81: * ..
! 82: * .. External Subroutines ..
! 83: EXTERNAL LSAME, ZLASET
! 84: LOGICAL LSAME
! 85: * ..
! 86: * .. Statement Functions ..
! 87: DOUBLE PRECISION CABS1
! 88: * ..
! 89: * .. Statement Function Definitions ..
! 90: CABS1( ZDUM ) = ABS( DBLE ( ZDUM ) ) + ABS( DIMAG ( ZDUM ) )
! 91: * ..
! 92: * .. Executable Statements ..
! 93: *
! 94: UPPER = LSAME( 'Upper', UPLO )
! 95: IF ( INFO.EQ.0 ) THEN
! 96: IF ( UPPER ) THEN
! 97: NCOLS = 1
! 98: ELSE
! 99: NCOLS = N
! 100: END IF
! 101: ELSE
! 102: NCOLS = INFO
! 103: END IF
! 104:
! 105: RPVGRW = 1.0D+0
! 106: DO I = 1, 2*N
! 107: WORK( I ) = 0.0D+0
! 108: END DO
! 109: *
! 110: * Find the max magnitude entry of each column of A. Compute the max
! 111: * for all N columns so we can apply the pivot permutation while
! 112: * looping below. Assume a full factorization is the common case.
! 113: *
! 114: IF ( UPPER ) THEN
! 115: DO J = 1, N
! 116: DO I = 1, J
! 117: WORK( N+I ) = MAX( CABS1( A( I, J ) ), WORK( N+I ) )
! 118: WORK( N+J ) = MAX( CABS1( A( I, J ) ), WORK( N+J ) )
! 119: END DO
! 120: END DO
! 121: ELSE
! 122: DO J = 1, N
! 123: DO I = J, N
! 124: WORK( N+I ) = MAX( CABS1( A( I, J ) ), WORK( N+I ) )
! 125: WORK( N+J ) = MAX( CABS1( A( I, J ) ), WORK( N+J ) )
! 126: END DO
! 127: END DO
! 128: END IF
! 129: *
! 130: * Now find the max magnitude entry of each column of U or L. Also
! 131: * permute the magnitudes of A above so they're in the same order as
! 132: * the factor.
! 133: *
! 134: * The iteration orders and permutations were copied from zsytrs.
! 135: * Calls to SSWAP would be severe overkill.
! 136: *
! 137: IF ( UPPER ) THEN
! 138: K = N
! 139: DO WHILE ( K .LT. NCOLS .AND. K.GT.0 )
! 140: IF ( IPIV( K ).GT.0 ) THEN
! 141: ! 1x1 pivot
! 142: KP = IPIV( K )
! 143: IF ( KP .NE. K ) THEN
! 144: TMP = WORK( N+K )
! 145: WORK( N+K ) = WORK( N+KP )
! 146: WORK( N+KP ) = TMP
! 147: END IF
! 148: DO I = 1, K
! 149: WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
! 150: END DO
! 151: K = K - 1
! 152: ELSE
! 153: ! 2x2 pivot
! 154: KP = -IPIV( K )
! 155: TMP = WORK( N+K-1 )
! 156: WORK( N+K-1 ) = WORK( N+KP )
! 157: WORK( N+KP ) = TMP
! 158: DO I = 1, K-1
! 159: WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
! 160: WORK( K-1 ) =
! 161: $ MAX( CABS1( AF( I, K-1 ) ), WORK( K-1 ) )
! 162: END DO
! 163: WORK( K ) = MAX( CABS1( AF( K, K ) ), WORK( K ) )
! 164: K = K - 2
! 165: END IF
! 166: END DO
! 167: K = NCOLS
! 168: DO WHILE ( K .LE. N )
! 169: IF ( IPIV( K ).GT.0 ) THEN
! 170: KP = IPIV( K )
! 171: IF ( KP .NE. K ) THEN
! 172: TMP = WORK( N+K )
! 173: WORK( N+K ) = WORK( N+KP )
! 174: WORK( N+KP ) = TMP
! 175: END IF
! 176: K = K + 1
! 177: ELSE
! 178: KP = -IPIV( K )
! 179: TMP = WORK( N+K )
! 180: WORK( N+K ) = WORK( N+KP )
! 181: WORK( N+KP ) = TMP
! 182: K = K + 2
! 183: END IF
! 184: END DO
! 185: ELSE
! 186: K = 1
! 187: DO WHILE ( K .LE. NCOLS )
! 188: IF ( IPIV( K ).GT.0 ) THEN
! 189: ! 1x1 pivot
! 190: KP = IPIV( K )
! 191: IF ( KP .NE. K ) THEN
! 192: TMP = WORK( N+K )
! 193: WORK( N+K ) = WORK( N+KP )
! 194: WORK( N+KP ) = TMP
! 195: END IF
! 196: DO I = K, N
! 197: WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
! 198: END DO
! 199: K = K + 1
! 200: ELSE
! 201: ! 2x2 pivot
! 202: KP = -IPIV( K )
! 203: TMP = WORK( N+K+1 )
! 204: WORK( N+K+1 ) = WORK( N+KP )
! 205: WORK( N+KP ) = TMP
! 206: DO I = K+1, N
! 207: WORK( K ) = MAX( CABS1( AF( I, K ) ), WORK( K ) )
! 208: WORK( K+1 ) =
! 209: $ MAX( CABS1( AF( I, K+1 ) ), WORK( K+1 ) )
! 210: END DO
! 211: WORK( K ) = MAX( CABS1( AF( K, K ) ), WORK( K ) )
! 212: K = K + 2
! 213: END IF
! 214: END DO
! 215: K = NCOLS
! 216: DO WHILE ( K .GE. 1 )
! 217: IF ( IPIV( K ).GT.0 ) THEN
! 218: KP = IPIV( K )
! 219: IF ( KP .NE. K ) THEN
! 220: TMP = WORK( N+K )
! 221: WORK( N+K ) = WORK( N+KP )
! 222: WORK( N+KP ) = TMP
! 223: END IF
! 224: K = K - 1
! 225: ELSE
! 226: KP = -IPIV( K )
! 227: TMP = WORK( N+K )
! 228: WORK( N+K ) = WORK( N+KP )
! 229: WORK( N+KP ) = TMP
! 230: K = K - 2
! 231: ENDIF
! 232: END DO
! 233: END IF
! 234: *
! 235: * Compute the *inverse* of the max element growth factor. Dividing
! 236: * by zero would imply the largest entry of the factor's column is
! 237: * zero. Than can happen when either the column of A is zero or
! 238: * massive pivots made the factor underflow to zero. Neither counts
! 239: * as growth in itself, so simply ignore terms with zero
! 240: * denominators.
! 241: *
! 242: IF ( UPPER ) THEN
! 243: DO I = NCOLS, N
! 244: UMAX = WORK( I )
! 245: AMAX = WORK( N+I )
! 246: IF ( UMAX /= 0.0D+0 ) THEN
! 247: RPVGRW = MIN( AMAX / UMAX, RPVGRW )
! 248: END IF
! 249: END DO
! 250: ELSE
! 251: DO I = 1, NCOLS
! 252: UMAX = WORK( I )
! 253: AMAX = WORK( N+I )
! 254: IF ( UMAX /= 0.0D+0 ) THEN
! 255: RPVGRW = MIN( AMAX / UMAX, RPVGRW )
! 256: END IF
! 257: END DO
! 258: END IF
! 259:
! 260: ZLA_SYRPVGRW = RPVGRW
! 261: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zla_syrpvgrw.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack