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