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