 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 DOUBLE PRECISION FUNCTION DLANSF( 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 2011 -- 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 A( 0: * ), WORK( 0: * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * DLANSF returns the value of the one norm, or the Frobenius norm, or 00023 * the infinity norm, or the element of largest absolute value of a 00024 * real symmetric matrix A in RFP format. 00025 * 00026 * Description 00027 * =========== 00028 * 00029 * DLANSF returns the value 00030 * 00031 * DLANSF = ( max(abs(A(i,j))), NORM = 'M' or 'm' 00032 * ( 00033 * ( norm1(A), NORM = '1', 'O' or 'o' 00034 * ( 00035 * ( normI(A), NORM = 'I' or 'i' 00036 * ( 00037 * ( normF(A), NORM = 'F', 'f', 'E' or 'e' 00038 * 00039 * where norm1 denotes the one norm of a matrix (maximum column sum), 00040 * normI denotes the infinity norm of a matrix (maximum row sum) and 00041 * normF denotes the Frobenius norm of a matrix (square root of sum of 00042 * squares). Note that max(abs(A(i,j))) is not a matrix norm. 00043 * 00044 * Arguments 00045 * ========= 00046 * 00047 * NORM (input) CHARACTER*1 00048 * Specifies the value to be returned in DLANSF as described 00049 * above. 00050 * 00051 * TRANSR (input) CHARACTER*1 00052 * Specifies whether the RFP format of A is normal or 00053 * transposed format. 00054 * = 'N': RFP format is Normal; 00055 * = 'T': RFP format is Transpose. 00056 * 00057 * UPLO (input) CHARACTER*1 00058 * On entry, UPLO specifies whether the RFP matrix A came from 00059 * an upper or lower triangular matrix as follows: 00060 * = 'U': RFP A came from an upper triangular matrix; 00061 * = 'L': RFP A came from a lower triangular matrix. 00062 * 00063 * N (input) INTEGER 00064 * The order of the matrix A. N >= 0. When N = 0, DLANSF is 00065 * set to zero. 00066 * 00067 * A (input) DOUBLE PRECISION array, dimension ( N*(N+1)/2 ); 00068 * On entry, the upper (if UPLO = 'U') or lower (if UPLO = 'L') 00069 * part of the symmetric matrix A stored in RFP format. See the 00070 * "Notes" below for more details. 00071 * Unchanged on exit. 00072 * 00073 * WORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)), 00074 * where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, 00075 * WORK is not referenced. 00076 * 00077 * Further Details 00078 * =============== 00079 * 00080 * We first consider Rectangular Full Packed (RFP) Format when N is 00081 * even. We give an example where N = 6. 00082 * 00083 * AP is Upper AP is Lower 00084 * 00085 * 00 01 02 03 04 05 00 00086 * 11 12 13 14 15 10 11 00087 * 22 23 24 25 20 21 22 00088 * 33 34 35 30 31 32 33 00089 * 44 45 40 41 42 43 44 00090 * 55 50 51 52 53 54 55 00091 * 00092 * 00093 * Let TRANSR = 'N'. RFP holds AP as follows: 00094 * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last 00095 * three columns of AP upper. The lower triangle A(4:6,0:2) consists of 00096 * the transpose of the first three columns of AP upper. 00097 * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first 00098 * three columns of AP lower. The upper triangle A(0:2,0:2) consists of 00099 * the transpose of the last three columns of AP lower. 00100 * This covers the case N even and TRANSR = 'N'. 00101 * 00102 * RFP A RFP A 00103 * 00104 * 03 04 05 33 43 53 00105 * 13 14 15 00 44 54 00106 * 23 24 25 10 11 55 00107 * 33 34 35 20 21 22 00108 * 00 44 45 30 31 32 00109 * 01 11 55 40 41 42 00110 * 02 12 22 50 51 52 00111 * 00112 * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the 00113 * transpose of RFP A above. One therefore gets: 00114 * 00115 * 00116 * RFP A RFP A 00117 * 00118 * 03 13 23 33 00 01 02 33 00 10 20 30 40 50 00119 * 04 14 24 34 44 11 12 43 44 11 21 31 41 51 00120 * 05 15 25 35 45 55 22 53 54 55 22 32 42 52 00121 * 00122 * 00123 * We then consider Rectangular Full Packed (RFP) Format when N is 00124 * odd. We give an example where N = 5. 00125 * 00126 * AP is Upper AP is Lower 00127 * 00128 * 00 01 02 03 04 00 00129 * 11 12 13 14 10 11 00130 * 22 23 24 20 21 22 00131 * 33 34 30 31 32 33 00132 * 44 40 41 42 43 44 00133 * 00134 * 00135 * Let TRANSR = 'N'. RFP holds AP as follows: 00136 * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last 00137 * three columns of AP upper. The lower triangle A(3:4,0:1) consists of 00138 * the transpose of the first two columns of AP upper. 00139 * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first 00140 * three columns of AP lower. The upper triangle A(0:1,1:2) consists of 00141 * the transpose of the last two columns of AP lower. 00142 * This covers the case N odd and TRANSR = 'N'. 00143 * 00144 * RFP A RFP A 00145 * 00146 * 02 03 04 00 33 43 00147 * 12 13 14 10 11 44 00148 * 22 23 24 20 21 22 00149 * 00 33 34 30 31 32 00150 * 01 11 44 40 41 42 00151 * 00152 * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the 00153 * transpose of RFP A above. One therefore gets: 00154 * 00155 * RFP A RFP A 00156 * 00157 * 02 12 22 00 01 00 10 20 30 40 50 00158 * 03 13 23 33 11 33 11 21 31 41 51 00159 * 04 14 24 34 44 43 44 22 32 42 52 00160 * 00161 * Reference 00162 * ========= 00163 * 00164 * ===================================================================== 00165 * 00166 * .. Parameters .. 00167 DOUBLE PRECISION ONE, ZERO 00168 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00169 * .. 00170 * .. Local Scalars .. 00171 INTEGER I, J, IFM, ILU, NOE, N1, K, L, LDA 00172 DOUBLE PRECISION SCALE, S, VALUE, AA 00173 * .. 00174 * .. External Functions .. 00175 LOGICAL LSAME 00176 INTEGER IDAMAX 00177 EXTERNAL LSAME, IDAMAX 00178 * .. 00179 * .. External Subroutines .. 00180 EXTERNAL DLASSQ 00181 * .. 00182 * .. Intrinsic Functions .. 00183 INTRINSIC ABS, MAX, SQRT 00184 * .. 00185 * .. Executable Statements .. 00186 * 00187 IF( N.EQ.0 ) THEN 00188 DLANSF = ZERO 00189 RETURN 00190 END IF 00191 * 00192 * set noe = 1 if n is odd. if n is even set noe=0 00193 * 00194 NOE = 1 00195 IF( MOD( N, 2 ).EQ.0 ) 00196 $ NOE = 0 00197 * 00198 * set ifm = 0 when form='T or 't' and 1 otherwise 00199 * 00200 IFM = 1 00201 IF( LSAME( TRANSR, 'T' ) ) 00202 $ IFM = 0 00203 * 00204 * set ilu = 0 when uplo='U or 'u' and 1 otherwise 00205 * 00206 ILU = 1 00207 IF( LSAME( UPLO, 'U' ) ) 00208 $ ILU = 0 00209 * 00210 * set lda = (n+1)/2 when ifm = 0 00211 * set lda = n when ifm = 1 and noe = 1 00212 * set lda = n+1 when ifm = 1 and noe = 0 00213 * 00214 IF( IFM.EQ.1 ) THEN 00215 IF( NOE.EQ.1 ) THEN 00216 LDA = N 00217 ELSE 00218 * noe=0 00219 LDA = N + 1 00220 END IF 00221 ELSE 00222 * ifm=0 00223 LDA = ( N+1 ) / 2 00224 END IF 00225 * 00226 IF( LSAME( NORM, 'M' ) ) THEN 00227 * 00228 * Find max(abs(A(i,j))). 00229 * 00230 K = ( N+1 ) / 2 00231 VALUE = ZERO 00232 IF( NOE.EQ.1 ) THEN 00233 * n is odd 00234 IF( IFM.EQ.1 ) THEN 00235 * A is n by k 00236 DO J = 0, K - 1 00237 DO I = 0, N - 1 00238 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00239 END DO 00240 END DO 00241 ELSE 00242 * xpose case; A is k by n 00243 DO J = 0, N - 1 00244 DO I = 0, K - 1 00245 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00246 END DO 00247 END DO 00248 END IF 00249 ELSE 00250 * n is even 00251 IF( IFM.EQ.1 ) THEN 00252 * A is n+1 by k 00253 DO J = 0, K - 1 00254 DO I = 0, N 00255 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00256 END DO 00257 END DO 00258 ELSE 00259 * xpose case; A is k by n+1 00260 DO J = 0, N 00261 DO I = 0, K - 1 00262 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00263 END DO 00264 END DO 00265 END IF 00266 END IF 00267 ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR. 00268 $ ( NORM.EQ.'1' ) ) THEN 00269 * 00270 * Find normI(A) ( = norm1(A), since A is symmetric). 00271 * 00272 IF( IFM.EQ.1 ) THEN 00273 K = N / 2 00274 IF( NOE.EQ.1 ) THEN 00275 * n is odd 00276 IF( ILU.EQ.0 ) THEN 00277 DO I = 0, K - 1 00278 WORK( I ) = ZERO 00279 END DO 00280 DO J = 0, K 00281 S = ZERO 00282 DO I = 0, K + J - 1 00283 AA = ABS( A( I+J*LDA ) ) 00284 * -> A(i,j+k) 00285 S = S + AA 00286 WORK( I ) = WORK( I ) + AA 00287 END DO 00288 AA = ABS( A( I+J*LDA ) ) 00289 * -> A(j+k,j+k) 00290 WORK( J+K ) = S + AA 00291 IF( I.EQ.K+K ) 00292 $ GO TO 10 00293 I = I + 1 00294 AA = ABS( A( I+J*LDA ) ) 00295 * -> A(j,j) 00296 WORK( J ) = WORK( J ) + AA 00297 S = ZERO 00298 DO L = J + 1, K - 1 00299 I = I + 1 00300 AA = ABS( A( I+J*LDA ) ) 00301 * -> A(l,j) 00302 S = S + AA 00303 WORK( L ) = WORK( L ) + AA 00304 END DO 00305 WORK( J ) = WORK( J ) + S 00306 END DO 00307 10 CONTINUE 00308 I = IDAMAX( N, WORK, 1 ) 00309 VALUE = WORK( I-1 ) 00310 ELSE 00311 * ilu = 1 00312 K = K + 1 00313 * k=(n+1)/2 for n odd and ilu=1 00314 DO I = K, N - 1 00315 WORK( I ) = ZERO 00316 END DO 00317 DO J = K - 1, 0, -1 00318 S = ZERO 00319 DO I = 0, J - 2 00320 AA = ABS( A( I+J*LDA ) ) 00321 * -> A(j+k,i+k) 00322 S = S + AA 00323 WORK( I+K ) = WORK( I+K ) + AA 00324 END DO 00325 IF( J.GT.0 ) THEN 00326 AA = ABS( A( I+J*LDA ) ) 00327 * -> A(j+k,j+k) 00328 S = S + AA 00329 WORK( I+K ) = WORK( I+K ) + S 00330 * i=j 00331 I = I + 1 00332 END IF 00333 AA = ABS( A( I+J*LDA ) ) 00334 * -> A(j,j) 00335 WORK( J ) = AA 00336 S = ZERO 00337 DO L = J + 1, N - 1 00338 I = I + 1 00339 AA = ABS( A( I+J*LDA ) ) 00340 * -> A(l,j) 00341 S = S + AA 00342 WORK( L ) = WORK( L ) + AA 00343 END DO 00344 WORK( J ) = WORK( J ) + S 00345 END DO 00346 I = IDAMAX( N, WORK, 1 ) 00347 VALUE = WORK( I-1 ) 00348 END IF 00349 ELSE 00350 * n is even 00351 IF( ILU.EQ.0 ) THEN 00352 DO I = 0, K - 1 00353 WORK( I ) = ZERO 00354 END DO 00355 DO J = 0, K - 1 00356 S = ZERO 00357 DO I = 0, K + J - 1 00358 AA = ABS( A( I+J*LDA ) ) 00359 * -> A(i,j+k) 00360 S = S + AA 00361 WORK( I ) = WORK( I ) + AA 00362 END DO 00363 AA = ABS( A( I+J*LDA ) ) 00364 * -> A(j+k,j+k) 00365 WORK( J+K ) = S + AA 00366 I = I + 1 00367 AA = ABS( A( I+J*LDA ) ) 00368 * -> A(j,j) 00369 WORK( J ) = WORK( J ) + AA 00370 S = ZERO 00371 DO L = J + 1, K - 1 00372 I = I + 1 00373 AA = ABS( A( I+J*LDA ) ) 00374 * -> A(l,j) 00375 S = S + AA 00376 WORK( L ) = WORK( L ) + AA 00377 END DO 00378 WORK( J ) = WORK( J ) + S 00379 END DO 00380 I = IDAMAX( N, WORK, 1 ) 00381 VALUE = WORK( I-1 ) 00382 ELSE 00383 * ilu = 1 00384 DO I = K, N - 1 00385 WORK( I ) = ZERO 00386 END DO 00387 DO J = K - 1, 0, -1 00388 S = ZERO 00389 DO I = 0, J - 1 00390 AA = ABS( A( I+J*LDA ) ) 00391 * -> A(j+k,i+k) 00392 S = S + AA 00393 WORK( I+K ) = WORK( I+K ) + AA 00394 END DO 00395 AA = ABS( A( I+J*LDA ) ) 00396 * -> A(j+k,j+k) 00397 S = S + AA 00398 WORK( I+K ) = WORK( I+K ) + S 00399 * i=j 00400 I = I + 1 00401 AA = ABS( A( I+J*LDA ) ) 00402 * -> A(j,j) 00403 WORK( J ) = AA 00404 S = ZERO 00405 DO L = J + 1, N - 1 00406 I = I + 1 00407 AA = ABS( A( I+J*LDA ) ) 00408 * -> A(l,j) 00409 S = S + AA 00410 WORK( L ) = WORK( L ) + AA 00411 END DO 00412 WORK( J ) = WORK( J ) + S 00413 END DO 00414 I = IDAMAX( N, WORK, 1 ) 00415 VALUE = WORK( I-1 ) 00416 END IF 00417 END IF 00418 ELSE 00419 * ifm=0 00420 K = N / 2 00421 IF( NOE.EQ.1 ) THEN 00422 * n is odd 00423 IF( ILU.EQ.0 ) THEN 00424 N1 = K 00425 * n/2 00426 K = K + 1 00427 * k is the row size and lda 00428 DO I = N1, N - 1 00429 WORK( I ) = ZERO 00430 END DO 00431 DO J = 0, N1 - 1 00432 S = ZERO 00433 DO I = 0, K - 1 00434 AA = ABS( A( I+J*LDA ) ) 00435 * A(j,n1+i) 00436 WORK( I+N1 ) = WORK( I+N1 ) + AA 00437 S = S + AA 00438 END DO 00439 WORK( J ) = S 00440 END DO 00441 * j=n1=k-1 is special 00442 S = ABS( A( 0+J*LDA ) ) 00443 * A(k-1,k-1) 00444 DO I = 1, K - 1 00445 AA = ABS( A( I+J*LDA ) ) 00446 * A(k-1,i+n1) 00447 WORK( I+N1 ) = WORK( I+N1 ) + AA 00448 S = S + AA 00449 END DO 00450 WORK( J ) = WORK( J ) + S 00451 DO J = K, N - 1 00452 S = ZERO 00453 DO I = 0, J - K - 1 00454 AA = ABS( A( I+J*LDA ) ) 00455 * A(i,j-k) 00456 WORK( I ) = WORK( I ) + AA 00457 S = S + AA 00458 END DO 00459 * i=j-k 00460 AA = ABS( A( I+J*LDA ) ) 00461 * A(j-k,j-k) 00462 S = S + AA 00463 WORK( J-K ) = WORK( J-K ) + S 00464 I = I + 1 00465 S = ABS( A( I+J*LDA ) ) 00466 * A(j,j) 00467 DO L = J + 1, N - 1 00468 I = I + 1 00469 AA = ABS( A( I+J*LDA ) ) 00470 * A(j,l) 00471 WORK( L ) = WORK( L ) + AA 00472 S = S + AA 00473 END DO 00474 WORK( J ) = WORK( J ) + S 00475 END DO 00476 I = IDAMAX( N, WORK, 1 ) 00477 VALUE = WORK( I-1 ) 00478 ELSE 00479 * ilu=1 00480 K = K + 1 00481 * k=(n+1)/2 for n odd and ilu=1 00482 DO I = K, N - 1 00483 WORK( I ) = ZERO 00484 END DO 00485 DO J = 0, K - 2 00486 * process 00487 S = ZERO 00488 DO I = 0, J - 1 00489 AA = ABS( A( I+J*LDA ) ) 00490 * A(j,i) 00491 WORK( I ) = WORK( I ) + AA 00492 S = S + AA 00493 END DO 00494 AA = ABS( A( I+J*LDA ) ) 00495 * i=j so process of A(j,j) 00496 S = S + AA 00497 WORK( J ) = S 00498 * is initialised here 00499 I = I + 1 00500 * i=j process A(j+k,j+k) 00501 AA = ABS( A( I+J*LDA ) ) 00502 S = AA 00503 DO L = K + J + 1, N - 1 00504 I = I + 1 00505 AA = ABS( A( I+J*LDA ) ) 00506 * A(l,k+j) 00507 S = S + AA 00508 WORK( L ) = WORK( L ) + AA 00509 END DO 00510 WORK( K+J ) = WORK( K+J ) + S 00511 END DO 00512 * j=k-1 is special :process col A(k-1,0:k-1) 00513 S = ZERO 00514 DO I = 0, K - 2 00515 AA = ABS( A( I+J*LDA ) ) 00516 * A(k,i) 00517 WORK( I ) = WORK( I ) + AA 00518 S = S + AA 00519 END DO 00520 * i=k-1 00521 AA = ABS( A( I+J*LDA ) ) 00522 * A(k-1,k-1) 00523 S = S + AA 00524 WORK( I ) = S 00525 * done with col j=k+1 00526 DO J = K, N - 1 00527 * process col j of A = A(j,0:k-1) 00528 S = ZERO 00529 DO I = 0, K - 1 00530 AA = ABS( A( I+J*LDA ) ) 00531 * A(j,i) 00532 WORK( I ) = WORK( I ) + AA 00533 S = S + AA 00534 END DO 00535 WORK( J ) = WORK( J ) + S 00536 END DO 00537 I = IDAMAX( N, WORK, 1 ) 00538 VALUE = WORK( I-1 ) 00539 END IF 00540 ELSE 00541 * n is even 00542 IF( ILU.EQ.0 ) THEN 00543 DO I = K, N - 1 00544 WORK( I ) = ZERO 00545 END DO 00546 DO J = 0, K - 1 00547 S = ZERO 00548 DO I = 0, K - 1 00549 AA = ABS( A( I+J*LDA ) ) 00550 * A(j,i+k) 00551 WORK( I+K ) = WORK( I+K ) + AA 00552 S = S + AA 00553 END DO 00554 WORK( J ) = S 00555 END DO 00556 * j=k 00557 AA = ABS( A( 0+J*LDA ) ) 00558 * A(k,k) 00559 S = AA 00560 DO I = 1, K - 1 00561 AA = ABS( A( I+J*LDA ) ) 00562 * A(k,k+i) 00563 WORK( I+K ) = WORK( I+K ) + AA 00564 S = S + AA 00565 END DO 00566 WORK( J ) = WORK( J ) + S 00567 DO J = K + 1, N - 1 00568 S = ZERO 00569 DO I = 0, J - 2 - K 00570 AA = ABS( A( I+J*LDA ) ) 00571 * A(i,j-k-1) 00572 WORK( I ) = WORK( I ) + AA 00573 S = S + AA 00574 END DO 00575 * i=j-1-k 00576 AA = ABS( A( I+J*LDA ) ) 00577 * A(j-k-1,j-k-1) 00578 S = S + AA 00579 WORK( J-K-1 ) = WORK( J-K-1 ) + S 00580 I = I + 1 00581 AA = ABS( A( I+J*LDA ) ) 00582 * A(j,j) 00583 S = AA 00584 DO L = J + 1, N - 1 00585 I = I + 1 00586 AA = ABS( A( I+J*LDA ) ) 00587 * A(j,l) 00588 WORK( L ) = WORK( L ) + AA 00589 S = S + AA 00590 END DO 00591 WORK( J ) = WORK( J ) + S 00592 END DO 00593 * j=n 00594 S = ZERO 00595 DO I = 0, K - 2 00596 AA = ABS( A( I+J*LDA ) ) 00597 * A(i,k-1) 00598 WORK( I ) = WORK( I ) + AA 00599 S = S + AA 00600 END DO 00601 * i=k-1 00602 AA = ABS( A( I+J*LDA ) ) 00603 * A(k-1,k-1) 00604 S = S + AA 00605 WORK( I ) = WORK( I ) + S 00606 I = IDAMAX( N, WORK, 1 ) 00607 VALUE = WORK( I-1 ) 00608 ELSE 00609 * ilu=1 00610 DO I = K, N - 1 00611 WORK( I ) = ZERO 00612 END DO 00613 * j=0 is special :process col A(k:n-1,k) 00614 S = ABS( A( 0 ) ) 00615 * A(k,k) 00616 DO I = 1, K - 1 00617 AA = ABS( A( I ) ) 00618 * A(k+i,k) 00619 WORK( I+K ) = WORK( I+K ) + AA 00620 S = S + AA 00621 END DO 00622 WORK( K ) = WORK( K ) + S 00623 DO J = 1, K - 1 00624 * process 00625 S = ZERO 00626 DO I = 0, J - 2 00627 AA = ABS( A( I+J*LDA ) ) 00628 * A(j-1,i) 00629 WORK( I ) = WORK( I ) + AA 00630 S = S + AA 00631 END DO 00632 AA = ABS( A( I+J*LDA ) ) 00633 * i=j-1 so process of A(j-1,j-1) 00634 S = S + AA 00635 WORK( J-1 ) = S 00636 * is initialised here 00637 I = I + 1 00638 * i=j process A(j+k,j+k) 00639 AA = ABS( A( I+J*LDA ) ) 00640 S = AA 00641 DO L = K + J + 1, N - 1 00642 I = I + 1 00643 AA = ABS( A( I+J*LDA ) ) 00644 * A(l,k+j) 00645 S = S + AA 00646 WORK( L ) = WORK( L ) + AA 00647 END DO 00648 WORK( K+J ) = WORK( K+J ) + S 00649 END DO 00650 * j=k is special :process col A(k,0:k-1) 00651 S = ZERO 00652 DO I = 0, K - 2 00653 AA = ABS( A( I+J*LDA ) ) 00654 * A(k,i) 00655 WORK( I ) = WORK( I ) + AA 00656 S = S + AA 00657 END DO 00658 * i=k-1 00659 AA = ABS( A( I+J*LDA ) ) 00660 * A(k-1,k-1) 00661 S = S + AA 00662 WORK( I ) = S 00663 * done with col j=k+1 00664 DO J = K + 1, N 00665 * process col j-1 of A = A(j-1,0:k-1) 00666 S = ZERO 00667 DO I = 0, K - 1 00668 AA = ABS( A( I+J*LDA ) ) 00669 * A(j-1,i) 00670 WORK( I ) = WORK( I ) + AA 00671 S = S + AA 00672 END DO 00673 WORK( J-1 ) = WORK( J-1 ) + S 00674 END DO 00675 I = IDAMAX( N, WORK, 1 ) 00676 VALUE = WORK( I-1 ) 00677 END IF 00678 END IF 00679 END IF 00680 ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN 00681 * 00682 * Find normF(A). 00683 * 00684 K = ( N+1 ) / 2 00685 SCALE = ZERO 00686 S = ONE 00687 IF( NOE.EQ.1 ) THEN 00688 * n is odd 00689 IF( IFM.EQ.1 ) THEN 00690 * A is normal 00691 IF( ILU.EQ.0 ) THEN 00692 * A is upper 00693 DO J = 0, K - 3 00694 CALL DLASSQ( K-J-2, A( K+J+1+J*LDA ), 1, SCALE, S ) 00695 * L at A(k,0) 00696 END DO 00697 DO J = 0, K - 1 00698 CALL DLASSQ( K+J-1, A( 0+J*LDA ), 1, SCALE, S ) 00699 * trap U at A(0,0) 00700 END DO 00701 S = S + S 00702 * double s for the off diagonal elements 00703 CALL DLASSQ( K-1, A( K ), LDA+1, SCALE, S ) 00704 * tri L at A(k,0) 00705 CALL DLASSQ( K, A( K-1 ), LDA+1, SCALE, S ) 00706 * tri U at A(k-1,0) 00707 ELSE 00708 * ilu=1 & A is lower 00709 DO J = 0, K - 1 00710 CALL DLASSQ( N-J-1, A( J+1+J*LDA ), 1, SCALE, S ) 00711 * trap L at A(0,0) 00712 END DO 00713 DO J = 0, K - 2 00714 CALL DLASSQ( J, A( 0+( 1+J )*LDA ), 1, SCALE, S ) 00715 * U at A(0,1) 00716 END DO 00717 S = S + S 00718 * double s for the off diagonal elements 00719 CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S ) 00720 * tri L at A(0,0) 00721 CALL DLASSQ( K-1, A( 0+LDA ), LDA+1, SCALE, S ) 00722 * tri U at A(0,1) 00723 END IF 00724 ELSE 00725 * A is xpose 00726 IF( ILU.EQ.0 ) THEN 00727 * A**T is upper 00728 DO J = 1, K - 2 00729 CALL DLASSQ( J, A( 0+( K+J )*LDA ), 1, SCALE, S ) 00730 * U at A(0,k) 00731 END DO 00732 DO J = 0, K - 2 00733 CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 00734 * k by k-1 rect. at A(0,0) 00735 END DO 00736 DO J = 0, K - 2 00737 CALL DLASSQ( K-J-1, A( J+1+( J+K-1 )*LDA ), 1, 00738 $ SCALE, S ) 00739 * L at A(0,k-1) 00740 END DO 00741 S = S + S 00742 * double s for the off diagonal elements 00743 CALL DLASSQ( K-1, A( 0+K*LDA ), LDA+1, SCALE, S ) 00744 * tri U at A(0,k) 00745 CALL DLASSQ( K, A( 0+( K-1 )*LDA ), LDA+1, SCALE, S ) 00746 * tri L at A(0,k-1) 00747 ELSE 00748 * A**T is lower 00749 DO J = 1, K - 1 00750 CALL DLASSQ( J, A( 0+J*LDA ), 1, SCALE, S ) 00751 * U at A(0,0) 00752 END DO 00753 DO J = K, N - 1 00754 CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 00755 * k by k-1 rect. at A(0,k) 00756 END DO 00757 DO J = 0, K - 3 00758 CALL DLASSQ( K-J-2, A( J+2+J*LDA ), 1, SCALE, S ) 00759 * L at A(1,0) 00760 END DO 00761 S = S + S 00762 * double s for the off diagonal elements 00763 CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S ) 00764 * tri U at A(0,0) 00765 CALL DLASSQ( K-1, A( 1 ), LDA+1, SCALE, S ) 00766 * tri L at A(1,0) 00767 END IF 00768 END IF 00769 ELSE 00770 * n is even 00771 IF( IFM.EQ.1 ) THEN 00772 * A is normal 00773 IF( ILU.EQ.0 ) THEN 00774 * A is upper 00775 DO J = 0, K - 2 00776 CALL DLASSQ( K-J-1, A( K+J+2+J*LDA ), 1, SCALE, S ) 00777 * L at A(k+1,0) 00778 END DO 00779 DO J = 0, K - 1 00780 CALL DLASSQ( K+J, A( 0+J*LDA ), 1, SCALE, S ) 00781 * trap U at A(0,0) 00782 END DO 00783 S = S + S 00784 * double s for the off diagonal elements 00785 CALL DLASSQ( K, A( K+1 ), LDA+1, SCALE, S ) 00786 * tri L at A(k+1,0) 00787 CALL DLASSQ( K, A( K ), LDA+1, SCALE, S ) 00788 * tri U at A(k,0) 00789 ELSE 00790 * ilu=1 & A is lower 00791 DO J = 0, K - 1 00792 CALL DLASSQ( N-J-1, A( J+2+J*LDA ), 1, SCALE, S ) 00793 * trap L at A(1,0) 00794 END DO 00795 DO J = 1, K - 1 00796 CALL DLASSQ( J, A( 0+J*LDA ), 1, SCALE, S ) 00797 * U at A(0,0) 00798 END DO 00799 S = S + S 00800 * double s for the off diagonal elements 00801 CALL DLASSQ( K, A( 1 ), LDA+1, SCALE, S ) 00802 * tri L at A(1,0) 00803 CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S ) 00804 * tri U at A(0,0) 00805 END IF 00806 ELSE 00807 * A is xpose 00808 IF( ILU.EQ.0 ) THEN 00809 * A**T is upper 00810 DO J = 1, K - 1 00811 CALL DLASSQ( J, A( 0+( K+1+J )*LDA ), 1, SCALE, S ) 00812 * U at A(0,k+1) 00813 END DO 00814 DO J = 0, K - 1 00815 CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 00816 * k by k rect. at A(0,0) 00817 END DO 00818 DO J = 0, K - 2 00819 CALL DLASSQ( K-J-1, A( J+1+( J+K )*LDA ), 1, SCALE, 00820 $ S ) 00821 * L at A(0,k) 00822 END DO 00823 S = S + S 00824 * double s for the off diagonal elements 00825 CALL DLASSQ( K, A( 0+( K+1 )*LDA ), LDA+1, SCALE, S ) 00826 * tri U at A(0,k+1) 00827 CALL DLASSQ( K, A( 0+K*LDA ), LDA+1, SCALE, S ) 00828 * tri L at A(0,k) 00829 ELSE 00830 * A**T is lower 00831 DO J = 1, K - 1 00832 CALL DLASSQ( J, A( 0+( J+1 )*LDA ), 1, SCALE, S ) 00833 * U at A(0,1) 00834 END DO 00835 DO J = K + 1, N 00836 CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 00837 * k by k rect. at A(0,k+1) 00838 END DO 00839 DO J = 0, K - 2 00840 CALL DLASSQ( K-J-1, A( J+1+J*LDA ), 1, SCALE, S ) 00841 * L at A(0,0) 00842 END DO 00843 S = S + S 00844 * double s for the off diagonal elements 00845 CALL DLASSQ( K, A( LDA ), LDA+1, SCALE, S ) 00846 * tri L at A(0,1) 00847 CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S ) 00848 * tri U at A(0,0) 00849 END IF 00850 END IF 00851 END IF 00852 VALUE = SCALE*SQRT( S ) 00853 END IF 00854 * 00855 DLANSF = VALUE 00856 RETURN 00857 * 00858 * End of DLANSF 00859 * 00860 END
1.7.3 
