 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 DOUBLE PRECISION FUNCTION ZLANHF( NORM, TRANSR, UPLO, N, A, WORK ) 00002 * 00003 * -- LAPACK routine (version 3.3.1) -- 00004 * 00005 * -- Contributed by Fred Gustavson of the IBM Watson Research Center -- 00006 * -- April 2009 -- 00007 * 00008 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00009 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00010 * 00011 * .. Scalar Arguments .. 00012 CHARACTER NORM, TRANSR, UPLO 00013 INTEGER N 00014 * .. 00015 * .. Array Arguments .. 00016 DOUBLE PRECISION WORK( 0: * ) 00017 COMPLEX*16 A( 0: * ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * ZLANHF returns the value of the one norm, or the Frobenius norm, or 00024 * the infinity norm, or the element of largest absolute value of a 00025 * complex Hermitian matrix A in RFP format. 00026 * 00027 * Description 00028 * =========== 00029 * 00030 * ZLANHF returns the value 00031 * 00032 * ZLANHF = ( max(abs(A(i,j))), NORM = 'M' or 'm' 00033 * ( 00034 * ( norm1(A), NORM = '1', 'O' or 'o' 00035 * ( 00036 * ( normI(A), NORM = 'I' or 'i' 00037 * ( 00038 * ( normF(A), NORM = 'F', 'f', 'E' or 'e' 00039 * 00040 * where norm1 denotes the one norm of a matrix (maximum column sum), 00041 * normI denotes the infinity norm of a matrix (maximum row sum) and 00042 * normF denotes the Frobenius norm of a matrix (square root of sum of 00043 * squares). Note that max(abs(A(i,j))) is not a matrix norm. 00044 * 00045 * Arguments 00046 * ========= 00047 * 00048 * NORM (input) CHARACTER 00049 * Specifies the value to be returned in ZLANHF as described 00050 * above. 00051 * 00052 * TRANSR (input) CHARACTER 00053 * Specifies whether the RFP format of A is normal or 00054 * conjugate-transposed format. 00055 * = 'N': RFP format is Normal 00056 * = 'C': RFP format is Conjugate-transposed 00057 * 00058 * UPLO (input) CHARACTER 00059 * On entry, UPLO specifies whether the RFP matrix A came from 00060 * an upper or lower triangular matrix as follows: 00061 * 00062 * UPLO = 'U' or 'u' RFP A came from an upper triangular 00063 * matrix 00064 * 00065 * UPLO = 'L' or 'l' RFP A came from a lower triangular 00066 * matrix 00067 * 00068 * N (input) INTEGER 00069 * The order of the matrix A. N >= 0. When N = 0, ZLANHF is 00070 * set to zero. 00071 * 00072 * A (input) COMPLEX*16 array, dimension ( N*(N+1)/2 ); 00073 * On entry, the matrix A in RFP Format. 00074 * RFP Format is described by TRANSR, UPLO and N as follows: 00075 * If TRANSR='N' then RFP A is (0:N,0:K-1) when N is even; 00076 * K=N/2. RFP A is (0:N-1,0:K) when N is odd; K=N/2. If 00077 * TRANSR = 'C' then RFP is the Conjugate-transpose of RFP A 00078 * as defined when TRANSR = 'N'. The contents of RFP A are 00079 * defined by UPLO as follows: If UPLO = 'U' the RFP A 00080 * contains the ( N*(N+1)/2 ) elements of upper packed A 00081 * either in normal or conjugate-transpose Format. If 00082 * UPLO = 'L' the RFP A contains the ( N*(N+1) /2 ) elements 00083 * of lower packed A either in normal or conjugate-transpose 00084 * Format. The LDA of RFP A is (N+1)/2 when TRANSR = 'C'. When 00085 * TRANSR is 'N' the LDA is N+1 when N is even and is N when 00086 * is odd. See the Note below for more details. 00087 * Unchanged on exit. 00088 * 00089 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK), 00090 * where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, 00091 * WORK is not referenced. 00092 * 00093 * Further Details 00094 * =============== 00095 * 00096 * We first consider Standard Packed Format when N is even. 00097 * We give an example where N = 6. 00098 * 00099 * AP is Upper AP is Lower 00100 * 00101 * 00 01 02 03 04 05 00 00102 * 11 12 13 14 15 10 11 00103 * 22 23 24 25 20 21 22 00104 * 33 34 35 30 31 32 33 00105 * 44 45 40 41 42 43 44 00106 * 55 50 51 52 53 54 55 00107 * 00108 * 00109 * Let TRANSR = 'N'. RFP holds AP as follows: 00110 * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last 00111 * three columns of AP upper. The lower triangle A(4:6,0:2) consists of 00112 * conjugate-transpose of the first three columns of AP upper. 00113 * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first 00114 * three columns of AP lower. The upper triangle A(0:2,0:2) consists of 00115 * conjugate-transpose of the last three columns of AP lower. 00116 * To denote conjugate we place -- above the element. This covers the 00117 * case N even and TRANSR = 'N'. 00118 * 00119 * RFP A RFP A 00120 * 00121 * -- -- -- 00122 * 03 04 05 33 43 53 00123 * -- -- 00124 * 13 14 15 00 44 54 00125 * -- 00126 * 23 24 25 10 11 55 00127 * 00128 * 33 34 35 20 21 22 00129 * -- 00130 * 00 44 45 30 31 32 00131 * -- -- 00132 * 01 11 55 40 41 42 00133 * -- -- -- 00134 * 02 12 22 50 51 52 00135 * 00136 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- 00137 * transpose of RFP A above. One therefore gets: 00138 * 00139 * 00140 * RFP A RFP A 00141 * 00142 * -- -- -- -- -- -- -- -- -- -- 00143 * 03 13 23 33 00 01 02 33 00 10 20 30 40 50 00144 * -- -- -- -- -- -- -- -- -- -- 00145 * 04 14 24 34 44 11 12 43 44 11 21 31 41 51 00146 * -- -- -- -- -- -- -- -- -- -- 00147 * 05 15 25 35 45 55 22 53 54 55 22 32 42 52 00148 * 00149 * 00150 * We next consider Standard Packed Format when N is odd. 00151 * We give an example where N = 5. 00152 * 00153 * AP is Upper AP is Lower 00154 * 00155 * 00 01 02 03 04 00 00156 * 11 12 13 14 10 11 00157 * 22 23 24 20 21 22 00158 * 33 34 30 31 32 33 00159 * 44 40 41 42 43 44 00160 * 00161 * 00162 * Let TRANSR = 'N'. RFP holds AP as follows: 00163 * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last 00164 * three columns of AP upper. The lower triangle A(3:4,0:1) consists of 00165 * conjugate-transpose of the first two columns of AP upper. 00166 * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first 00167 * three columns of AP lower. The upper triangle A(0:1,1:2) consists of 00168 * conjugate-transpose of the last two columns of AP lower. 00169 * To denote conjugate we place -- above the element. This covers the 00170 * case N odd and TRANSR = 'N'. 00171 * 00172 * RFP A RFP A 00173 * 00174 * -- -- 00175 * 02 03 04 00 33 43 00176 * -- 00177 * 12 13 14 10 11 44 00178 * 00179 * 22 23 24 20 21 22 00180 * -- 00181 * 00 33 34 30 31 32 00182 * -- -- 00183 * 01 11 44 40 41 42 00184 * 00185 * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- 00186 * transpose of RFP A above. One therefore gets: 00187 * 00188 * 00189 * RFP A RFP A 00190 * 00191 * -- -- -- -- -- -- -- -- -- 00192 * 02 12 22 00 01 00 10 20 30 40 50 00193 * -- -- -- -- -- -- -- -- -- 00194 * 03 13 23 33 11 33 11 21 31 41 51 00195 * -- -- -- -- -- -- -- -- -- 00196 * 04 14 24 34 44 43 44 22 32 42 52 00197 * 00198 * ===================================================================== 00199 * 00200 * .. Parameters .. 00201 DOUBLE PRECISION ONE, ZERO 00202 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00203 * .. 00204 * .. Local Scalars .. 00205 INTEGER I, J, IFM, ILU, NOE, N1, K, L, LDA 00206 DOUBLE PRECISION SCALE, S, VALUE, AA 00207 * .. 00208 * .. External Functions .. 00209 LOGICAL LSAME 00210 INTEGER IDAMAX 00211 EXTERNAL LSAME, IDAMAX 00212 * .. 00213 * .. External Subroutines .. 00214 EXTERNAL ZLASSQ 00215 * .. 00216 * .. Intrinsic Functions .. 00217 INTRINSIC ABS, DBLE, MAX, SQRT 00218 * .. 00219 * .. Executable Statements .. 00220 * 00221 IF( N.EQ.0 ) THEN 00222 ZLANHF = ZERO 00223 RETURN 00224 END IF 00225 * 00226 * set noe = 1 if n is odd. if n is even set noe=0 00227 * 00228 NOE = 1 00229 IF( MOD( N, 2 ).EQ.0 ) 00230 $ NOE = 0 00231 * 00232 * set ifm = 0 when form='C' or 'c' and 1 otherwise 00233 * 00234 IFM = 1 00235 IF( LSAME( TRANSR, 'C' ) ) 00236 $ IFM = 0 00237 * 00238 * set ilu = 0 when uplo='U or 'u' and 1 otherwise 00239 * 00240 ILU = 1 00241 IF( LSAME( UPLO, 'U' ) ) 00242 $ ILU = 0 00243 * 00244 * set lda = (n+1)/2 when ifm = 0 00245 * set lda = n when ifm = 1 and noe = 1 00246 * set lda = n+1 when ifm = 1 and noe = 0 00247 * 00248 IF( IFM.EQ.1 ) THEN 00249 IF( NOE.EQ.1 ) THEN 00250 LDA = N 00251 ELSE 00252 * noe=0 00253 LDA = N + 1 00254 END IF 00255 ELSE 00256 * ifm=0 00257 LDA = ( N+1 ) / 2 00258 END IF 00259 * 00260 IF( LSAME( NORM, 'M' ) ) THEN 00261 * 00262 * Find max(abs(A(i,j))). 00263 * 00264 K = ( N+1 ) / 2 00265 VALUE = ZERO 00266 IF( NOE.EQ.1 ) THEN 00267 * n is odd & n = k + k - 1 00268 IF( IFM.EQ.1 ) THEN 00269 * A is n by k 00270 IF( ILU.EQ.1 ) THEN 00271 * uplo ='L' 00272 J = 0 00273 * -> L(0,0) 00274 VALUE = MAX( VALUE, ABS( DBLE( A( J+J*LDA ) ) ) ) 00275 DO I = 1, N - 1 00276 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00277 END DO 00278 DO J = 1, K - 1 00279 DO I = 0, J - 2 00280 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00281 END DO 00282 I = J - 1 00283 * L(k+j,k+j) 00284 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00285 I = J 00286 * -> L(j,j) 00287 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00288 DO I = J + 1, N - 1 00289 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00290 END DO 00291 END DO 00292 ELSE 00293 * uplo = 'U' 00294 DO J = 0, K - 2 00295 DO I = 0, K + J - 2 00296 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00297 END DO 00298 I = K + J - 1 00299 * -> U(i,i) 00300 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00301 I = I + 1 00302 * =k+j; i -> U(j,j) 00303 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00304 DO I = K + J + 1, N - 1 00305 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00306 END DO 00307 END DO 00308 DO I = 0, N - 2 00309 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00310 * j=k-1 00311 END DO 00312 * i=n-1 -> U(n-1,n-1) 00313 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00314 END IF 00315 ELSE 00316 * xpose case; A is k by n 00317 IF( ILU.EQ.1 ) THEN 00318 * uplo ='L' 00319 DO J = 0, K - 2 00320 DO I = 0, J - 1 00321 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00322 END DO 00323 I = J 00324 * L(i,i) 00325 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00326 I = J + 1 00327 * L(j+k,j+k) 00328 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00329 DO I = J + 2, K - 1 00330 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00331 END DO 00332 END DO 00333 J = K - 1 00334 DO I = 0, K - 2 00335 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00336 END DO 00337 I = K - 1 00338 * -> L(i,i) is at A(i,j) 00339 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00340 DO J = K, N - 1 00341 DO I = 0, K - 1 00342 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00343 END DO 00344 END DO 00345 ELSE 00346 * uplo = 'U' 00347 DO J = 0, K - 2 00348 DO I = 0, K - 1 00349 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00350 END DO 00351 END DO 00352 J = K - 1 00353 * -> U(j,j) is at A(0,j) 00354 VALUE = MAX( VALUE, ABS( DBLE( A( 0+J*LDA ) ) ) ) 00355 DO I = 1, K - 1 00356 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00357 END DO 00358 DO J = K, N - 1 00359 DO I = 0, J - K - 1 00360 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00361 END DO 00362 I = J - K 00363 * -> U(i,i) at A(i,j) 00364 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00365 I = J - K + 1 00366 * U(j,j) 00367 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00368 DO I = J - K + 2, K - 1 00369 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00370 END DO 00371 END DO 00372 END IF 00373 END IF 00374 ELSE 00375 * n is even & k = n/2 00376 IF( IFM.EQ.1 ) THEN 00377 * A is n+1 by k 00378 IF( ILU.EQ.1 ) THEN 00379 * uplo ='L' 00380 J = 0 00381 * -> L(k,k) & j=1 -> L(0,0) 00382 VALUE = MAX( VALUE, ABS( DBLE( A( J+J*LDA ) ) ) ) 00383 VALUE = MAX( VALUE, ABS( DBLE( A( J+1+J*LDA ) ) ) ) 00384 DO I = 2, N 00385 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00386 END DO 00387 DO J = 1, K - 1 00388 DO I = 0, J - 1 00389 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00390 END DO 00391 I = J 00392 * L(k+j,k+j) 00393 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00394 I = J + 1 00395 * -> L(j,j) 00396 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00397 DO I = J + 2, N 00398 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00399 END DO 00400 END DO 00401 ELSE 00402 * uplo = 'U' 00403 DO J = 0, K - 2 00404 DO I = 0, K + J - 1 00405 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00406 END DO 00407 I = K + J 00408 * -> U(i,i) 00409 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00410 I = I + 1 00411 * =k+j+1; i -> U(j,j) 00412 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00413 DO I = K + J + 2, N 00414 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00415 END DO 00416 END DO 00417 DO I = 0, N - 2 00418 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00419 * j=k-1 00420 END DO 00421 * i=n-1 -> U(n-1,n-1) 00422 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00423 I = N 00424 * -> U(k-1,k-1) 00425 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00426 END IF 00427 ELSE 00428 * xpose case; A is k by n+1 00429 IF( ILU.EQ.1 ) THEN 00430 * uplo ='L' 00431 J = 0 00432 * -> L(k,k) at A(0,0) 00433 VALUE = MAX( VALUE, ABS( DBLE( A( J+J*LDA ) ) ) ) 00434 DO I = 1, K - 1 00435 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00436 END DO 00437 DO J = 1, K - 1 00438 DO I = 0, J - 2 00439 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00440 END DO 00441 I = J - 1 00442 * L(i,i) 00443 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00444 I = J 00445 * L(j+k,j+k) 00446 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00447 DO I = J + 1, K - 1 00448 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00449 END DO 00450 END DO 00451 J = K 00452 DO I = 0, K - 2 00453 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00454 END DO 00455 I = K - 1 00456 * -> L(i,i) is at A(i,j) 00457 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00458 DO J = K + 1, N 00459 DO I = 0, K - 1 00460 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00461 END DO 00462 END DO 00463 ELSE 00464 * uplo = 'U' 00465 DO J = 0, K - 1 00466 DO I = 0, K - 1 00467 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00468 END DO 00469 END DO 00470 J = K 00471 * -> U(j,j) is at A(0,j) 00472 VALUE = MAX( VALUE, ABS( DBLE( A( 0+J*LDA ) ) ) ) 00473 DO I = 1, K - 1 00474 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00475 END DO 00476 DO J = K + 1, N - 1 00477 DO I = 0, J - K - 2 00478 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00479 END DO 00480 I = J - K - 1 00481 * -> U(i,i) at A(i,j) 00482 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00483 I = J - K 00484 * U(j,j) 00485 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00486 DO I = J - K + 1, K - 1 00487 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00488 END DO 00489 END DO 00490 J = N 00491 DO I = 0, K - 2 00492 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00493 END DO 00494 I = K - 1 00495 * U(k,k) at A(i,j) 00496 VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 00497 END IF 00498 END IF 00499 END IF 00500 ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR. 00501 $ ( NORM.EQ.'1' ) ) THEN 00502 * 00503 * Find normI(A) ( = norm1(A), since A is Hermitian). 00504 * 00505 IF( IFM.EQ.1 ) THEN 00506 * A is 'N' 00507 K = N / 2 00508 IF( NOE.EQ.1 ) THEN 00509 * n is odd & A is n by (n+1)/2 00510 IF( ILU.EQ.0 ) THEN 00511 * uplo = 'U' 00512 DO I = 0, K - 1 00513 WORK( I ) = ZERO 00514 END DO 00515 DO J = 0, K 00516 S = ZERO 00517 DO I = 0, K + J - 1 00518 AA = ABS( A( I+J*LDA ) ) 00519 * -> A(i,j+k) 00520 S = S + AA 00521 WORK( I ) = WORK( I ) + AA 00522 END DO 00523 AA = ABS( DBLE( A( I+J*LDA ) ) ) 00524 * -> A(j+k,j+k) 00525 WORK( J+K ) = S + AA 00526 IF( I.EQ.K+K ) 00527 $ GO TO 10 00528 I = I + 1 00529 AA = ABS( DBLE( A( I+J*LDA ) ) ) 00530 * -> A(j,j) 00531 WORK( J ) = WORK( J ) + AA 00532 S = ZERO 00533 DO L = J + 1, K - 1 00534 I = I + 1 00535 AA = ABS( A( I+J*LDA ) ) 00536 * -> A(l,j) 00537 S = S + AA 00538 WORK( L ) = WORK( L ) + AA 00539 END DO 00540 WORK( J ) = WORK( J ) + S 00541 END DO 00542 10 CONTINUE 00543 I = IDAMAX( N, WORK, 1 ) 00544 VALUE = WORK( I-1 ) 00545 ELSE 00546 * ilu = 1 & uplo = 'L' 00547 K = K + 1 00548 * k=(n+1)/2 for n odd and ilu=1 00549 DO I = K, N - 1 00550 WORK( I ) = ZERO 00551 END DO 00552 DO J = K - 1, 0, -1 00553 S = ZERO 00554 DO I = 0, J - 2 00555 AA = ABS( A( I+J*LDA ) ) 00556 * -> A(j+k,i+k) 00557 S = S + AA 00558 WORK( I+K ) = WORK( I+K ) + AA 00559 END DO 00560 IF( J.GT.0 ) THEN 00561 AA = ABS( DBLE( A( I+J*LDA ) ) ) 00562 * -> A(j+k,j+k) 00563 S = S + AA 00564 WORK( I+K ) = WORK( I+K ) + S 00565 * i=j 00566 I = I + 1 00567 END IF 00568 AA = ABS( DBLE( A( I+J*LDA ) ) ) 00569 * -> A(j,j) 00570 WORK( J ) = AA 00571 S = ZERO 00572 DO L = J + 1, N - 1 00573 I = I + 1 00574 AA = ABS( A( I+J*LDA ) ) 00575 * -> A(l,j) 00576 S = S + AA 00577 WORK( L ) = WORK( L ) + AA 00578 END DO 00579 WORK( J ) = WORK( J ) + S 00580 END DO 00581 I = IDAMAX( N, WORK, 1 ) 00582 VALUE = WORK( I-1 ) 00583 END IF 00584 ELSE 00585 * n is even & A is n+1 by k = n/2 00586 IF( ILU.EQ.0 ) THEN 00587 * uplo = 'U' 00588 DO I = 0, K - 1 00589 WORK( I ) = ZERO 00590 END DO 00591 DO J = 0, K - 1 00592 S = ZERO 00593 DO I = 0, K + J - 1 00594 AA = ABS( A( I+J*LDA ) ) 00595 * -> A(i,j+k) 00596 S = S + AA 00597 WORK( I ) = WORK( I ) + AA 00598 END DO 00599 AA = ABS( DBLE( A( I+J*LDA ) ) ) 00600 * -> A(j+k,j+k) 00601 WORK( J+K ) = S + AA 00602 I = I + 1 00603 AA = ABS( DBLE( A( I+J*LDA ) ) ) 00604 * -> A(j,j) 00605 WORK( J ) = WORK( J ) + AA 00606 S = ZERO 00607 DO L = J + 1, K - 1 00608 I = I + 1 00609 AA = ABS( A( I+J*LDA ) ) 00610 * -> A(l,j) 00611 S = S + AA 00612 WORK( L ) = WORK( L ) + AA 00613 END DO 00614 WORK( J ) = WORK( J ) + S 00615 END DO 00616 I = IDAMAX( N, WORK, 1 ) 00617 VALUE = WORK( I-1 ) 00618 ELSE 00619 * ilu = 1 & uplo = 'L' 00620 DO I = K, N - 1 00621 WORK( I ) = ZERO 00622 END DO 00623 DO J = K - 1, 0, -1 00624 S = ZERO 00625 DO I = 0, J - 1 00626 AA = ABS( A( I+J*LDA ) ) 00627 * -> A(j+k,i+k) 00628 S = S + AA 00629 WORK( I+K ) = WORK( I+K ) + AA 00630 END DO 00631 AA = ABS( DBLE( A( I+J*LDA ) ) ) 00632 * -> A(j+k,j+k) 00633 S = S + AA 00634 WORK( I+K ) = WORK( I+K ) + S 00635 * i=j 00636 I = I + 1 00637 AA = ABS( DBLE( A( I+J*LDA ) ) ) 00638 * -> A(j,j) 00639 WORK( J ) = AA 00640 S = ZERO 00641 DO L = J + 1, N - 1 00642 I = I + 1 00643 AA = ABS( A( I+J*LDA ) ) 00644 * -> A(l,j) 00645 S = S + AA 00646 WORK( L ) = WORK( L ) + AA 00647 END DO 00648 WORK( J ) = WORK( J ) + S 00649 END DO 00650 I = IDAMAX( N, WORK, 1 ) 00651 VALUE = WORK( I-1 ) 00652 END IF 00653 END IF 00654 ELSE 00655 * ifm=0 00656 K = N / 2 00657 IF( NOE.EQ.1 ) THEN 00658 * n is odd & A is (n+1)/2 by n 00659 IF( ILU.EQ.0 ) THEN 00660 * uplo = 'U' 00661 N1 = K 00662 * n/2 00663 K = K + 1 00664 * k is the row size and lda 00665 DO I = N1, N - 1 00666 WORK( I ) = ZERO 00667 END DO 00668 DO J = 0, N1 - 1 00669 S = ZERO 00670 DO I = 0, K - 1 00671 AA = ABS( A( I+J*LDA ) ) 00672 * A(j,n1+i) 00673 WORK( I+N1 ) = WORK( I+N1 ) + AA 00674 S = S + AA 00675 END DO 00676 WORK( J ) = S 00677 END DO 00678 * j=n1=k-1 is special 00679 S = ABS( DBLE( A( 0+J*LDA ) ) ) 00680 * A(k-1,k-1) 00681 DO I = 1, K - 1 00682 AA = ABS( A( I+J*LDA ) ) 00683 * A(k-1,i+n1) 00684 WORK( I+N1 ) = WORK( I+N1 ) + AA 00685 S = S + AA 00686 END DO 00687 WORK( J ) = WORK( J ) + S 00688 DO J = K, N - 1 00689 S = ZERO 00690 DO I = 0, J - K - 1 00691 AA = ABS( A( I+J*LDA ) ) 00692 * A(i,j-k) 00693 WORK( I ) = WORK( I ) + AA 00694 S = S + AA 00695 END DO 00696 * i=j-k 00697 AA = ABS( DBLE( A( I+J*LDA ) ) ) 00698 * A(j-k,j-k) 00699 S = S + AA 00700 WORK( J-K ) = WORK( J-K ) + S 00701 I = I + 1 00702 S = ABS( DBLE( A( I+J*LDA ) ) ) 00703 * A(j,j) 00704 DO L = J + 1, N - 1 00705 I = I + 1 00706 AA = ABS( A( I+J*LDA ) ) 00707 * A(j,l) 00708 WORK( L ) = WORK( L ) + AA 00709 S = S + AA 00710 END DO 00711 WORK( J ) = WORK( J ) + S 00712 END DO 00713 I = IDAMAX( N, WORK, 1 ) 00714 VALUE = WORK( I-1 ) 00715 ELSE 00716 * ilu=1 & uplo = 'L' 00717 K = K + 1 00718 * k=(n+1)/2 for n odd and ilu=1 00719 DO I = K, N - 1 00720 WORK( I ) = ZERO 00721 END DO 00722 DO J = 0, K - 2 00723 * process 00724 S = ZERO 00725 DO I = 0, J - 1 00726 AA = ABS( A( I+J*LDA ) ) 00727 * A(j,i) 00728 WORK( I ) = WORK( I ) + AA 00729 S = S + AA 00730 END DO 00731 AA = ABS( DBLE( A( I+J*LDA ) ) ) 00732 * i=j so process of A(j,j) 00733 S = S + AA 00734 WORK( J ) = S 00735 * is initialised here 00736 I = I + 1 00737 * i=j process A(j+k,j+k) 00738 AA = ABS( DBLE( A( I+J*LDA ) ) ) 00739 S = AA 00740 DO L = K + J + 1, N - 1 00741 I = I + 1 00742 AA = ABS( A( I+J*LDA ) ) 00743 * A(l,k+j) 00744 S = S + AA 00745 WORK( L ) = WORK( L ) + AA 00746 END DO 00747 WORK( K+J ) = WORK( K+J ) + S 00748 END DO 00749 * j=k-1 is special :process col A(k-1,0:k-1) 00750 S = ZERO 00751 DO I = 0, K - 2 00752 AA = ABS( A( I+J*LDA ) ) 00753 * A(k,i) 00754 WORK( I ) = WORK( I ) + AA 00755 S = S + AA 00756 END DO 00757 * i=k-1 00758 AA = ABS( DBLE( A( I+J*LDA ) ) ) 00759 * A(k-1,k-1) 00760 S = S + AA 00761 WORK( I ) = S 00762 * done with col j=k+1 00763 DO J = K, N - 1 00764 * process col j of A = A(j,0:k-1) 00765 S = ZERO 00766 DO I = 0, K - 1 00767 AA = ABS( A( I+J*LDA ) ) 00768 * A(j,i) 00769 WORK( I ) = WORK( I ) + AA 00770 S = S + AA 00771 END DO 00772 WORK( J ) = WORK( J ) + S 00773 END DO 00774 I = IDAMAX( N, WORK, 1 ) 00775 VALUE = WORK( I-1 ) 00776 END IF 00777 ELSE 00778 * n is even & A is k=n/2 by n+1 00779 IF( ILU.EQ.0 ) THEN 00780 * uplo = 'U' 00781 DO I = K, N - 1 00782 WORK( I ) = ZERO 00783 END DO 00784 DO J = 0, K - 1 00785 S = ZERO 00786 DO I = 0, K - 1 00787 AA = ABS( A( I+J*LDA ) ) 00788 * A(j,i+k) 00789 WORK( I+K ) = WORK( I+K ) + AA 00790 S = S + AA 00791 END DO 00792 WORK( J ) = S 00793 END DO 00794 * j=k 00795 AA = ABS( DBLE( A( 0+J*LDA ) ) ) 00796 * A(k,k) 00797 S = AA 00798 DO I = 1, K - 1 00799 AA = ABS( A( I+J*LDA ) ) 00800 * A(k,k+i) 00801 WORK( I+K ) = WORK( I+K ) + AA 00802 S = S + AA 00803 END DO 00804 WORK( J ) = WORK( J ) + S 00805 DO J = K + 1, N - 1 00806 S = ZERO 00807 DO I = 0, J - 2 - K 00808 AA = ABS( A( I+J*LDA ) ) 00809 * A(i,j-k-1) 00810 WORK( I ) = WORK( I ) + AA 00811 S = S + AA 00812 END DO 00813 * i=j-1-k 00814 AA = ABS( DBLE( A( I+J*LDA ) ) ) 00815 * A(j-k-1,j-k-1) 00816 S = S + AA 00817 WORK( J-K-1 ) = WORK( J-K-1 ) + S 00818 I = I + 1 00819 AA = ABS( DBLE( A( I+J*LDA ) ) ) 00820 * A(j,j) 00821 S = AA 00822 DO L = J + 1, N - 1 00823 I = I + 1 00824 AA = ABS( A( I+J*LDA ) ) 00825 * A(j,l) 00826 WORK( L ) = WORK( L ) + AA 00827 S = S + AA 00828 END DO 00829 WORK( J ) = WORK( J ) + S 00830 END DO 00831 * j=n 00832 S = ZERO 00833 DO I = 0, K - 2 00834 AA = ABS( A( I+J*LDA ) ) 00835 * A(i,k-1) 00836 WORK( I ) = WORK( I ) + AA 00837 S = S + AA 00838 END DO 00839 * i=k-1 00840 AA = ABS( DBLE( A( I+J*LDA ) ) ) 00841 * A(k-1,k-1) 00842 S = S + AA 00843 WORK( I ) = WORK( I ) + S 00844 I = IDAMAX( N, WORK, 1 ) 00845 VALUE = WORK( I-1 ) 00846 ELSE 00847 * ilu=1 & uplo = 'L' 00848 DO I = K, N - 1 00849 WORK( I ) = ZERO 00850 END DO 00851 * j=0 is special :process col A(k:n-1,k) 00852 S = ABS( DBLE( A( 0 ) ) ) 00853 * A(k,k) 00854 DO I = 1, K - 1 00855 AA = ABS( A( I ) ) 00856 * A(k+i,k) 00857 WORK( I+K ) = WORK( I+K ) + AA 00858 S = S + AA 00859 END DO 00860 WORK( K ) = WORK( K ) + S 00861 DO J = 1, K - 1 00862 * process 00863 S = ZERO 00864 DO I = 0, J - 2 00865 AA = ABS( A( I+J*LDA ) ) 00866 * A(j-1,i) 00867 WORK( I ) = WORK( I ) + AA 00868 S = S + AA 00869 END DO 00870 AA = ABS( DBLE( A( I+J*LDA ) ) ) 00871 * i=j-1 so process of A(j-1,j-1) 00872 S = S + AA 00873 WORK( J-1 ) = S 00874 * is initialised here 00875 I = I + 1 00876 * i=j process A(j+k,j+k) 00877 AA = ABS( DBLE( A( I+J*LDA ) ) ) 00878 S = AA 00879 DO L = K + J + 1, N - 1 00880 I = I + 1 00881 AA = ABS( A( I+J*LDA ) ) 00882 * A(l,k+j) 00883 S = S + AA 00884 WORK( L ) = WORK( L ) + AA 00885 END DO 00886 WORK( K+J ) = WORK( K+J ) + S 00887 END DO 00888 * j=k is special :process col A(k,0:k-1) 00889 S = ZERO 00890 DO I = 0, K - 2 00891 AA = ABS( A( I+J*LDA ) ) 00892 * A(k,i) 00893 WORK( I ) = WORK( I ) + AA 00894 S = S + AA 00895 END DO 00896 * 00897 * i=k-1 00898 AA = ABS( DBLE( A( I+J*LDA ) ) ) 00899 * A(k-1,k-1) 00900 S = S + AA 00901 WORK( I ) = S 00902 * done with col j=k+1 00903 DO J = K + 1, N 00904 * 00905 * process col j-1 of A = A(j-1,0:k-1) 00906 S = ZERO 00907 DO I = 0, K - 1 00908 AA = ABS( A( I+J*LDA ) ) 00909 * A(j-1,i) 00910 WORK( I ) = WORK( I ) + AA 00911 S = S + AA 00912 END DO 00913 WORK( J-1 ) = WORK( J-1 ) + S 00914 END DO 00915 I = IDAMAX( N, WORK, 1 ) 00916 VALUE = WORK( I-1 ) 00917 END IF 00918 END IF 00919 END IF 00920 ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN 00921 * 00922 * Find normF(A). 00923 * 00924 K = ( N+1 ) / 2 00925 SCALE = ZERO 00926 S = ONE 00927 IF( NOE.EQ.1 ) THEN 00928 * n is odd 00929 IF( IFM.EQ.1 ) THEN 00930 * A is normal & A is n by k 00931 IF( ILU.EQ.0 ) THEN 00932 * A is upper 00933 DO J = 0, K - 3 00934 CALL ZLASSQ( K-J-2, A( K+J+1+J*LDA ), 1, SCALE, S ) 00935 * L at A(k,0) 00936 END DO 00937 DO J = 0, K - 1 00938 CALL ZLASSQ( K+J-1, A( 0+J*LDA ), 1, SCALE, S ) 00939 * trap U at A(0,0) 00940 END DO 00941 S = S + S 00942 * double s for the off diagonal elements 00943 L = K - 1 00944 * -> U(k,k) at A(k-1,0) 00945 DO I = 0, K - 2 00946 AA = DBLE( A( L ) ) 00947 * U(k+i,k+i) 00948 IF( AA.NE.ZERO ) THEN 00949 IF( SCALE.LT.AA ) THEN 00950 S = ONE + S*( SCALE / AA )**2 00951 SCALE = AA 00952 ELSE 00953 S = S + ( AA / SCALE )**2 00954 END IF 00955 END IF 00956 AA = DBLE( A( L+1 ) ) 00957 * U(i,i) 00958 IF( AA.NE.ZERO ) THEN 00959 IF( SCALE.LT.AA ) THEN 00960 S = ONE + S*( SCALE / AA )**2 00961 SCALE = AA 00962 ELSE 00963 S = S + ( AA / SCALE )**2 00964 END IF 00965 END IF 00966 L = L + LDA + 1 00967 END DO 00968 AA = DBLE( A( L ) ) 00969 * U(n-1,n-1) 00970 IF( AA.NE.ZERO ) THEN 00971 IF( SCALE.LT.AA ) THEN 00972 S = ONE + S*( SCALE / AA )**2 00973 SCALE = AA 00974 ELSE 00975 S = S + ( AA / SCALE )**2 00976 END IF 00977 END IF 00978 ELSE 00979 * ilu=1 & A is lower 00980 DO J = 0, K - 1 00981 CALL ZLASSQ( N-J-1, A( J+1+J*LDA ), 1, SCALE, S ) 00982 * trap L at A(0,0) 00983 END DO 00984 DO J = 1, K - 2 00985 CALL ZLASSQ( J, A( 0+( 1+J )*LDA ), 1, SCALE, S ) 00986 * U at A(0,1) 00987 END DO 00988 S = S + S 00989 * double s for the off diagonal elements 00990 AA = DBLE( A( 0 ) ) 00991 * L(0,0) at A(0,0) 00992 IF( AA.NE.ZERO ) THEN 00993 IF( SCALE.LT.AA ) THEN 00994 S = ONE + S*( SCALE / AA )**2 00995 SCALE = AA 00996 ELSE 00997 S = S + ( AA / SCALE )**2 00998 END IF 00999 END IF 01000 L = LDA 01001 * -> L(k,k) at A(0,1) 01002 DO I = 1, K - 1 01003 AA = DBLE( A( L ) ) 01004 * L(k-1+i,k-1+i) 01005 IF( AA.NE.ZERO ) THEN 01006 IF( SCALE.LT.AA ) THEN 01007 S = ONE + S*( SCALE / AA )**2 01008 SCALE = AA 01009 ELSE 01010 S = S + ( AA / SCALE )**2 01011 END IF 01012 END IF 01013 AA = DBLE( A( L+1 ) ) 01014 * L(i,i) 01015 IF( AA.NE.ZERO ) THEN 01016 IF( SCALE.LT.AA ) THEN 01017 S = ONE + S*( SCALE / AA )**2 01018 SCALE = AA 01019 ELSE 01020 S = S + ( AA / SCALE )**2 01021 END IF 01022 END IF 01023 L = L + LDA + 1 01024 END DO 01025 END IF 01026 ELSE 01027 * A is xpose & A is k by n 01028 IF( ILU.EQ.0 ) THEN 01029 * A**H is upper 01030 DO J = 1, K - 2 01031 CALL ZLASSQ( J, A( 0+( K+J )*LDA ), 1, SCALE, S ) 01032 * U at A(0,k) 01033 END DO 01034 DO J = 0, K - 2 01035 CALL ZLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 01036 * k by k-1 rect. at A(0,0) 01037 END DO 01038 DO J = 0, K - 2 01039 CALL ZLASSQ( K-J-1, A( J+1+( J+K-1 )*LDA ), 1, 01040 $ SCALE, S ) 01041 * L at A(0,k-1) 01042 END DO 01043 S = S + S 01044 * double s for the off diagonal elements 01045 L = 0 + K*LDA - LDA 01046 * -> U(k-1,k-1) at A(0,k-1) 01047 AA = DBLE( A( L ) ) 01048 * U(k-1,k-1) 01049 IF( AA.NE.ZERO ) THEN 01050 IF( SCALE.LT.AA ) THEN 01051 S = ONE + S*( SCALE / AA )**2 01052 SCALE = AA 01053 ELSE 01054 S = S + ( AA / SCALE )**2 01055 END IF 01056 END IF 01057 L = L + LDA 01058 * -> U(0,0) at A(0,k) 01059 DO J = K, N - 1 01060 AA = DBLE( A( L ) ) 01061 * -> U(j-k,j-k) 01062 IF( AA.NE.ZERO ) THEN 01063 IF( SCALE.LT.AA ) THEN 01064 S = ONE + S*( SCALE / AA )**2 01065 SCALE = AA 01066 ELSE 01067 S = S + ( AA / SCALE )**2 01068 END IF 01069 END IF 01070 AA = DBLE( A( L+1 ) ) 01071 * -> U(j,j) 01072 IF( AA.NE.ZERO ) THEN 01073 IF( SCALE.LT.AA ) THEN 01074 S = ONE + S*( SCALE / AA )**2 01075 SCALE = AA 01076 ELSE 01077 S = S + ( AA / SCALE )**2 01078 END IF 01079 END IF 01080 L = L + LDA + 1 01081 END DO 01082 ELSE 01083 * A**H is lower 01084 DO J = 1, K - 1 01085 CALL ZLASSQ( J, A( 0+J*LDA ), 1, SCALE, S ) 01086 * U at A(0,0) 01087 END DO 01088 DO J = K, N - 1 01089 CALL ZLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 01090 * k by k-1 rect. at A(0,k) 01091 END DO 01092 DO J = 0, K - 3 01093 CALL ZLASSQ( K-J-2, A( J+2+J*LDA ), 1, SCALE, S ) 01094 * L at A(1,0) 01095 END DO 01096 S = S + S 01097 * double s for the off diagonal elements 01098 L = 0 01099 * -> L(0,0) at A(0,0) 01100 DO I = 0, K - 2 01101 AA = DBLE( A( L ) ) 01102 * L(i,i) 01103 IF( AA.NE.ZERO ) THEN 01104 IF( SCALE.LT.AA ) THEN 01105 S = ONE + S*( SCALE / AA )**2 01106 SCALE = AA 01107 ELSE 01108 S = S + ( AA / SCALE )**2 01109 END IF 01110 END IF 01111 AA = DBLE( A( L+1 ) ) 01112 * L(k+i,k+i) 01113 IF( AA.NE.ZERO ) THEN 01114 IF( SCALE.LT.AA ) THEN 01115 S = ONE + S*( SCALE / AA )**2 01116 SCALE = AA 01117 ELSE 01118 S = S + ( AA / SCALE )**2 01119 END IF 01120 END IF 01121 L = L + LDA + 1 01122 END DO 01123 * L-> k-1 + (k-1)*lda or L(k-1,k-1) at A(k-1,k-1) 01124 AA = DBLE( A( L ) ) 01125 * L(k-1,k-1) at A(k-1,k-1) 01126 IF( AA.NE.ZERO ) THEN 01127 IF( SCALE.LT.AA ) THEN 01128 S = ONE + S*( SCALE / AA )**2 01129 SCALE = AA 01130 ELSE 01131 S = S + ( AA / SCALE )**2 01132 END IF 01133 END IF 01134 END IF 01135 END IF 01136 ELSE 01137 * n is even 01138 IF( IFM.EQ.1 ) THEN 01139 * A is normal 01140 IF( ILU.EQ.0 ) THEN 01141 * A is upper 01142 DO J = 0, K - 2 01143 CALL ZLASSQ( K-J-1, A( K+J+2+J*LDA ), 1, SCALE, S ) 01144 * L at A(k+1,0) 01145 END DO 01146 DO J = 0, K - 1 01147 CALL ZLASSQ( K+J, A( 0+J*LDA ), 1, SCALE, S ) 01148 * trap U at A(0,0) 01149 END DO 01150 S = S + S 01151 * double s for the off diagonal elements 01152 L = K 01153 * -> U(k,k) at A(k,0) 01154 DO I = 0, K - 1 01155 AA = DBLE( A( L ) ) 01156 * U(k+i,k+i) 01157 IF( AA.NE.ZERO ) THEN 01158 IF( SCALE.LT.AA ) THEN 01159 S = ONE + S*( SCALE / AA )**2 01160 SCALE = AA 01161 ELSE 01162 S = S + ( AA / SCALE )**2 01163 END IF 01164 END IF 01165 AA = DBLE( A( L+1 ) ) 01166 * U(i,i) 01167 IF( AA.NE.ZERO ) THEN 01168 IF( SCALE.LT.AA ) THEN 01169 S = ONE + S*( SCALE / AA )**2 01170 SCALE = AA 01171 ELSE 01172 S = S + ( AA / SCALE )**2 01173 END IF 01174 END IF 01175 L = L + LDA + 1 01176 END DO 01177 ELSE 01178 * ilu=1 & A is lower 01179 DO J = 0, K - 1 01180 CALL ZLASSQ( N-J-1, A( J+2+J*LDA ), 1, SCALE, S ) 01181 * trap L at A(1,0) 01182 END DO 01183 DO J = 1, K - 1 01184 CALL ZLASSQ( J, A( 0+J*LDA ), 1, SCALE, S ) 01185 * U at A(0,0) 01186 END DO 01187 S = S + S 01188 * double s for the off diagonal elements 01189 L = 0 01190 * -> L(k,k) at A(0,0) 01191 DO I = 0, K - 1 01192 AA = DBLE( A( L ) ) 01193 * L(k-1+i,k-1+i) 01194 IF( AA.NE.ZERO ) THEN 01195 IF( SCALE.LT.AA ) THEN 01196 S = ONE + S*( SCALE / AA )**2 01197 SCALE = AA 01198 ELSE 01199 S = S + ( AA / SCALE )**2 01200 END IF 01201 END IF 01202 AA = DBLE( A( L+1 ) ) 01203 * L(i,i) 01204 IF( AA.NE.ZERO ) THEN 01205 IF( SCALE.LT.AA ) THEN 01206 S = ONE + S*( SCALE / AA )**2 01207 SCALE = AA 01208 ELSE 01209 S = S + ( AA / SCALE )**2 01210 END IF 01211 END IF 01212 L = L + LDA + 1 01213 END DO 01214 END IF 01215 ELSE 01216 * A is xpose 01217 IF( ILU.EQ.0 ) THEN 01218 * A**H is upper 01219 DO J = 1, K - 1 01220 CALL ZLASSQ( J, A( 0+( K+1+J )*LDA ), 1, SCALE, S ) 01221 * U at A(0,k+1) 01222 END DO 01223 DO J = 0, K - 1 01224 CALL ZLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 01225 * k by k rect. at A(0,0) 01226 END DO 01227 DO J = 0, K - 2 01228 CALL ZLASSQ( K-J-1, A( J+1+( J+K )*LDA ), 1, SCALE, 01229 $ S ) 01230 * L at A(0,k) 01231 END DO 01232 S = S + S 01233 * double s for the off diagonal elements 01234 L = 0 + K*LDA 01235 * -> U(k,k) at A(0,k) 01236 AA = DBLE( A( L ) ) 01237 * U(k,k) 01238 IF( AA.NE.ZERO ) THEN 01239 IF( SCALE.LT.AA ) THEN 01240 S = ONE + S*( SCALE / AA )**2 01241 SCALE = AA 01242 ELSE 01243 S = S + ( AA / SCALE )**2 01244 END IF 01245 END IF 01246 L = L + LDA 01247 * -> U(0,0) at A(0,k+1) 01248 DO J = K + 1, N - 1 01249 AA = DBLE( A( L ) ) 01250 * -> U(j-k-1,j-k-1) 01251 IF( AA.NE.ZERO ) THEN 01252 IF( SCALE.LT.AA ) THEN 01253 S = ONE + S*( SCALE / AA )**2 01254 SCALE = AA 01255 ELSE 01256 S = S + ( AA / SCALE )**2 01257 END IF 01258 END IF 01259 AA = DBLE( A( L+1 ) ) 01260 * -> U(j,j) 01261 IF( AA.NE.ZERO ) THEN 01262 IF( SCALE.LT.AA ) THEN 01263 S = ONE + S*( SCALE / AA )**2 01264 SCALE = AA 01265 ELSE 01266 S = S + ( AA / SCALE )**2 01267 END IF 01268 END IF 01269 L = L + LDA + 1 01270 END DO 01271 * L=k-1+n*lda 01272 * -> U(k-1,k-1) at A(k-1,n) 01273 AA = DBLE( A( L ) ) 01274 * U(k,k) 01275 IF( AA.NE.ZERO ) THEN 01276 IF( SCALE.LT.AA ) THEN 01277 S = ONE + S*( SCALE / AA )**2 01278 SCALE = AA 01279 ELSE 01280 S = S + ( AA / SCALE )**2 01281 END IF 01282 END IF 01283 ELSE 01284 * A**H is lower 01285 DO J = 1, K - 1 01286 CALL ZLASSQ( J, A( 0+( J+1 )*LDA ), 1, SCALE, S ) 01287 * U at A(0,1) 01288 END DO 01289 DO J = K + 1, N 01290 CALL ZLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 01291 * k by k rect. at A(0,k+1) 01292 END DO 01293 DO J = 0, K - 2 01294 CALL ZLASSQ( K-J-1, A( J+1+J*LDA ), 1, SCALE, S ) 01295 * L at A(0,0) 01296 END DO 01297 S = S + S 01298 * double s for the off diagonal elements 01299 L = 0 01300 * -> L(k,k) at A(0,0) 01301 AA = DBLE( A( L ) ) 01302 * L(k,k) at A(0,0) 01303 IF( AA.NE.ZERO ) THEN 01304 IF( SCALE.LT.AA ) THEN 01305 S = ONE + S*( SCALE / AA )**2 01306 SCALE = AA 01307 ELSE 01308 S = S + ( AA / SCALE )**2 01309 END IF 01310 END IF 01311 L = LDA 01312 * -> L(0,0) at A(0,1) 01313 DO I = 0, K - 2 01314 AA = DBLE( A( L ) ) 01315 * L(i,i) 01316 IF( AA.NE.ZERO ) THEN 01317 IF( SCALE.LT.AA ) THEN 01318 S = ONE + S*( SCALE / AA )**2 01319 SCALE = AA 01320 ELSE 01321 S = S + ( AA / SCALE )**2 01322 END IF 01323 END IF 01324 AA = DBLE( A( L+1 ) ) 01325 * L(k+i+1,k+i+1) 01326 IF( AA.NE.ZERO ) THEN 01327 IF( SCALE.LT.AA ) THEN 01328 S = ONE + S*( SCALE / AA )**2 01329 SCALE = AA 01330 ELSE 01331 S = S + ( AA / SCALE )**2 01332 END IF 01333 END IF 01334 L = L + LDA + 1 01335 END DO 01336 * L-> k - 1 + k*lda or L(k-1,k-1) at A(k-1,k) 01337 AA = DBLE( A( L ) ) 01338 * L(k-1,k-1) at A(k-1,k) 01339 IF( AA.NE.ZERO ) THEN 01340 IF( SCALE.LT.AA ) THEN 01341 S = ONE + S*( SCALE / AA )**2 01342 SCALE = AA 01343 ELSE 01344 S = S + ( AA / SCALE )**2 01345 END IF 01346 END IF 01347 END IF 01348 END IF 01349 END IF 01350 VALUE = SCALE*SQRT( S ) 01351 END IF 01352 * 01353 ZLANHF = VALUE 01354 RETURN 01355 * 01356 * End of ZLANHF 01357 * 01358 END
1.7.3 
