Annotation of rpl/lapack/lapack/dlansf.f, revision 1.1
1.1 ! bertrand 1: DOUBLE PRECISION FUNCTION DLANSF( NORM, TRANSR, UPLO, N, A, WORK )
! 2: *
! 3: * -- LAPACK routine (version 3.2.2) --
! 4: *
! 5: * -- Contributed by Fred Gustavson of the IBM Watson Research Center --
! 6: * -- June 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
! 48: * Specifies the value to be returned in DLANSF as described
! 49: * above.
! 50: *
! 51: * TRANSR (input) CHARACTER
! 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
! 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlansf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack