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