![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLANSF 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLANSF + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlansf.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlansf.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlansf.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * DOUBLE PRECISION FUNCTION DLANSF( NORM, TRANSR, UPLO, N, A, WORK ) 00022 * 00023 * .. Scalar Arguments .. 00024 * CHARACTER NORM, TRANSR, UPLO 00025 * INTEGER N 00026 * .. 00027 * .. Array Arguments .. 00028 * DOUBLE PRECISION A( 0: * ), WORK( 0: * ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> DLANSF 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 DLANSF 00043 *> \verbatim 00044 *> 00045 *> DLANSF = ( 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 DLANSF 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, DLANSF is 00091 *> set to zero. 00092 *> \endverbatim 00093 *> 00094 *> \param[in] A 00095 *> \verbatim 00096 *> A is DOUBLE PRECISION 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 DOUBLE PRECISION 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 doubleOTHERcomputational 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 DOUBLE PRECISION FUNCTION DLANSF( 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 DOUBLE PRECISION A( 0: * ), WORK( 0: * ) 00223 * .. 00224 * 00225 * ===================================================================== 00226 * 00227 * .. Parameters .. 00228 DOUBLE PRECISION ONE, ZERO 00229 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00230 * .. 00231 * .. Local Scalars .. 00232 INTEGER I, J, IFM, ILU, NOE, N1, K, L, LDA 00233 DOUBLE PRECISION SCALE, S, VALUE, AA 00234 * .. 00235 * .. External Functions .. 00236 LOGICAL LSAME 00237 INTEGER IDAMAX 00238 EXTERNAL LSAME, IDAMAX 00239 * .. 00240 * .. External Subroutines .. 00241 EXTERNAL DLASSQ 00242 * .. 00243 * .. Intrinsic Functions .. 00244 INTRINSIC ABS, MAX, SQRT 00245 * .. 00246 * .. Executable Statements .. 00247 * 00248 IF( N.EQ.0 ) THEN 00249 DLANSF = ZERO 00250 RETURN 00251 ELSE IF( N.EQ.1 ) THEN 00252 DLANSF = ABS( A(0) ) 00253 RETURN 00254 END IF 00255 * 00256 * set noe = 1 if n is odd. if n is even set noe=0 00257 * 00258 NOE = 1 00259 IF( MOD( N, 2 ).EQ.0 ) 00260 $ NOE = 0 00261 * 00262 * set ifm = 0 when form='T or 't' and 1 otherwise 00263 * 00264 IFM = 1 00265 IF( LSAME( TRANSR, 'T' ) ) 00266 $ IFM = 0 00267 * 00268 * set ilu = 0 when uplo='U or 'u' and 1 otherwise 00269 * 00270 ILU = 1 00271 IF( LSAME( UPLO, 'U' ) ) 00272 $ ILU = 0 00273 * 00274 * set lda = (n+1)/2 when ifm = 0 00275 * set lda = n when ifm = 1 and noe = 1 00276 * set lda = n+1 when ifm = 1 and noe = 0 00277 * 00278 IF( IFM.EQ.1 ) THEN 00279 IF( NOE.EQ.1 ) THEN 00280 LDA = N 00281 ELSE 00282 * noe=0 00283 LDA = N + 1 00284 END IF 00285 ELSE 00286 * ifm=0 00287 LDA = ( N+1 ) / 2 00288 END IF 00289 * 00290 IF( LSAME( NORM, 'M' ) ) THEN 00291 * 00292 * Find max(abs(A(i,j))). 00293 * 00294 K = ( N+1 ) / 2 00295 VALUE = ZERO 00296 IF( NOE.EQ.1 ) THEN 00297 * n is odd 00298 IF( IFM.EQ.1 ) THEN 00299 * A is n by k 00300 DO J = 0, K - 1 00301 DO I = 0, N - 1 00302 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00303 END DO 00304 END DO 00305 ELSE 00306 * xpose case; A is k by n 00307 DO J = 0, N - 1 00308 DO I = 0, K - 1 00309 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00310 END DO 00311 END DO 00312 END IF 00313 ELSE 00314 * n is even 00315 IF( IFM.EQ.1 ) THEN 00316 * A is n+1 by k 00317 DO J = 0, K - 1 00318 DO I = 0, N 00319 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00320 END DO 00321 END DO 00322 ELSE 00323 * xpose case; A is k by n+1 00324 DO J = 0, N 00325 DO I = 0, K - 1 00326 VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 00327 END DO 00328 END DO 00329 END IF 00330 END IF 00331 ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR. 00332 $ ( NORM.EQ.'1' ) ) THEN 00333 * 00334 * Find normI(A) ( = norm1(A), since A is symmetric). 00335 * 00336 IF( IFM.EQ.1 ) THEN 00337 K = N / 2 00338 IF( NOE.EQ.1 ) THEN 00339 * n is odd 00340 IF( ILU.EQ.0 ) THEN 00341 DO I = 0, K - 1 00342 WORK( I ) = ZERO 00343 END DO 00344 DO J = 0, K 00345 S = ZERO 00346 DO I = 0, K + J - 1 00347 AA = ABS( A( I+J*LDA ) ) 00348 * -> A(i,j+k) 00349 S = S + AA 00350 WORK( I ) = WORK( I ) + AA 00351 END DO 00352 AA = ABS( A( I+J*LDA ) ) 00353 * -> A(j+k,j+k) 00354 WORK( J+K ) = S + AA 00355 IF( I.EQ.K+K ) 00356 $ GO TO 10 00357 I = I + 1 00358 AA = ABS( A( I+J*LDA ) ) 00359 * -> A(j,j) 00360 WORK( J ) = WORK( J ) + AA 00361 S = ZERO 00362 DO L = J + 1, K - 1 00363 I = I + 1 00364 AA = ABS( A( I+J*LDA ) ) 00365 * -> A(l,j) 00366 S = S + AA 00367 WORK( L ) = WORK( L ) + AA 00368 END DO 00369 WORK( J ) = WORK( J ) + S 00370 END DO 00371 10 CONTINUE 00372 I = IDAMAX( N, WORK, 1 ) 00373 VALUE = WORK( I-1 ) 00374 ELSE 00375 * ilu = 1 00376 K = K + 1 00377 * k=(n+1)/2 for n odd and ilu=1 00378 DO I = K, N - 1 00379 WORK( I ) = ZERO 00380 END DO 00381 DO J = K - 1, 0, -1 00382 S = ZERO 00383 DO I = 0, J - 2 00384 AA = ABS( A( I+J*LDA ) ) 00385 * -> A(j+k,i+k) 00386 S = S + AA 00387 WORK( I+K ) = WORK( I+K ) + AA 00388 END DO 00389 IF( J.GT.0 ) THEN 00390 AA = ABS( A( I+J*LDA ) ) 00391 * -> A(j+k,j+k) 00392 S = S + AA 00393 WORK( I+K ) = WORK( I+K ) + S 00394 * i=j 00395 I = I + 1 00396 END IF 00397 AA = ABS( A( I+J*LDA ) ) 00398 * -> A(j,j) 00399 WORK( J ) = AA 00400 S = ZERO 00401 DO L = J + 1, N - 1 00402 I = I + 1 00403 AA = ABS( A( I+J*LDA ) ) 00404 * -> A(l,j) 00405 S = S + AA 00406 WORK( L ) = WORK( L ) + AA 00407 END DO 00408 WORK( J ) = WORK( J ) + S 00409 END DO 00410 I = IDAMAX( N, WORK, 1 ) 00411 VALUE = WORK( I-1 ) 00412 END IF 00413 ELSE 00414 * n is even 00415 IF( ILU.EQ.0 ) THEN 00416 DO I = 0, K - 1 00417 WORK( I ) = ZERO 00418 END DO 00419 DO J = 0, K - 1 00420 S = ZERO 00421 DO I = 0, K + J - 1 00422 AA = ABS( A( I+J*LDA ) ) 00423 * -> A(i,j+k) 00424 S = S + AA 00425 WORK( I ) = WORK( I ) + AA 00426 END DO 00427 AA = ABS( A( I+J*LDA ) ) 00428 * -> A(j+k,j+k) 00429 WORK( J+K ) = S + AA 00430 I = I + 1 00431 AA = ABS( A( I+J*LDA ) ) 00432 * -> A(j,j) 00433 WORK( J ) = WORK( J ) + AA 00434 S = ZERO 00435 DO L = J + 1, K - 1 00436 I = I + 1 00437 AA = ABS( A( I+J*LDA ) ) 00438 * -> A(l,j) 00439 S = S + AA 00440 WORK( L ) = WORK( L ) + AA 00441 END DO 00442 WORK( J ) = WORK( J ) + S 00443 END DO 00444 I = IDAMAX( N, WORK, 1 ) 00445 VALUE = WORK( I-1 ) 00446 ELSE 00447 * ilu = 1 00448 DO I = K, N - 1 00449 WORK( I ) = ZERO 00450 END DO 00451 DO J = K - 1, 0, -1 00452 S = ZERO 00453 DO I = 0, J - 1 00454 AA = ABS( A( I+J*LDA ) ) 00455 * -> A(j+k,i+k) 00456 S = S + AA 00457 WORK( I+K ) = WORK( I+K ) + AA 00458 END DO 00459 AA = ABS( A( I+J*LDA ) ) 00460 * -> A(j+k,j+k) 00461 S = S + AA 00462 WORK( I+K ) = WORK( I+K ) + S 00463 * i=j 00464 I = I + 1 00465 AA = ABS( A( I+J*LDA ) ) 00466 * -> A(j,j) 00467 WORK( J ) = AA 00468 S = ZERO 00469 DO L = J + 1, N - 1 00470 I = I + 1 00471 AA = ABS( A( I+J*LDA ) ) 00472 * -> A(l,j) 00473 S = S + AA 00474 WORK( L ) = WORK( L ) + AA 00475 END DO 00476 WORK( J ) = WORK( J ) + S 00477 END DO 00478 I = IDAMAX( N, WORK, 1 ) 00479 VALUE = WORK( I-1 ) 00480 END IF 00481 END IF 00482 ELSE 00483 * ifm=0 00484 K = N / 2 00485 IF( NOE.EQ.1 ) THEN 00486 * n is odd 00487 IF( ILU.EQ.0 ) THEN 00488 N1 = K 00489 * n/2 00490 K = K + 1 00491 * k is the row size and lda 00492 DO I = N1, N - 1 00493 WORK( I ) = ZERO 00494 END DO 00495 DO J = 0, N1 - 1 00496 S = ZERO 00497 DO I = 0, K - 1 00498 AA = ABS( A( I+J*LDA ) ) 00499 * A(j,n1+i) 00500 WORK( I+N1 ) = WORK( I+N1 ) + AA 00501 S = S + AA 00502 END DO 00503 WORK( J ) = S 00504 END DO 00505 * j=n1=k-1 is special 00506 S = ABS( A( 0+J*LDA ) ) 00507 * A(k-1,k-1) 00508 DO I = 1, K - 1 00509 AA = ABS( A( I+J*LDA ) ) 00510 * A(k-1,i+n1) 00511 WORK( I+N1 ) = WORK( I+N1 ) + AA 00512 S = S + AA 00513 END DO 00514 WORK( J ) = WORK( J ) + S 00515 DO J = K, N - 1 00516 S = ZERO 00517 DO I = 0, J - K - 1 00518 AA = ABS( A( I+J*LDA ) ) 00519 * A(i,j-k) 00520 WORK( I ) = WORK( I ) + AA 00521 S = S + AA 00522 END DO 00523 * i=j-k 00524 AA = ABS( A( I+J*LDA ) ) 00525 * A(j-k,j-k) 00526 S = S + AA 00527 WORK( J-K ) = WORK( J-K ) + S 00528 I = I + 1 00529 S = ABS( A( I+J*LDA ) ) 00530 * A(j,j) 00531 DO L = J + 1, N - 1 00532 I = I + 1 00533 AA = ABS( A( I+J*LDA ) ) 00534 * A(j,l) 00535 WORK( L ) = WORK( L ) + AA 00536 S = S + AA 00537 END DO 00538 WORK( J ) = WORK( J ) + S 00539 END DO 00540 I = IDAMAX( N, WORK, 1 ) 00541 VALUE = WORK( I-1 ) 00542 ELSE 00543 * ilu=1 00544 K = K + 1 00545 * k=(n+1)/2 for n odd and ilu=1 00546 DO I = K, N - 1 00547 WORK( I ) = ZERO 00548 END DO 00549 DO J = 0, K - 2 00550 * process 00551 S = ZERO 00552 DO I = 0, J - 1 00553 AA = ABS( A( I+J*LDA ) ) 00554 * A(j,i) 00555 WORK( I ) = WORK( I ) + AA 00556 S = S + AA 00557 END DO 00558 AA = ABS( A( I+J*LDA ) ) 00559 * i=j so process of A(j,j) 00560 S = S + AA 00561 WORK( J ) = S 00562 * is initialised here 00563 I = I + 1 00564 * i=j process A(j+k,j+k) 00565 AA = ABS( A( I+J*LDA ) ) 00566 S = AA 00567 DO L = K + J + 1, N - 1 00568 I = I + 1 00569 AA = ABS( A( I+J*LDA ) ) 00570 * A(l,k+j) 00571 S = S + AA 00572 WORK( L ) = WORK( L ) + AA 00573 END DO 00574 WORK( K+J ) = WORK( K+J ) + S 00575 END DO 00576 * j=k-1 is special :process col A(k-1,0:k-1) 00577 S = ZERO 00578 DO I = 0, K - 2 00579 AA = ABS( A( I+J*LDA ) ) 00580 * A(k,i) 00581 WORK( I ) = WORK( I ) + AA 00582 S = S + AA 00583 END DO 00584 * i=k-1 00585 AA = ABS( A( I+J*LDA ) ) 00586 * A(k-1,k-1) 00587 S = S + AA 00588 WORK( I ) = S 00589 * done with col j=k+1 00590 DO J = K, N - 1 00591 * process col j of A = A(j,0:k-1) 00592 S = ZERO 00593 DO I = 0, K - 1 00594 AA = ABS( A( I+J*LDA ) ) 00595 * A(j,i) 00596 WORK( I ) = WORK( I ) + AA 00597 S = S + AA 00598 END DO 00599 WORK( J ) = WORK( J ) + S 00600 END DO 00601 I = IDAMAX( N, WORK, 1 ) 00602 VALUE = WORK( I-1 ) 00603 END IF 00604 ELSE 00605 * n is even 00606 IF( ILU.EQ.0 ) THEN 00607 DO I = K, N - 1 00608 WORK( I ) = ZERO 00609 END DO 00610 DO J = 0, K - 1 00611 S = ZERO 00612 DO I = 0, K - 1 00613 AA = ABS( A( I+J*LDA ) ) 00614 * A(j,i+k) 00615 WORK( I+K ) = WORK( I+K ) + AA 00616 S = S + AA 00617 END DO 00618 WORK( J ) = S 00619 END DO 00620 * j=k 00621 AA = ABS( A( 0+J*LDA ) ) 00622 * A(k,k) 00623 S = AA 00624 DO I = 1, K - 1 00625 AA = ABS( A( I+J*LDA ) ) 00626 * A(k,k+i) 00627 WORK( I+K ) = WORK( I+K ) + AA 00628 S = S + AA 00629 END DO 00630 WORK( J ) = WORK( J ) + S 00631 DO J = K + 1, N - 1 00632 S = ZERO 00633 DO I = 0, J - 2 - K 00634 AA = ABS( A( I+J*LDA ) ) 00635 * A(i,j-k-1) 00636 WORK( I ) = WORK( I ) + AA 00637 S = S + AA 00638 END DO 00639 * i=j-1-k 00640 AA = ABS( A( I+J*LDA ) ) 00641 * A(j-k-1,j-k-1) 00642 S = S + AA 00643 WORK( J-K-1 ) = WORK( J-K-1 ) + S 00644 I = I + 1 00645 AA = ABS( A( I+J*LDA ) ) 00646 * A(j,j) 00647 S = AA 00648 DO L = J + 1, N - 1 00649 I = I + 1 00650 AA = ABS( A( I+J*LDA ) ) 00651 * A(j,l) 00652 WORK( L ) = WORK( L ) + AA 00653 S = S + AA 00654 END DO 00655 WORK( J ) = WORK( J ) + S 00656 END DO 00657 * j=n 00658 S = ZERO 00659 DO I = 0, K - 2 00660 AA = ABS( A( I+J*LDA ) ) 00661 * A(i,k-1) 00662 WORK( I ) = WORK( I ) + AA 00663 S = S + AA 00664 END DO 00665 * i=k-1 00666 AA = ABS( A( I+J*LDA ) ) 00667 * A(k-1,k-1) 00668 S = S + AA 00669 WORK( I ) = WORK( I ) + S 00670 I = IDAMAX( N, WORK, 1 ) 00671 VALUE = WORK( I-1 ) 00672 ELSE 00673 * ilu=1 00674 DO I = K, N - 1 00675 WORK( I ) = ZERO 00676 END DO 00677 * j=0 is special :process col A(k:n-1,k) 00678 S = ABS( A( 0 ) ) 00679 * A(k,k) 00680 DO I = 1, K - 1 00681 AA = ABS( A( I ) ) 00682 * A(k+i,k) 00683 WORK( I+K ) = WORK( I+K ) + AA 00684 S = S + AA 00685 END DO 00686 WORK( K ) = WORK( K ) + S 00687 DO J = 1, K - 1 00688 * process 00689 S = ZERO 00690 DO I = 0, J - 2 00691 AA = ABS( A( I+J*LDA ) ) 00692 * A(j-1,i) 00693 WORK( I ) = WORK( I ) + AA 00694 S = S + AA 00695 END DO 00696 AA = ABS( A( I+J*LDA ) ) 00697 * i=j-1 so process of A(j-1,j-1) 00698 S = S + AA 00699 WORK( J-1 ) = S 00700 * is initialised here 00701 I = I + 1 00702 * i=j process A(j+k,j+k) 00703 AA = ABS( A( I+J*LDA ) ) 00704 S = AA 00705 DO L = K + J + 1, N - 1 00706 I = I + 1 00707 AA = ABS( A( I+J*LDA ) ) 00708 * A(l,k+j) 00709 S = S + AA 00710 WORK( L ) = WORK( L ) + AA 00711 END DO 00712 WORK( K+J ) = WORK( K+J ) + S 00713 END DO 00714 * j=k is special :process col A(k,0:k-1) 00715 S = ZERO 00716 DO I = 0, K - 2 00717 AA = ABS( A( I+J*LDA ) ) 00718 * A(k,i) 00719 WORK( I ) = WORK( I ) + AA 00720 S = S + AA 00721 END DO 00722 * i=k-1 00723 AA = ABS( A( I+J*LDA ) ) 00724 * A(k-1,k-1) 00725 S = S + AA 00726 WORK( I ) = S 00727 * done with col j=k+1 00728 DO J = K + 1, N 00729 * process col j-1 of A = A(j-1,0:k-1) 00730 S = ZERO 00731 DO I = 0, K - 1 00732 AA = ABS( A( I+J*LDA ) ) 00733 * A(j-1,i) 00734 WORK( I ) = WORK( I ) + AA 00735 S = S + AA 00736 END DO 00737 WORK( J-1 ) = WORK( J-1 ) + S 00738 END DO 00739 I = IDAMAX( N, WORK, 1 ) 00740 VALUE = WORK( I-1 ) 00741 END IF 00742 END IF 00743 END IF 00744 ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN 00745 * 00746 * Find normF(A). 00747 * 00748 K = ( N+1 ) / 2 00749 SCALE = ZERO 00750 S = ONE 00751 IF( NOE.EQ.1 ) THEN 00752 * n is odd 00753 IF( IFM.EQ.1 ) THEN 00754 * A is normal 00755 IF( ILU.EQ.0 ) THEN 00756 * A is upper 00757 DO J = 0, K - 3 00758 CALL DLASSQ( K-J-2, A( K+J+1+J*LDA ), 1, SCALE, S ) 00759 * L at A(k,0) 00760 END DO 00761 DO J = 0, K - 1 00762 CALL DLASSQ( K+J-1, A( 0+J*LDA ), 1, SCALE, S ) 00763 * trap U at A(0,0) 00764 END DO 00765 S = S + S 00766 * double s for the off diagonal elements 00767 CALL DLASSQ( K-1, A( K ), LDA+1, SCALE, S ) 00768 * tri L at A(k,0) 00769 CALL DLASSQ( K, A( K-1 ), LDA+1, SCALE, S ) 00770 * tri U at A(k-1,0) 00771 ELSE 00772 * ilu=1 & A is lower 00773 DO J = 0, K - 1 00774 CALL DLASSQ( N-J-1, A( J+1+J*LDA ), 1, SCALE, S ) 00775 * trap L at A(0,0) 00776 END DO 00777 DO J = 0, K - 2 00778 CALL DLASSQ( J, A( 0+( 1+J )*LDA ), 1, SCALE, S ) 00779 * U at A(0,1) 00780 END DO 00781 S = S + S 00782 * double s for the off diagonal elements 00783 CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S ) 00784 * tri L at A(0,0) 00785 CALL DLASSQ( K-1, A( 0+LDA ), LDA+1, SCALE, S ) 00786 * tri U at A(0,1) 00787 END IF 00788 ELSE 00789 * A is xpose 00790 IF( ILU.EQ.0 ) THEN 00791 * A**T is upper 00792 DO J = 1, K - 2 00793 CALL DLASSQ( J, A( 0+( K+J )*LDA ), 1, SCALE, S ) 00794 * U at A(0,k) 00795 END DO 00796 DO J = 0, K - 2 00797 CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 00798 * k by k-1 rect. at A(0,0) 00799 END DO 00800 DO J = 0, K - 2 00801 CALL DLASSQ( K-J-1, A( J+1+( J+K-1 )*LDA ), 1, 00802 $ SCALE, S ) 00803 * L at A(0,k-1) 00804 END DO 00805 S = S + S 00806 * double s for the off diagonal elements 00807 CALL DLASSQ( K-1, A( 0+K*LDA ), LDA+1, SCALE, S ) 00808 * tri U at A(0,k) 00809 CALL DLASSQ( K, A( 0+( K-1 )*LDA ), LDA+1, SCALE, S ) 00810 * tri L at A(0,k-1) 00811 ELSE 00812 * A**T is lower 00813 DO J = 1, K - 1 00814 CALL DLASSQ( J, A( 0+J*LDA ), 1, SCALE, S ) 00815 * U at A(0,0) 00816 END DO 00817 DO J = K, N - 1 00818 CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 00819 * k by k-1 rect. at A(0,k) 00820 END DO 00821 DO J = 0, K - 3 00822 CALL DLASSQ( K-J-2, A( J+2+J*LDA ), 1, SCALE, S ) 00823 * L at A(1,0) 00824 END DO 00825 S = S + S 00826 * double s for the off diagonal elements 00827 CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S ) 00828 * tri U at A(0,0) 00829 CALL DLASSQ( K-1, A( 1 ), LDA+1, SCALE, S ) 00830 * tri L at A(1,0) 00831 END IF 00832 END IF 00833 ELSE 00834 * n is even 00835 IF( IFM.EQ.1 ) THEN 00836 * A is normal 00837 IF( ILU.EQ.0 ) THEN 00838 * A is upper 00839 DO J = 0, K - 2 00840 CALL DLASSQ( K-J-1, A( K+J+2+J*LDA ), 1, SCALE, S ) 00841 * L at A(k+1,0) 00842 END DO 00843 DO J = 0, K - 1 00844 CALL DLASSQ( K+J, A( 0+J*LDA ), 1, SCALE, S ) 00845 * trap U at A(0,0) 00846 END DO 00847 S = S + S 00848 * double s for the off diagonal elements 00849 CALL DLASSQ( K, A( K+1 ), LDA+1, SCALE, S ) 00850 * tri L at A(k+1,0) 00851 CALL DLASSQ( K, A( K ), LDA+1, SCALE, S ) 00852 * tri U at A(k,0) 00853 ELSE 00854 * ilu=1 & A is lower 00855 DO J = 0, K - 1 00856 CALL DLASSQ( N-J-1, A( J+2+J*LDA ), 1, SCALE, S ) 00857 * trap L at A(1,0) 00858 END DO 00859 DO J = 1, K - 1 00860 CALL DLASSQ( J, A( 0+J*LDA ), 1, SCALE, S ) 00861 * U at A(0,0) 00862 END DO 00863 S = S + S 00864 * double s for the off diagonal elements 00865 CALL DLASSQ( K, A( 1 ), LDA+1, SCALE, S ) 00866 * tri L at A(1,0) 00867 CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S ) 00868 * tri U at A(0,0) 00869 END IF 00870 ELSE 00871 * A is xpose 00872 IF( ILU.EQ.0 ) THEN 00873 * A**T is upper 00874 DO J = 1, K - 1 00875 CALL DLASSQ( J, A( 0+( K+1+J )*LDA ), 1, SCALE, S ) 00876 * U at A(0,k+1) 00877 END DO 00878 DO J = 0, K - 1 00879 CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 00880 * k by k rect. at A(0,0) 00881 END DO 00882 DO J = 0, K - 2 00883 CALL DLASSQ( K-J-1, A( J+1+( J+K )*LDA ), 1, SCALE, 00884 $ S ) 00885 * L at A(0,k) 00886 END DO 00887 S = S + S 00888 * double s for the off diagonal elements 00889 CALL DLASSQ( K, A( 0+( K+1 )*LDA ), LDA+1, SCALE, S ) 00890 * tri U at A(0,k+1) 00891 CALL DLASSQ( K, A( 0+K*LDA ), LDA+1, SCALE, S ) 00892 * tri L at A(0,k) 00893 ELSE 00894 * A**T is lower 00895 DO J = 1, K - 1 00896 CALL DLASSQ( J, A( 0+( J+1 )*LDA ), 1, SCALE, S ) 00897 * U at A(0,1) 00898 END DO 00899 DO J = K + 1, N 00900 CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 00901 * k by k rect. at A(0,k+1) 00902 END DO 00903 DO J = 0, K - 2 00904 CALL DLASSQ( K-J-1, A( J+1+J*LDA ), 1, SCALE, S ) 00905 * L at A(0,0) 00906 END DO 00907 S = S + S 00908 * double s for the off diagonal elements 00909 CALL DLASSQ( K, A( LDA ), LDA+1, SCALE, S ) 00910 * tri L at A(0,1) 00911 CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S ) 00912 * tri U at A(0,0) 00913 END IF 00914 END IF 00915 END IF 00916 VALUE = SCALE*SQRT( S ) 00917 END IF 00918 * 00919 DLANSF = VALUE 00920 RETURN 00921 * 00922 * End of DLANSF 00923 * 00924 END