 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 REAL FUNCTION SLANSF( 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 REAL A( 0: * ), WORK( 0: * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * SLANSF 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 * SLANSF returns the value 00030 * 00031 * SLANSF = ( 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 SLANSF 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, SLANSF is 00065 * set to zero. 00066 * 00067 * A (input) REAL 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) REAL 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 * .. 00167 * .. Parameters .. 00168 REAL ONE, ZERO 00169 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00170 * .. 00171 * .. Local Scalars .. 00172 INTEGER I, J, IFM, ILU, NOE, N1, K, L, LDA 00173 REAL SCALE, S, VALUE, AA 00174 * .. 00175 * .. External Functions .. 00176 LOGICAL LSAME 00177 INTEGER ISAMAX 00178 EXTERNAL LSAME, ISAMAX 00179 * .. 00180 * .. External Subroutines .. 00181 EXTERNAL SLASSQ 00182 * .. 00183 * .. Intrinsic Functions .. 00184 INTRINSIC ABS, MAX, SQRT 00185 * .. 00186 * .. Executable Statements .. 00187 * 00188 IF( N.EQ.0 ) THEN 00189 SLANSF = ZERO 00190 RETURN 00191 END IF 00192 * 00193 * set noe = 1 if n is odd. if n is even set noe=0 00194 * 00195 NOE = 1 00196 IF( MOD( N, 2 ).EQ.0 ) 00197 $ NOE = 0 00198 * 00199 * set ifm = 0 when form='T or 't' and 1 otherwise 00200 * 00201 IFM = 1 00202 IF( LSAME( TRANSR, 'T' ) ) 00203 $ IFM = 0 00204 * 00205 * set ilu = 0 when uplo='U or 'u' and 1 otherwise 00206 * 00207 ILU = 1 00208 IF( LSAME( UPLO, 'U' ) ) 00209 $ ILU = 0 00210 * 00211 * set lda = (n+1)/2 when ifm = 0 00212 * set lda = n when ifm = 1 and noe = 1 00213 * set lda = n+1 when ifm = 1 and noe = 0 00214 * 00215 IF( IFM.EQ.1 ) THEN 00216 IF( NOE.EQ.1 ) THEN 00217 LDA = N 00218 ELSE 00219 * noe=0 00220 LDA = N + 1 00221 END IF 00222 ELSE 00223 * ifm=0 00224 LDA = ( N+1 ) / 2 00225 END IF 00226 * 00227 IF( LSAME( NORM, 'M' ) ) THEN 00228 * 00229 * Find max(abs(A(i,j))). 00230 * 00231 K = ( N+1 ) / 2 00232 VALUE = ZERO 00233 IF( NOE.EQ.1 ) THEN 00234 * n is odd 00235 IF( IFM.EQ.1 ) THEN 00236 * A is n by k 00237 DO J = 0, K - 1 00238 DO I = 0, N - 1 00239 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00240 END DO 00241 END DO 00242 ELSE 00243 * xpose case; A is k by n 00244 DO J = 0, N - 1 00245 DO I = 0, K - 1 00246 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00247 END DO 00248 END DO 00249 END IF 00250 ELSE 00251 * n is even 00252 IF( IFM.EQ.1 ) THEN 00253 * A is n+1 by k 00254 DO J = 0, K - 1 00255 DO I = 0, N 00256 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00257 END DO 00258 END DO 00259 ELSE 00260 * xpose case; A is k by n+1 00261 DO J = 0, N 00262 DO I = 0, K - 1 00263 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00264 END DO 00265 END DO 00266 END IF 00267 END IF 00268 ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR. 00269 $ ( NORM.EQ.'1' ) ) THEN 00270 * 00271 * Find normI(A) ( = norm1(A), since A is symmetric). 00272 * 00273 IF( IFM.EQ.1 ) THEN 00274 K = N / 2 00275 IF( NOE.EQ.1 ) THEN 00276 * n is odd 00277 IF( ILU.EQ.0 ) THEN 00278 DO I = 0, K - 1 00279 WORK( I ) = ZERO 00280 END DO 00281 DO J = 0, K 00282 S = ZERO 00283 DO I = 0, K + J - 1 00284 AA = ABS( A( I+J*LDA ) ) 00285 * -> A(i,j+k) 00286 S = S + AA 00287 WORK( I ) = WORK( I ) + AA 00288 END DO 00289 AA = ABS( A( I+J*LDA ) ) 00290 * -> A(j+k,j+k) 00291 WORK( J+K ) = S + AA 00292 IF( I.EQ.K+K ) 00293 $ GO TO 10 00294 I = I + 1 00295 AA = ABS( A( I+J*LDA ) ) 00296 * -> A(j,j) 00297 WORK( J ) = WORK( J ) + AA 00298 S = ZERO 00299 DO L = J + 1, K - 1 00300 I = I + 1 00301 AA = ABS( A( I+J*LDA ) ) 00302 * -> A(l,j) 00303 S = S + AA 00304 WORK( L ) = WORK( L ) + AA 00305 END DO 00306 WORK( J ) = WORK( J ) + S 00307 END DO 00308 10 CONTINUE 00309 I = ISAMAX( N, WORK, 1 ) 00310 VALUE = WORK( I-1 ) 00311 ELSE 00312 * ilu = 1 00313 K = K + 1 00314 * k=(n+1)/2 for n odd and ilu=1 00315 DO I = K, N - 1 00316 WORK( I ) = ZERO 00317 END DO 00318 DO J = K - 1, 0, -1 00319 S = ZERO 00320 DO I = 0, J - 2 00321 AA = ABS( A( I+J*LDA ) ) 00322 * -> A(j+k,i+k) 00323 S = S + AA 00324 WORK( I+K ) = WORK( I+K ) + AA 00325 END DO 00326 IF( J.GT.0 ) THEN 00327 AA = ABS( A( I+J*LDA ) ) 00328 * -> A(j+k,j+k) 00329 S = S + AA 00330 WORK( I+K ) = WORK( I+K ) + S 00331 * i=j 00332 I = I + 1 00333 END IF 00334 AA = ABS( A( I+J*LDA ) ) 00335 * -> A(j,j) 00336 WORK( J ) = AA 00337 S = ZERO 00338 DO L = J + 1, N - 1 00339 I = I + 1 00340 AA = ABS( A( I+J*LDA ) ) 00341 * -> A(l,j) 00342 S = S + AA 00343 WORK( L ) = WORK( L ) + AA 00344 END DO 00345 WORK( J ) = WORK( J ) + S 00346 END DO 00347 I = ISAMAX( N, WORK, 1 ) 00348 VALUE = WORK( I-1 ) 00349 END IF 00350 ELSE 00351 * n is even 00352 IF( ILU.EQ.0 ) THEN 00353 DO I = 0, K - 1 00354 WORK( I ) = ZERO 00355 END DO 00356 DO J = 0, K - 1 00357 S = ZERO 00358 DO I = 0, K + J - 1 00359 AA = ABS( A( I+J*LDA ) ) 00360 * -> A(i,j+k) 00361 S = S + AA 00362 WORK( I ) = WORK( I ) + AA 00363 END DO 00364 AA = ABS( A( I+J*LDA ) ) 00365 * -> A(j+k,j+k) 00366 WORK( J+K ) = S + AA 00367 I = I + 1 00368 AA = ABS( A( I+J*LDA ) ) 00369 * -> A(j,j) 00370 WORK( J ) = WORK( J ) + AA 00371 S = ZERO 00372 DO L = J + 1, K - 1 00373 I = I + 1 00374 AA = ABS( A( I+J*LDA ) ) 00375 * -> A(l,j) 00376 S = S + AA 00377 WORK( L ) = WORK( L ) + AA 00378 END DO 00379 WORK( J ) = WORK( J ) + S 00380 END DO 00381 I = ISAMAX( N, WORK, 1 ) 00382 VALUE = WORK( I-1 ) 00383 ELSE 00384 * ilu = 1 00385 DO I = K, N - 1 00386 WORK( I ) = ZERO 00387 END DO 00388 DO J = K - 1, 0, -1 00389 S = ZERO 00390 DO I = 0, J - 1 00391 AA = ABS( A( I+J*LDA ) ) 00392 * -> A(j+k,i+k) 00393 S = S + AA 00394 WORK( I+K ) = WORK( I+K ) + AA 00395 END DO 00396 AA = ABS( A( I+J*LDA ) ) 00397 * -> A(j+k,j+k) 00398 S = S + AA 00399 WORK( I+K ) = WORK( I+K ) + S 00400 * i=j 00401 I = I + 1 00402 AA = ABS( A( I+J*LDA ) ) 00403 * -> A(j,j) 00404 WORK( J ) = AA 00405 S = ZERO 00406 DO L = J + 1, N - 1 00407 I = I + 1 00408 AA = ABS( A( I+J*LDA ) ) 00409 * -> A(l,j) 00410 S = S + AA 00411 WORK( L ) = WORK( L ) + AA 00412 END DO 00413 WORK( J ) = WORK( J ) + S 00414 END DO 00415 I = ISAMAX( N, WORK, 1 ) 00416 VALUE = WORK( I-1 ) 00417 END IF 00418 END IF 00419 ELSE 00420 * ifm=0 00421 K = N / 2 00422 IF( NOE.EQ.1 ) THEN 00423 * n is odd 00424 IF( ILU.EQ.0 ) THEN 00425 N1 = K 00426 * n/2 00427 K = K + 1 00428 * k is the row size and lda 00429 DO I = N1, N - 1 00430 WORK( I ) = ZERO 00431 END DO 00432 DO J = 0, N1 - 1 00433 S = ZERO 00434 DO I = 0, K - 1 00435 AA = ABS( A( I+J*LDA ) ) 00436 * A(j,n1+i) 00437 WORK( I+N1 ) = WORK( I+N1 ) + AA 00438 S = S + AA 00439 END DO 00440 WORK( J ) = S 00441 END DO 00442 * j=n1=k-1 is special 00443 S = ABS( A( 0+J*LDA ) ) 00444 * A(k-1,k-1) 00445 DO I = 1, K - 1 00446 AA = ABS( A( I+J*LDA ) ) 00447 * A(k-1,i+n1) 00448 WORK( I+N1 ) = WORK( I+N1 ) + AA 00449 S = S + AA 00450 END DO 00451 WORK( J ) = WORK( J ) + S 00452 DO J = K, N - 1 00453 S = ZERO 00454 DO I = 0, J - K - 1 00455 AA = ABS( A( I+J*LDA ) ) 00456 * A(i,j-k) 00457 WORK( I ) = WORK( I ) + AA 00458 S = S + AA 00459 END DO 00460 * i=j-k 00461 AA = ABS( A( I+J*LDA ) ) 00462 * A(j-k,j-k) 00463 S = S + AA 00464 WORK( J-K ) = WORK( J-K ) + S 00465 I = I + 1 00466 S = ABS( A( I+J*LDA ) ) 00467 * A(j,j) 00468 DO L = J + 1, N - 1 00469 I = I + 1 00470 AA = ABS( A( I+J*LDA ) ) 00471 * A(j,l) 00472 WORK( L ) = WORK( L ) + AA 00473 S = S + AA 00474 END DO 00475 WORK( J ) = WORK( J ) + S 00476 END DO 00477 I = ISAMAX( N, WORK, 1 ) 00478 VALUE = WORK( I-1 ) 00479 ELSE 00480 * ilu=1 00481 K = K + 1 00482 * k=(n+1)/2 for n odd and ilu=1 00483 DO I = K, N - 1 00484 WORK( I ) = ZERO 00485 END DO 00486 DO J = 0, K - 2 00487 * process 00488 S = ZERO 00489 DO I = 0, J - 1 00490 AA = ABS( A( I+J*LDA ) ) 00491 * A(j,i) 00492 WORK( I ) = WORK( I ) + AA 00493 S = S + AA 00494 END DO 00495 AA = ABS( A( I+J*LDA ) ) 00496 * i=j so process of A(j,j) 00497 S = S + AA 00498 WORK( J ) = S 00499 * is initialised here 00500 I = I + 1 00501 * i=j process A(j+k,j+k) 00502 AA = ABS( A( I+J*LDA ) ) 00503 S = AA 00504 DO L = K + J + 1, N - 1 00505 I = I + 1 00506 AA = ABS( A( I+J*LDA ) ) 00507 * A(l,k+j) 00508 S = S + AA 00509 WORK( L ) = WORK( L ) + AA 00510 END DO 00511 WORK( K+J ) = WORK( K+J ) + S 00512 END DO 00513 * j=k-1 is special :process col A(k-1,0:k-1) 00514 S = ZERO 00515 DO I = 0, K - 2 00516 AA = ABS( A( I+J*LDA ) ) 00517 * A(k,i) 00518 WORK( I ) = WORK( I ) + AA 00519 S = S + AA 00520 END DO 00521 * i=k-1 00522 AA = ABS( A( I+J*LDA ) ) 00523 * A(k-1,k-1) 00524 S = S + AA 00525 WORK( I ) = S 00526 * done with col j=k+1 00527 DO J = K, N - 1 00528 * process col j of A = A(j,0:k-1) 00529 S = ZERO 00530 DO I = 0, K - 1 00531 AA = ABS( A( I+J*LDA ) ) 00532 * A(j,i) 00533 WORK( I ) = WORK( I ) + AA 00534 S = S + AA 00535 END DO 00536 WORK( J ) = WORK( J ) + S 00537 END DO 00538 I = ISAMAX( N, WORK, 1 ) 00539 VALUE = WORK( I-1 ) 00540 END IF 00541 ELSE 00542 * n is even 00543 IF( ILU.EQ.0 ) THEN 00544 DO I = K, N - 1 00545 WORK( I ) = ZERO 00546 END DO 00547 DO J = 0, K - 1 00548 S = ZERO 00549 DO I = 0, K - 1 00550 AA = ABS( A( I+J*LDA ) ) 00551 * A(j,i+k) 00552 WORK( I+K ) = WORK( I+K ) + AA 00553 S = S + AA 00554 END DO 00555 WORK( J ) = S 00556 END DO 00557 * j=k 00558 AA = ABS( A( 0+J*LDA ) ) 00559 * A(k,k) 00560 S = AA 00561 DO I = 1, K - 1 00562 AA = ABS( A( I+J*LDA ) ) 00563 * A(k,k+i) 00564 WORK( I+K ) = WORK( I+K ) + AA 00565 S = S + AA 00566 END DO 00567 WORK( J ) = WORK( J ) + S 00568 DO J = K + 1, N - 1 00569 S = ZERO 00570 DO I = 0, J - 2 - K 00571 AA = ABS( A( I+J*LDA ) ) 00572 * A(i,j-k-1) 00573 WORK( I ) = WORK( I ) + AA 00574 S = S + AA 00575 END DO 00576 * i=j-1-k 00577 AA = ABS( A( I+J*LDA ) ) 00578 * A(j-k-1,j-k-1) 00579 S = S + AA 00580 WORK( J-K-1 ) = WORK( J-K-1 ) + S 00581 I = I + 1 00582 AA = ABS( A( I+J*LDA ) ) 00583 * A(j,j) 00584 S = AA 00585 DO L = J + 1, N - 1 00586 I = I + 1 00587 AA = ABS( A( I+J*LDA ) ) 00588 * A(j,l) 00589 WORK( L ) = WORK( L ) + AA 00590 S = S + AA 00591 END DO 00592 WORK( J ) = WORK( J ) + S 00593 END DO 00594 * j=n 00595 S = ZERO 00596 DO I = 0, K - 2 00597 AA = ABS( A( I+J*LDA ) ) 00598 * A(i,k-1) 00599 WORK( I ) = WORK( I ) + AA 00600 S = S + AA 00601 END DO 00602 * i=k-1 00603 AA = ABS( A( I+J*LDA ) ) 00604 * A(k-1,k-1) 00605 S = S + AA 00606 WORK( I ) = WORK( I ) + S 00607 I = ISAMAX( N, WORK, 1 ) 00608 VALUE = WORK( I-1 ) 00609 ELSE 00610 * ilu=1 00611 DO I = K, N - 1 00612 WORK( I ) = ZERO 00613 END DO 00614 * j=0 is special :process col A(k:n-1,k) 00615 S = ABS( A( 0 ) ) 00616 * A(k,k) 00617 DO I = 1, K - 1 00618 AA = ABS( A( I ) ) 00619 * A(k+i,k) 00620 WORK( I+K ) = WORK( I+K ) + AA 00621 S = S + AA 00622 END DO 00623 WORK( K ) = WORK( K ) + S 00624 DO J = 1, K - 1 00625 * process 00626 S = ZERO 00627 DO I = 0, J - 2 00628 AA = ABS( A( I+J*LDA ) ) 00629 * A(j-1,i) 00630 WORK( I ) = WORK( I ) + AA 00631 S = S + AA 00632 END DO 00633 AA = ABS( A( I+J*LDA ) ) 00634 * i=j-1 so process of A(j-1,j-1) 00635 S = S + AA 00636 WORK( J-1 ) = S 00637 * is initialised here 00638 I = I + 1 00639 * i=j process A(j+k,j+k) 00640 AA = ABS( A( I+J*LDA ) ) 00641 S = AA 00642 DO L = K + J + 1, N - 1 00643 I = I + 1 00644 AA = ABS( A( I+J*LDA ) ) 00645 * A(l,k+j) 00646 S = S + AA 00647 WORK( L ) = WORK( L ) + AA 00648 END DO 00649 WORK( K+J ) = WORK( K+J ) + S 00650 END DO 00651 * j=k is special :process col A(k,0:k-1) 00652 S = ZERO 00653 DO I = 0, K - 2 00654 AA = ABS( A( I+J*LDA ) ) 00655 * A(k,i) 00656 WORK( I ) = WORK( I ) + AA 00657 S = S + AA 00658 END DO 00659 * i=k-1 00660 AA = ABS( A( I+J*LDA ) ) 00661 * A(k-1,k-1) 00662 S = S + AA 00663 WORK( I ) = S 00664 * done with col j=k+1 00665 DO J = K + 1, N 00666 * process col j-1 of A = A(j-1,0:k-1) 00667 S = ZERO 00668 DO I = 0, K - 1 00669 AA = ABS( A( I+J*LDA ) ) 00670 * A(j-1,i) 00671 WORK( I ) = WORK( I ) + AA 00672 S = S + AA 00673 END DO 00674 WORK( J-1 ) = WORK( J-1 ) + S 00675 END DO 00676 I = ISAMAX( N, WORK, 1 ) 00677 VALUE = WORK( I-1 ) 00678 END IF 00679 END IF 00680 END IF 00681 ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN 00682 * 00683 * Find normF(A). 00684 * 00685 K = ( N+1 ) / 2 00686 SCALE = ZERO 00687 S = ONE 00688 IF( NOE.EQ.1 ) THEN 00689 * n is odd 00690 IF( IFM.EQ.1 ) THEN 00691 * A is normal 00692 IF( ILU.EQ.0 ) THEN 00693 * A is upper 00694 DO J = 0, K - 3 00695 CALL SLASSQ( K-J-2, A( K+J+1+J*LDA ), 1, SCALE, S ) 00696 * L at A(k,0) 00697 END DO 00698 DO J = 0, K - 1 00699 CALL SLASSQ( K+J-1, A( 0+J*LDA ), 1, SCALE, S ) 00700 * trap U at A(0,0) 00701 END DO 00702 S = S + S 00703 * double s for the off diagonal elements 00704 CALL SLASSQ( K-1, A( K ), LDA+1, SCALE, S ) 00705 * tri L at A(k,0) 00706 CALL SLASSQ( K, A( K-1 ), LDA+1, SCALE, S ) 00707 * tri U at A(k-1,0) 00708 ELSE 00709 * ilu=1 & A is lower 00710 DO J = 0, K - 1 00711 CALL SLASSQ( N-J-1, A( J+1+J*LDA ), 1, SCALE, S ) 00712 * trap L at A(0,0) 00713 END DO 00714 DO J = 0, K - 2 00715 CALL SLASSQ( J, A( 0+( 1+J )*LDA ), 1, SCALE, S ) 00716 * U at A(0,1) 00717 END DO 00718 S = S + S 00719 * double s for the off diagonal elements 00720 CALL SLASSQ( K, A( 0 ), LDA+1, SCALE, S ) 00721 * tri L at A(0,0) 00722 CALL SLASSQ( K-1, A( 0+LDA ), LDA+1, SCALE, S ) 00723 * tri U at A(0,1) 00724 END IF 00725 ELSE 00726 * A is xpose 00727 IF( ILU.EQ.0 ) THEN 00728 * A**T is upper 00729 DO J = 1, K - 2 00730 CALL SLASSQ( J, A( 0+( K+J )*LDA ), 1, SCALE, S ) 00731 * U at A(0,k) 00732 END DO 00733 DO J = 0, K - 2 00734 CALL SLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 00735 * k by k-1 rect. at A(0,0) 00736 END DO 00737 DO J = 0, K - 2 00738 CALL SLASSQ( K-J-1, A( J+1+( J+K-1 )*LDA ), 1, 00739 $ SCALE, S ) 00740 * L at A(0,k-1) 00741 END DO 00742 S = S + S 00743 * double s for the off diagonal elements 00744 CALL SLASSQ( K-1, A( 0+K*LDA ), LDA+1, SCALE, S ) 00745 * tri U at A(0,k) 00746 CALL SLASSQ( K, A( 0+( K-1 )*LDA ), LDA+1, SCALE, S ) 00747 * tri L at A(0,k-1) 00748 ELSE 00749 * A**T is lower 00750 DO J = 1, K - 1 00751 CALL SLASSQ( J, A( 0+J*LDA ), 1, SCALE, S ) 00752 * U at A(0,0) 00753 END DO 00754 DO J = K, N - 1 00755 CALL SLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 00756 * k by k-1 rect. at A(0,k) 00757 END DO 00758 DO J = 0, K - 3 00759 CALL SLASSQ( K-J-2, A( J+2+J*LDA ), 1, SCALE, S ) 00760 * L at A(1,0) 00761 END DO 00762 S = S + S 00763 * double s for the off diagonal elements 00764 CALL SLASSQ( K, A( 0 ), LDA+1, SCALE, S ) 00765 * tri U at A(0,0) 00766 CALL SLASSQ( K-1, A( 1 ), LDA+1, SCALE, S ) 00767 * tri L at A(1,0) 00768 END IF 00769 END IF 00770 ELSE 00771 * n is even 00772 IF( IFM.EQ.1 ) THEN 00773 * A is normal 00774 IF( ILU.EQ.0 ) THEN 00775 * A is upper 00776 DO J = 0, K - 2 00777 CALL SLASSQ( K-J-1, A( K+J+2+J*LDA ), 1, SCALE, S ) 00778 * L at A(k+1,0) 00779 END DO 00780 DO J = 0, K - 1 00781 CALL SLASSQ( K+J, A( 0+J*LDA ), 1, SCALE, S ) 00782 * trap U at A(0,0) 00783 END DO 00784 S = S + S 00785 * double s for the off diagonal elements 00786 CALL SLASSQ( K, A( K+1 ), LDA+1, SCALE, S ) 00787 * tri L at A(k+1,0) 00788 CALL SLASSQ( K, A( K ), LDA+1, SCALE, S ) 00789 * tri U at A(k,0) 00790 ELSE 00791 * ilu=1 & A is lower 00792 DO J = 0, K - 1 00793 CALL SLASSQ( N-J-1, A( J+2+J*LDA ), 1, SCALE, S ) 00794 * trap L at A(1,0) 00795 END DO 00796 DO J = 1, K - 1 00797 CALL SLASSQ( J, A( 0+J*LDA ), 1, SCALE, S ) 00798 * U at A(0,0) 00799 END DO 00800 S = S + S 00801 * double s for the off diagonal elements 00802 CALL SLASSQ( K, A( 1 ), LDA+1, SCALE, S ) 00803 * tri L at A(1,0) 00804 CALL SLASSQ( K, A( 0 ), LDA+1, SCALE, S ) 00805 * tri U at A(0,0) 00806 END IF 00807 ELSE 00808 * A is xpose 00809 IF( ILU.EQ.0 ) THEN 00810 * A**T is upper 00811 DO J = 1, K - 1 00812 CALL SLASSQ( J, A( 0+( K+1+J )*LDA ), 1, SCALE, S ) 00813 * U at A(0,k+1) 00814 END DO 00815 DO J = 0, K - 1 00816 CALL SLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 00817 * k by k rect. at A(0,0) 00818 END DO 00819 DO J = 0, K - 2 00820 CALL SLASSQ( K-J-1, A( J+1+( J+K )*LDA ), 1, SCALE, 00821 $ S ) 00822 * L at A(0,k) 00823 END DO 00824 S = S + S 00825 * double s for the off diagonal elements 00826 CALL SLASSQ( K, A( 0+( K+1 )*LDA ), LDA+1, SCALE, S ) 00827 * tri U at A(0,k+1) 00828 CALL SLASSQ( K, A( 0+K*LDA ), LDA+1, SCALE, S ) 00829 * tri L at A(0,k) 00830 ELSE 00831 * A**T is lower 00832 DO J = 1, K - 1 00833 CALL SLASSQ( J, A( 0+( J+1 )*LDA ), 1, SCALE, S ) 00834 * U at A(0,1) 00835 END DO 00836 DO J = K + 1, N 00837 CALL SLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 00838 * k by k rect. at A(0,k+1) 00839 END DO 00840 DO J = 0, K - 2 00841 CALL SLASSQ( K-J-1, A( J+1+J*LDA ), 1, SCALE, S ) 00842 * L at A(0,0) 00843 END DO 00844 S = S + S 00845 * double s for the off diagonal elements 00846 CALL SLASSQ( K, A( LDA ), LDA+1, SCALE, S ) 00847 * tri L at A(0,1) 00848 CALL SLASSQ( K, A( 0 ), LDA+1, SCALE, S ) 00849 * tri U at A(0,0) 00850 END IF 00851 END IF 00852 END IF 00853 VALUE = SCALE*SQRT( S ) 00854 END IF 00855 * 00856 SLANSF = VALUE 00857 RETURN 00858 * 00859 * End of SLANSF 00860 * 00861 END
1.7.3 
