![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: DOUBLE PRECISION FUNCTION DLANSF( NORM, TRANSR, UPLO, N, A, WORK ) 2: * 3: * -- LAPACK routine (version 3.3.0) -- 4: * 5: * -- Contributed by Fred Gustavson of the IBM Watson Research Center -- 6: * November 2010 7: * 8: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 9: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 10: * 11: * .. Scalar Arguments .. 12: CHARACTER NORM, TRANSR, UPLO 13: INTEGER N 14: * .. 15: * .. Array Arguments .. 16: DOUBLE PRECISION A( 0: * ), WORK( 0: * ) 17: * .. 18: * 19: * Purpose 20: * ======= 21: * 22: * DLANSF returns the value of the one norm, or the Frobenius norm, or 23: * the infinity norm, or the element of largest absolute value of a 24: * real symmetric matrix A in RFP format. 25: * 26: * Description 27: * =========== 28: * 29: * DLANSF returns the value 30: * 31: * DLANSF = ( max(abs(A(i,j))), NORM = 'M' or 'm' 32: * ( 33: * ( norm1(A), NORM = '1', 'O' or 'o' 34: * ( 35: * ( normI(A), NORM = 'I' or 'i' 36: * ( 37: * ( normF(A), NORM = 'F', 'f', 'E' or 'e' 38: * 39: * where norm1 denotes the one norm of a matrix (maximum column sum), 40: * normI denotes the infinity norm of a matrix (maximum row sum) and 41: * normF denotes the Frobenius norm of a matrix (square root of sum of 42: * squares). Note that max(abs(A(i,j))) is not a matrix norm. 43: * 44: * Arguments 45: * ========= 46: * 47: * NORM (input) CHARACTER*1 48: * Specifies the value to be returned in DLANSF as described 49: * above. 50: * 51: * TRANSR (input) CHARACTER*1 52: * Specifies whether the RFP format of A is normal or 53: * transposed format. 54: * = 'N': RFP format is Normal; 55: * = 'T': RFP format is Transpose. 56: * 57: * UPLO (input) CHARACTER*1 58: * On entry, UPLO specifies whether the RFP matrix A came from 59: * an upper or lower triangular matrix as follows: 60: * = 'U': RFP A came from an upper triangular matrix; 61: * = 'L': RFP A came from a lower triangular matrix. 62: * 63: * N (input) INTEGER 64: * The order of the matrix A. N >= 0. When N = 0, DLANSF is 65: * set to zero. 66: * 67: * A (input) DOUBLE PRECISION array, dimension ( N*(N+1)/2 ); 68: * On entry, the upper (if UPLO = 'U') or lower (if UPLO = 'L') 69: * part of the symmetric matrix A stored in RFP format. See the 70: * "Notes" below for more details. 71: * Unchanged on exit. 72: * 73: * WORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)), 74: * where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, 75: * WORK is not referenced. 76: * 77: * Further Details 78: * =============== 79: * 80: * We first consider Rectangular Full Packed (RFP) Format when N is 81: * even. We give an example where N = 6. 82: * 83: * AP is Upper AP is Lower 84: * 85: * 00 01 02 03 04 05 00 86: * 11 12 13 14 15 10 11 87: * 22 23 24 25 20 21 22 88: * 33 34 35 30 31 32 33 89: * 44 45 40 41 42 43 44 90: * 55 50 51 52 53 54 55 91: * 92: * 93: * Let TRANSR = 'N'. RFP holds AP as follows: 94: * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last 95: * three columns of AP upper. The lower triangle A(4:6,0:2) consists of 96: * the transpose of the first three columns of AP upper. 97: * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first 98: * three columns of AP lower. The upper triangle A(0:2,0:2) consists of 99: * the transpose of the last three columns of AP lower. 100: * This covers the case N even and TRANSR = 'N'. 101: * 102: * RFP A RFP A 103: * 104: * 03 04 05 33 43 53 105: * 13 14 15 00 44 54 106: * 23 24 25 10 11 55 107: * 33 34 35 20 21 22 108: * 00 44 45 30 31 32 109: * 01 11 55 40 41 42 110: * 02 12 22 50 51 52 111: * 112: * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the 113: * transpose of RFP A above. One therefore gets: 114: * 115: * 116: * RFP A RFP A 117: * 118: * 03 13 23 33 00 01 02 33 00 10 20 30 40 50 119: * 04 14 24 34 44 11 12 43 44 11 21 31 41 51 120: * 05 15 25 35 45 55 22 53 54 55 22 32 42 52 121: * 122: * 123: * We then consider Rectangular Full Packed (RFP) Format when N is 124: * odd. We give an example where N = 5. 125: * 126: * AP is Upper AP is Lower 127: * 128: * 00 01 02 03 04 00 129: * 11 12 13 14 10 11 130: * 22 23 24 20 21 22 131: * 33 34 30 31 32 33 132: * 44 40 41 42 43 44 133: * 134: * 135: * Let TRANSR = 'N'. RFP holds AP as follows: 136: * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last 137: * three columns of AP upper. The lower triangle A(3:4,0:1) consists of 138: * the transpose of the first two columns of AP upper. 139: * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first 140: * three columns of AP lower. The upper triangle A(0:1,1:2) consists of 141: * the transpose of the last two columns of AP lower. 142: * This covers the case N odd and TRANSR = 'N'. 143: * 144: * RFP A RFP A 145: * 146: * 02 03 04 00 33 43 147: * 12 13 14 10 11 44 148: * 22 23 24 20 21 22 149: * 00 33 34 30 31 32 150: * 01 11 44 40 41 42 151: * 152: * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the 153: * transpose of RFP A above. One therefore gets: 154: * 155: * RFP A RFP A 156: * 157: * 02 12 22 00 01 00 10 20 30 40 50 158: * 03 13 23 33 11 33 11 21 31 41 51 159: * 04 14 24 34 44 43 44 22 32 42 52 160: * 161: * Reference 162: * ========= 163: * 164: * ===================================================================== 165: * 166: * .. Parameters .. 167: DOUBLE PRECISION ONE, ZERO 168: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 169: * .. 170: * .. Local Scalars .. 171: INTEGER I, J, IFM, ILU, NOE, N1, K, L, LDA 172: DOUBLE PRECISION SCALE, S, VALUE, AA 173: * .. 174: * .. External Functions .. 175: LOGICAL LSAME 176: INTEGER IDAMAX 177: EXTERNAL LSAME, IDAMAX 178: * .. 179: * .. External Subroutines .. 180: EXTERNAL DLASSQ 181: * .. 182: * .. Intrinsic Functions .. 183: INTRINSIC ABS, MAX, SQRT 184: * .. 185: * .. Executable Statements .. 186: * 187: IF( N.EQ.0 ) THEN 188: DLANSF = ZERO 189: RETURN 190: END IF 191: * 192: * set noe = 1 if n is odd. if n is even set noe=0 193: * 194: NOE = 1 195: IF( MOD( N, 2 ).EQ.0 ) 196: + NOE = 0 197: * 198: * set ifm = 0 when form='T or 't' and 1 otherwise 199: * 200: IFM = 1 201: IF( LSAME( TRANSR, 'T' ) ) 202: + IFM = 0 203: * 204: * set ilu = 0 when uplo='U or 'u' and 1 otherwise 205: * 206: ILU = 1 207: IF( LSAME( UPLO, 'U' ) ) 208: + ILU = 0 209: * 210: * set lda = (n+1)/2 when ifm = 0 211: * set lda = n when ifm = 1 and noe = 1 212: * set lda = n+1 when ifm = 1 and noe = 0 213: * 214: IF( IFM.EQ.1 ) THEN 215: IF( NOE.EQ.1 ) THEN 216: LDA = N 217: ELSE 218: * noe=0 219: LDA = N + 1 220: END IF 221: ELSE 222: * ifm=0 223: LDA = ( N+1 ) / 2 224: END IF 225: * 226: IF( LSAME( NORM, 'M' ) ) THEN 227: * 228: * Find max(abs(A(i,j))). 229: * 230: K = ( N+1 ) / 2 231: VALUE = ZERO 232: IF( NOE.EQ.1 ) THEN 233: * n is odd 234: IF( IFM.EQ.1 ) THEN 235: * A is n by k 236: DO J = 0, K - 1 237: DO I = 0, N - 1 238: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 239: END DO 240: END DO 241: ELSE 242: * xpose case; A is k by n 243: DO J = 0, N - 1 244: DO I = 0, K - 1 245: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 246: END DO 247: END DO 248: END IF 249: ELSE 250: * n is even 251: IF( IFM.EQ.1 ) THEN 252: * A is n+1 by k 253: DO J = 0, K - 1 254: DO I = 0, N 255: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 256: END DO 257: END DO 258: ELSE 259: * xpose case; A is k by n+1 260: DO J = 0, N 261: DO I = 0, K - 1 262: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 263: END DO 264: END DO 265: END IF 266: END IF 267: ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR. 268: + ( NORM.EQ.'1' ) ) THEN 269: * 270: * Find normI(A) ( = norm1(A), since A is symmetric). 271: * 272: IF( IFM.EQ.1 ) THEN 273: K = N / 2 274: IF( NOE.EQ.1 ) THEN 275: * n is odd 276: IF( ILU.EQ.0 ) THEN 277: DO I = 0, K - 1 278: WORK( I ) = ZERO 279: END DO 280: DO J = 0, K 281: S = ZERO 282: DO I = 0, K + J - 1 283: AA = ABS( A( I+J*LDA ) ) 284: * -> A(i,j+k) 285: S = S + AA 286: WORK( I ) = WORK( I ) + AA 287: END DO 288: AA = ABS( A( I+J*LDA ) ) 289: * -> A(j+k,j+k) 290: WORK( J+K ) = S + AA 291: IF( I.EQ.K+K ) 292: + GO TO 10 293: I = I + 1 294: AA = ABS( A( I+J*LDA ) ) 295: * -> A(j,j) 296: WORK( J ) = WORK( J ) + AA 297: S = ZERO 298: DO L = J + 1, K - 1 299: I = I + 1 300: AA = ABS( A( I+J*LDA ) ) 301: * -> A(l,j) 302: S = S + AA 303: WORK( L ) = WORK( L ) + AA 304: END DO 305: WORK( J ) = WORK( J ) + S 306: END DO 307: 10 CONTINUE 308: I = IDAMAX( N, WORK, 1 ) 309: VALUE = WORK( I-1 ) 310: ELSE 311: * ilu = 1 312: K = K + 1 313: * k=(n+1)/2 for n odd and ilu=1 314: DO I = K, N - 1 315: WORK( I ) = ZERO 316: END DO 317: DO J = K - 1, 0, -1 318: S = ZERO 319: DO I = 0, J - 2 320: AA = ABS( A( I+J*LDA ) ) 321: * -> A(j+k,i+k) 322: S = S + AA 323: WORK( I+K ) = WORK( I+K ) + AA 324: END DO 325: IF( J.GT.0 ) THEN 326: AA = ABS( A( I+J*LDA ) ) 327: * -> A(j+k,j+k) 328: S = S + AA 329: WORK( I+K ) = WORK( I+K ) + S 330: * i=j 331: I = I + 1 332: END IF 333: AA = ABS( A( I+J*LDA ) ) 334: * -> A(j,j) 335: WORK( J ) = AA 336: S = ZERO 337: DO L = J + 1, N - 1 338: I = I + 1 339: AA = ABS( A( I+J*LDA ) ) 340: * -> A(l,j) 341: S = S + AA 342: WORK( L ) = WORK( L ) + AA 343: END DO 344: WORK( J ) = WORK( J ) + S 345: END DO 346: I = IDAMAX( N, WORK, 1 ) 347: VALUE = WORK( I-1 ) 348: END IF 349: ELSE 350: * n is even 351: IF( ILU.EQ.0 ) THEN 352: DO I = 0, K - 1 353: WORK( I ) = ZERO 354: END DO 355: DO J = 0, K - 1 356: S = ZERO 357: DO I = 0, K + J - 1 358: AA = ABS( A( I+J*LDA ) ) 359: * -> A(i,j+k) 360: S = S + AA 361: WORK( I ) = WORK( I ) + AA 362: END DO 363: AA = ABS( A( I+J*LDA ) ) 364: * -> A(j+k,j+k) 365: WORK( J+K ) = S + AA 366: I = I + 1 367: AA = ABS( A( I+J*LDA ) ) 368: * -> A(j,j) 369: WORK( J ) = WORK( J ) + AA 370: S = ZERO 371: DO L = J + 1, K - 1 372: I = I + 1 373: AA = ABS( A( I+J*LDA ) ) 374: * -> A(l,j) 375: S = S + AA 376: WORK( L ) = WORK( L ) + AA 377: END DO 378: WORK( J ) = WORK( J ) + S 379: END DO 380: I = IDAMAX( N, WORK, 1 ) 381: VALUE = WORK( I-1 ) 382: ELSE 383: * ilu = 1 384: DO I = K, N - 1 385: WORK( I ) = ZERO 386: END DO 387: DO J = K - 1, 0, -1 388: S = ZERO 389: DO I = 0, J - 1 390: AA = ABS( A( I+J*LDA ) ) 391: * -> A(j+k,i+k) 392: S = S + AA 393: WORK( I+K ) = WORK( I+K ) + AA 394: END DO 395: AA = ABS( A( I+J*LDA ) ) 396: * -> A(j+k,j+k) 397: S = S + AA 398: WORK( I+K ) = WORK( I+K ) + S 399: * i=j 400: I = I + 1 401: AA = ABS( A( I+J*LDA ) ) 402: * -> A(j,j) 403: WORK( J ) = AA 404: S = ZERO 405: DO L = J + 1, N - 1 406: I = I + 1 407: AA = ABS( A( I+J*LDA ) ) 408: * -> A(l,j) 409: S = S + AA 410: WORK( L ) = WORK( L ) + AA 411: END DO 412: WORK( J ) = WORK( J ) + S 413: END DO 414: I = IDAMAX( N, WORK, 1 ) 415: VALUE = WORK( I-1 ) 416: END IF 417: END IF 418: ELSE 419: * ifm=0 420: K = N / 2 421: IF( NOE.EQ.1 ) THEN 422: * n is odd 423: IF( ILU.EQ.0 ) THEN 424: N1 = K 425: * n/2 426: K = K + 1 427: * k is the row size and lda 428: DO I = N1, N - 1 429: WORK( I ) = ZERO 430: END DO 431: DO J = 0, N1 - 1 432: S = ZERO 433: DO I = 0, K - 1 434: AA = ABS( A( I+J*LDA ) ) 435: * A(j,n1+i) 436: WORK( I+N1 ) = WORK( I+N1 ) + AA 437: S = S + AA 438: END DO 439: WORK( J ) = S 440: END DO 441: * j=n1=k-1 is special 442: S = ABS( A( 0+J*LDA ) ) 443: * A(k-1,k-1) 444: DO I = 1, K - 1 445: AA = ABS( A( I+J*LDA ) ) 446: * A(k-1,i+n1) 447: WORK( I+N1 ) = WORK( I+N1 ) + AA 448: S = S + AA 449: END DO 450: WORK( J ) = WORK( J ) + S 451: DO J = K, N - 1 452: S = ZERO 453: DO I = 0, J - K - 1 454: AA = ABS( A( I+J*LDA ) ) 455: * A(i,j-k) 456: WORK( I ) = WORK( I ) + AA 457: S = S + AA 458: END DO 459: * i=j-k 460: AA = ABS( A( I+J*LDA ) ) 461: * A(j-k,j-k) 462: S = S + AA 463: WORK( J-K ) = WORK( J-K ) + S 464: I = I + 1 465: S = ABS( A( I+J*LDA ) ) 466: * A(j,j) 467: DO L = J + 1, N - 1 468: I = I + 1 469: AA = ABS( A( I+J*LDA ) ) 470: * A(j,l) 471: WORK( L ) = WORK( L ) + AA 472: S = S + AA 473: END DO 474: WORK( J ) = WORK( J ) + S 475: END DO 476: I = IDAMAX( N, WORK, 1 ) 477: VALUE = WORK( I-1 ) 478: ELSE 479: * ilu=1 480: K = K + 1 481: * k=(n+1)/2 for n odd and ilu=1 482: DO I = K, N - 1 483: WORK( I ) = ZERO 484: END DO 485: DO J = 0, K - 2 486: * process 487: S = ZERO 488: DO I = 0, J - 1 489: AA = ABS( A( I+J*LDA ) ) 490: * A(j,i) 491: WORK( I ) = WORK( I ) + AA 492: S = S + AA 493: END DO 494: AA = ABS( A( I+J*LDA ) ) 495: * i=j so process of A(j,j) 496: S = S + AA 497: WORK( J ) = S 498: * is initialised here 499: I = I + 1 500: * i=j process A(j+k,j+k) 501: AA = ABS( A( I+J*LDA ) ) 502: S = AA 503: DO L = K + J + 1, N - 1 504: I = I + 1 505: AA = ABS( A( I+J*LDA ) ) 506: * A(l,k+j) 507: S = S + AA 508: WORK( L ) = WORK( L ) + AA 509: END DO 510: WORK( K+J ) = WORK( K+J ) + S 511: END DO 512: * j=k-1 is special :process col A(k-1,0:k-1) 513: S = ZERO 514: DO I = 0, K - 2 515: AA = ABS( A( I+J*LDA ) ) 516: * A(k,i) 517: WORK( I ) = WORK( I ) + AA 518: S = S + AA 519: END DO 520: * i=k-1 521: AA = ABS( A( I+J*LDA ) ) 522: * A(k-1,k-1) 523: S = S + AA 524: WORK( I ) = S 525: * done with col j=k+1 526: DO J = K, N - 1 527: * process col j of A = A(j,0:k-1) 528: S = ZERO 529: DO I = 0, K - 1 530: AA = ABS( A( I+J*LDA ) ) 531: * A(j,i) 532: WORK( I ) = WORK( I ) + AA 533: S = S + AA 534: END DO 535: WORK( J ) = WORK( J ) + S 536: END DO 537: I = IDAMAX( N, WORK, 1 ) 538: VALUE = WORK( I-1 ) 539: END IF 540: ELSE 541: * n is even 542: IF( ILU.EQ.0 ) THEN 543: DO I = K, N - 1 544: WORK( I ) = ZERO 545: END DO 546: DO J = 0, K - 1 547: S = ZERO 548: DO I = 0, K - 1 549: AA = ABS( A( I+J*LDA ) ) 550: * A(j,i+k) 551: WORK( I+K ) = WORK( I+K ) + AA 552: S = S + AA 553: END DO 554: WORK( J ) = S 555: END DO 556: * j=k 557: AA = ABS( A( 0+J*LDA ) ) 558: * A(k,k) 559: S = AA 560: DO I = 1, K - 1 561: AA = ABS( A( I+J*LDA ) ) 562: * A(k,k+i) 563: WORK( I+K ) = WORK( I+K ) + AA 564: S = S + AA 565: END DO 566: WORK( J ) = WORK( J ) + S 567: DO J = K + 1, N - 1 568: S = ZERO 569: DO I = 0, J - 2 - K 570: AA = ABS( A( I+J*LDA ) ) 571: * A(i,j-k-1) 572: WORK( I ) = WORK( I ) + AA 573: S = S + AA 574: END DO 575: * i=j-1-k 576: AA = ABS( A( I+J*LDA ) ) 577: * A(j-k-1,j-k-1) 578: S = S + AA 579: WORK( J-K-1 ) = WORK( J-K-1 ) + S 580: I = I + 1 581: AA = ABS( A( I+J*LDA ) ) 582: * A(j,j) 583: S = AA 584: DO L = J + 1, N - 1 585: I = I + 1 586: AA = ABS( A( I+J*LDA ) ) 587: * A(j,l) 588: WORK( L ) = WORK( L ) + AA 589: S = S + AA 590: END DO 591: WORK( J ) = WORK( J ) + S 592: END DO 593: * j=n 594: S = ZERO 595: DO I = 0, K - 2 596: AA = ABS( A( I+J*LDA ) ) 597: * A(i,k-1) 598: WORK( I ) = WORK( I ) + AA 599: S = S + AA 600: END DO 601: * i=k-1 602: AA = ABS( A( I+J*LDA ) ) 603: * A(k-1,k-1) 604: S = S + AA 605: WORK( I ) = WORK( I ) + S 606: I = IDAMAX( N, WORK, 1 ) 607: VALUE = WORK( I-1 ) 608: ELSE 609: * ilu=1 610: DO I = K, N - 1 611: WORK( I ) = ZERO 612: END DO 613: * j=0 is special :process col A(k:n-1,k) 614: S = ABS( A( 0 ) ) 615: * A(k,k) 616: DO I = 1, K - 1 617: AA = ABS( A( I ) ) 618: * A(k+i,k) 619: WORK( I+K ) = WORK( I+K ) + AA 620: S = S + AA 621: END DO 622: WORK( K ) = WORK( K ) + S 623: DO J = 1, K - 1 624: * process 625: S = ZERO 626: DO I = 0, J - 2 627: AA = ABS( A( I+J*LDA ) ) 628: * A(j-1,i) 629: WORK( I ) = WORK( I ) + AA 630: S = S + AA 631: END DO 632: AA = ABS( A( I+J*LDA ) ) 633: * i=j-1 so process of A(j-1,j-1) 634: S = S + AA 635: WORK( J-1 ) = S 636: * is initialised here 637: I = I + 1 638: * i=j process A(j+k,j+k) 639: AA = ABS( A( I+J*LDA ) ) 640: S = AA 641: DO L = K + J + 1, N - 1 642: I = I + 1 643: AA = ABS( A( I+J*LDA ) ) 644: * A(l,k+j) 645: S = S + AA 646: WORK( L ) = WORK( L ) + AA 647: END DO 648: WORK( K+J ) = WORK( K+J ) + S 649: END DO 650: * j=k is special :process col A(k,0:k-1) 651: S = ZERO 652: DO I = 0, K - 2 653: AA = ABS( A( I+J*LDA ) ) 654: * A(k,i) 655: WORK( I ) = WORK( I ) + AA 656: S = S + AA 657: END DO 658: * i=k-1 659: AA = ABS( A( I+J*LDA ) ) 660: * A(k-1,k-1) 661: S = S + AA 662: WORK( I ) = S 663: * done with col j=k+1 664: DO J = K + 1, N 665: * process col j-1 of A = A(j-1,0:k-1) 666: S = ZERO 667: DO I = 0, K - 1 668: AA = ABS( A( I+J*LDA ) ) 669: * A(j-1,i) 670: WORK( I ) = WORK( I ) + AA 671: S = S + AA 672: END DO 673: WORK( J-1 ) = WORK( J-1 ) + S 674: END DO 675: I = IDAMAX( N, WORK, 1 ) 676: VALUE = WORK( I-1 ) 677: END IF 678: END IF 679: END IF 680: ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN 681: * 682: * Find normF(A). 683: * 684: K = ( N+1 ) / 2 685: SCALE = ZERO 686: S = ONE 687: IF( NOE.EQ.1 ) THEN 688: * n is odd 689: IF( IFM.EQ.1 ) THEN 690: * A is normal 691: IF( ILU.EQ.0 ) THEN 692: * A is upper 693: DO J = 0, K - 3 694: CALL DLASSQ( K-J-2, A( K+J+1+J*LDA ), 1, SCALE, S ) 695: * L at A(k,0) 696: END DO 697: DO J = 0, K - 1 698: CALL DLASSQ( K+J-1, A( 0+J*LDA ), 1, SCALE, S ) 699: * trap U at A(0,0) 700: END DO 701: S = S + S 702: * double s for the off diagonal elements 703: CALL DLASSQ( K-1, A( K ), LDA+1, SCALE, S ) 704: * tri L at A(k,0) 705: CALL DLASSQ( K, A( K-1 ), LDA+1, SCALE, S ) 706: * tri U at A(k-1,0) 707: ELSE 708: * ilu=1 & A is lower 709: DO J = 0, K - 1 710: CALL DLASSQ( N-J-1, A( J+1+J*LDA ), 1, SCALE, S ) 711: * trap L at A(0,0) 712: END DO 713: DO J = 0, K - 2 714: CALL DLASSQ( J, A( 0+( 1+J )*LDA ), 1, SCALE, S ) 715: * U at A(0,1) 716: END DO 717: S = S + S 718: * double s for the off diagonal elements 719: CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S ) 720: * tri L at A(0,0) 721: CALL DLASSQ( K-1, A( 0+LDA ), LDA+1, SCALE, S ) 722: * tri U at A(0,1) 723: END IF 724: ELSE 725: * A is xpose 726: IF( ILU.EQ.0 ) THEN 727: * A' is upper 728: DO J = 1, K - 2 729: CALL DLASSQ( J, A( 0+( K+J )*LDA ), 1, SCALE, S ) 730: * U at A(0,k) 731: END DO 732: DO J = 0, K - 2 733: CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 734: * k by k-1 rect. at A(0,0) 735: END DO 736: DO J = 0, K - 2 737: CALL DLASSQ( K-J-1, A( J+1+( J+K-1 )*LDA ), 1, 738: + SCALE, S ) 739: * L at A(0,k-1) 740: END DO 741: S = S + S 742: * double s for the off diagonal elements 743: CALL DLASSQ( K-1, A( 0+K*LDA ), LDA+1, SCALE, S ) 744: * tri U at A(0,k) 745: CALL DLASSQ( K, A( 0+( K-1 )*LDA ), LDA+1, SCALE, S ) 746: * tri L at A(0,k-1) 747: ELSE 748: * A' is lower 749: DO J = 1, K - 1 750: CALL DLASSQ( J, A( 0+J*LDA ), 1, SCALE, S ) 751: * U at A(0,0) 752: END DO 753: DO J = K, N - 1 754: CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 755: * k by k-1 rect. at A(0,k) 756: END DO 757: DO J = 0, K - 3 758: CALL DLASSQ( K-J-2, A( J+2+J*LDA ), 1, SCALE, S ) 759: * L at A(1,0) 760: END DO 761: S = S + S 762: * double s for the off diagonal elements 763: CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S ) 764: * tri U at A(0,0) 765: CALL DLASSQ( K-1, A( 1 ), LDA+1, SCALE, S ) 766: * tri L at A(1,0) 767: END IF 768: END IF 769: ELSE 770: * n is even 771: IF( IFM.EQ.1 ) THEN 772: * A is normal 773: IF( ILU.EQ.0 ) THEN 774: * A is upper 775: DO J = 0, K - 2 776: CALL DLASSQ( K-J-1, A( K+J+2+J*LDA ), 1, SCALE, S ) 777: * L at A(k+1,0) 778: END DO 779: DO J = 0, K - 1 780: CALL DLASSQ( K+J, A( 0+J*LDA ), 1, SCALE, S ) 781: * trap U at A(0,0) 782: END DO 783: S = S + S 784: * double s for the off diagonal elements 785: CALL DLASSQ( K, A( K+1 ), LDA+1, SCALE, S ) 786: * tri L at A(k+1,0) 787: CALL DLASSQ( K, A( K ), LDA+1, SCALE, S ) 788: * tri U at A(k,0) 789: ELSE 790: * ilu=1 & A is lower 791: DO J = 0, K - 1 792: CALL DLASSQ( N-J-1, A( J+2+J*LDA ), 1, SCALE, S ) 793: * trap L at A(1,0) 794: END DO 795: DO J = 1, K - 1 796: CALL DLASSQ( J, A( 0+J*LDA ), 1, SCALE, S ) 797: * U at A(0,0) 798: END DO 799: S = S + S 800: * double s for the off diagonal elements 801: CALL DLASSQ( K, A( 1 ), LDA+1, SCALE, S ) 802: * tri L at A(1,0) 803: CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S ) 804: * tri U at A(0,0) 805: END IF 806: ELSE 807: * A is xpose 808: IF( ILU.EQ.0 ) THEN 809: * A' is upper 810: DO J = 1, K - 1 811: CALL DLASSQ( J, A( 0+( K+1+J )*LDA ), 1, SCALE, S ) 812: * U at A(0,k+1) 813: END DO 814: DO J = 0, K - 1 815: CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 816: * k by k rect. at A(0,0) 817: END DO 818: DO J = 0, K - 2 819: CALL DLASSQ( K-J-1, A( J+1+( J+K )*LDA ), 1, SCALE, 820: + S ) 821: * L at A(0,k) 822: END DO 823: S = S + S 824: * double s for the off diagonal elements 825: CALL DLASSQ( K, A( 0+( K+1 )*LDA ), LDA+1, SCALE, S ) 826: * tri U at A(0,k+1) 827: CALL DLASSQ( K, A( 0+K*LDA ), LDA+1, SCALE, S ) 828: * tri L at A(0,k) 829: ELSE 830: * A' is lower 831: DO J = 1, K - 1 832: CALL DLASSQ( J, A( 0+( J+1 )*LDA ), 1, SCALE, S ) 833: * U at A(0,1) 834: END DO 835: DO J = K + 1, N 836: CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 837: * k by k rect. at A(0,k+1) 838: END DO 839: DO J = 0, K - 2 840: CALL DLASSQ( K-J-1, A( J+1+J*LDA ), 1, SCALE, S ) 841: * L at A(0,0) 842: END DO 843: S = S + S 844: * double s for the off diagonal elements 845: CALL DLASSQ( K, A( LDA ), LDA+1, SCALE, S ) 846: * tri L at A(0,1) 847: CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S ) 848: * tri U at A(0,0) 849: END IF 850: END IF 851: END IF 852: VALUE = SCALE*SQRT( S ) 853: END IF 854: * 855: DLANSF = VALUE 856: RETURN 857: * 858: * End of DLANSF 859: * 860: END