![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: DOUBLE PRECISION FUNCTION ZLANHF( NORM, TRANSR, UPLO, N, A, WORK ) 2: * 3: * -- LAPACK routine (version 3.2.1) -- 4: * 5: * -- Contributed by Fred Gustavson of the IBM Watson Research Center -- 6: * -- April 2009 -- 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 WORK( 0: * ) 17: COMPLEX*16 A( 0: * ) 18: * .. 19: * 20: * Purpose 21: * ======= 22: * 23: * ZLANHF returns the value of the one norm, or the Frobenius norm, or 24: * the infinity norm, or the element of largest absolute value of a 25: * complex Hermitian matrix A in RFP format. 26: * 27: * Description 28: * =========== 29: * 30: * ZLANHF returns the value 31: * 32: * ZLANHF = ( max(abs(A(i,j))), NORM = 'M' or 'm' 33: * ( 34: * ( norm1(A), NORM = '1', 'O' or 'o' 35: * ( 36: * ( normI(A), NORM = 'I' or 'i' 37: * ( 38: * ( normF(A), NORM = 'F', 'f', 'E' or 'e' 39: * 40: * where norm1 denotes the one norm of a matrix (maximum column sum), 41: * normI denotes the infinity norm of a matrix (maximum row sum) and 42: * normF denotes the Frobenius norm of a matrix (square root of sum of 43: * squares). Note that max(abs(A(i,j))) is not a matrix norm. 44: * 45: * Arguments 46: * ========= 47: * 48: * NORM (input) CHARACTER 49: * Specifies the value to be returned in ZLANHF as described 50: * above. 51: * 52: * TRANSR (input) CHARACTER 53: * Specifies whether the RFP format of A is normal or 54: * conjugate-transposed format. 55: * = 'N': RFP format is Normal 56: * = 'C': RFP format is Conjugate-transposed 57: * 58: * UPLO (input) CHARACTER 59: * On entry, UPLO specifies whether the RFP matrix A came from 60: * an upper or lower triangular matrix as follows: 61: * 62: * UPLO = 'U' or 'u' RFP A came from an upper triangular 63: * matrix 64: * 65: * UPLO = 'L' or 'l' RFP A came from a lower triangular 66: * matrix 67: * 68: * N (input) INTEGER 69: * The order of the matrix A. N >= 0. When N = 0, ZLANHF is 70: * set to zero. 71: * 72: * A (input) COMPLEX*16 array, dimension ( N*(N+1)/2 ); 73: * On entry, the matrix A in RFP Format. 74: * RFP Format is described by TRANSR, UPLO and N as follows: 75: * If TRANSR='N' then RFP A is (0:N,0:K-1) when N is even; 76: * K=N/2. RFP A is (0:N-1,0:K) when N is odd; K=N/2. If 77: * TRANSR = 'C' then RFP is the Conjugate-transpose of RFP A 78: * as defined when TRANSR = 'N'. The contents of RFP A are 79: * defined by UPLO as follows: If UPLO = 'U' the RFP A 80: * contains the ( N*(N+1)/2 ) elements of upper packed A 81: * either in normal or conjugate-transpose Format. If 82: * UPLO = 'L' the RFP A contains the ( N*(N+1) /2 ) elements 83: * of lower packed A either in normal or conjugate-transpose 84: * Format. The LDA of RFP A is (N+1)/2 when TRANSR = 'C'. When 85: * TRANSR is 'N' the LDA is N+1 when N is even and is N when 86: * is odd. See the Note below for more details. 87: * Unchanged on exit. 88: * 89: * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK), 90: * where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, 91: * WORK is not referenced. 92: * 93: * Further Details 94: * =============== 95: * 96: * We first consider Standard Packed Format when N is even. 97: * We give an example where N = 6. 98: * 99: * AP is Upper AP is Lower 100: * 101: * 00 01 02 03 04 05 00 102: * 11 12 13 14 15 10 11 103: * 22 23 24 25 20 21 22 104: * 33 34 35 30 31 32 33 105: * 44 45 40 41 42 43 44 106: * 55 50 51 52 53 54 55 107: * 108: * 109: * Let TRANSR = 'N'. RFP holds AP as follows: 110: * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last 111: * three columns of AP upper. The lower triangle A(4:6,0:2) consists of 112: * conjugate-transpose of the first three columns of AP upper. 113: * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first 114: * three columns of AP lower. The upper triangle A(0:2,0:2) consists of 115: * conjugate-transpose of the last three columns of AP lower. 116: * To denote conjugate we place -- above the element. This covers the 117: * case N even and TRANSR = 'N'. 118: * 119: * RFP A RFP A 120: * 121: * -- -- -- 122: * 03 04 05 33 43 53 123: * -- -- 124: * 13 14 15 00 44 54 125: * -- 126: * 23 24 25 10 11 55 127: * 128: * 33 34 35 20 21 22 129: * -- 130: * 00 44 45 30 31 32 131: * -- -- 132: * 01 11 55 40 41 42 133: * -- -- -- 134: * 02 12 22 50 51 52 135: * 136: * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- 137: * transpose of RFP A above. One therefore gets: 138: * 139: * 140: * RFP A RFP A 141: * 142: * -- -- -- -- -- -- -- -- -- -- 143: * 03 13 23 33 00 01 02 33 00 10 20 30 40 50 144: * -- -- -- -- -- -- -- -- -- -- 145: * 04 14 24 34 44 11 12 43 44 11 21 31 41 51 146: * -- -- -- -- -- -- -- -- -- -- 147: * 05 15 25 35 45 55 22 53 54 55 22 32 42 52 148: * 149: * 150: * We next consider Standard Packed Format when N is odd. 151: * We give an example where N = 5. 152: * 153: * AP is Upper AP is Lower 154: * 155: * 00 01 02 03 04 00 156: * 11 12 13 14 10 11 157: * 22 23 24 20 21 22 158: * 33 34 30 31 32 33 159: * 44 40 41 42 43 44 160: * 161: * 162: * Let TRANSR = 'N'. RFP holds AP as follows: 163: * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last 164: * three columns of AP upper. The lower triangle A(3:4,0:1) consists of 165: * conjugate-transpose of the first two columns of AP upper. 166: * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first 167: * three columns of AP lower. The upper triangle A(0:1,1:2) consists of 168: * conjugate-transpose of the last two columns of AP lower. 169: * To denote conjugate we place -- above the element. This covers the 170: * case N odd and TRANSR = 'N'. 171: * 172: * RFP A RFP A 173: * 174: * -- -- 175: * 02 03 04 00 33 43 176: * -- 177: * 12 13 14 10 11 44 178: * 179: * 22 23 24 20 21 22 180: * -- 181: * 00 33 34 30 31 32 182: * -- -- 183: * 01 11 44 40 41 42 184: * 185: * Now let TRANSR = 'C'. RFP A in both UPLO cases is just the conjugate- 186: * transpose of RFP A above. One therefore gets: 187: * 188: * 189: * RFP A RFP A 190: * 191: * -- -- -- -- -- -- -- -- -- 192: * 02 12 22 00 01 00 10 20 30 40 50 193: * -- -- -- -- -- -- -- -- -- 194: * 03 13 23 33 11 33 11 21 31 41 51 195: * -- -- -- -- -- -- -- -- -- 196: * 04 14 24 34 44 43 44 22 32 42 52 197: * 198: * ===================================================================== 199: * 200: * .. Parameters .. 201: DOUBLE PRECISION ONE, ZERO 202: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 203: * .. 204: * .. Local Scalars .. 205: INTEGER I, J, IFM, ILU, NOE, N1, K, L, LDA 206: DOUBLE PRECISION SCALE, S, VALUE, AA 207: * .. 208: * .. External Functions .. 209: LOGICAL LSAME 210: INTEGER IDAMAX 211: EXTERNAL LSAME, IDAMAX 212: * .. 213: * .. External Subroutines .. 214: EXTERNAL ZLASSQ 215: * .. 216: * .. Intrinsic Functions .. 217: INTRINSIC ABS, DBLE, MAX, SQRT 218: * .. 219: * .. Executable Statements .. 220: * 221: IF( N.EQ.0 ) THEN 222: ZLANHF = ZERO 223: RETURN 224: END IF 225: * 226: * set noe = 1 if n is odd. if n is even set noe=0 227: * 228: NOE = 1 229: IF( MOD( N, 2 ).EQ.0 ) 230: + NOE = 0 231: * 232: * set ifm = 0 when form='C' or 'c' and 1 otherwise 233: * 234: IFM = 1 235: IF( LSAME( TRANSR, 'C' ) ) 236: + IFM = 0 237: * 238: * set ilu = 0 when uplo='U or 'u' and 1 otherwise 239: * 240: ILU = 1 241: IF( LSAME( UPLO, 'U' ) ) 242: + ILU = 0 243: * 244: * set lda = (n+1)/2 when ifm = 0 245: * set lda = n when ifm = 1 and noe = 1 246: * set lda = n+1 when ifm = 1 and noe = 0 247: * 248: IF( IFM.EQ.1 ) THEN 249: IF( NOE.EQ.1 ) THEN 250: LDA = N 251: ELSE 252: * noe=0 253: LDA = N + 1 254: END IF 255: ELSE 256: * ifm=0 257: LDA = ( N+1 ) / 2 258: END IF 259: * 260: IF( LSAME( NORM, 'M' ) ) THEN 261: * 262: * Find max(abs(A(i,j))). 263: * 264: K = ( N+1 ) / 2 265: VALUE = ZERO 266: IF( NOE.EQ.1 ) THEN 267: * n is odd & n = k + k - 1 268: IF( IFM.EQ.1 ) THEN 269: * A is n by k 270: IF( ILU.EQ.1 ) THEN 271: * uplo ='L' 272: J = 0 273: * -> L(0,0) 274: VALUE = MAX( VALUE, ABS( DBLE( A( J+J*LDA ) ) ) ) 275: DO I = 1, N - 1 276: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 277: END DO 278: DO J = 1, K - 1 279: DO I = 0, J - 2 280: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 281: END DO 282: I = J - 1 283: * L(k+j,k+j) 284: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 285: I = J 286: * -> L(j,j) 287: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 288: DO I = J + 1, N - 1 289: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 290: END DO 291: END DO 292: ELSE 293: * uplo = 'U' 294: DO J = 0, K - 2 295: DO I = 0, K + J - 2 296: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 297: END DO 298: I = K + J - 1 299: * -> U(i,i) 300: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 301: I = I + 1 302: * =k+j; i -> U(j,j) 303: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 304: DO I = K + J + 1, N - 1 305: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 306: END DO 307: END DO 308: DO I = 0, N - 2 309: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 310: * j=k-1 311: END DO 312: * i=n-1 -> U(n-1,n-1) 313: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 314: END IF 315: ELSE 316: * xpose case; A is k by n 317: IF( ILU.EQ.1 ) THEN 318: * uplo ='L' 319: DO J = 0, K - 2 320: DO I = 0, J - 1 321: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 322: END DO 323: I = J 324: * L(i,i) 325: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 326: I = J + 1 327: * L(j+k,j+k) 328: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 329: DO I = J + 2, K - 1 330: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 331: END DO 332: END DO 333: J = K - 1 334: DO I = 0, K - 2 335: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 336: END DO 337: I = K - 1 338: * -> L(i,i) is at A(i,j) 339: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 340: DO J = K, N - 1 341: DO I = 0, K - 1 342: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 343: END DO 344: END DO 345: ELSE 346: * uplo = 'U' 347: DO J = 0, K - 2 348: DO I = 0, K - 1 349: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 350: END DO 351: END DO 352: J = K - 1 353: * -> U(j,j) is at A(0,j) 354: VALUE = MAX( VALUE, ABS( DBLE( A( 0+J*LDA ) ) ) ) 355: DO I = 1, K - 1 356: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 357: END DO 358: DO J = K, N - 1 359: DO I = 0, J - K - 1 360: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 361: END DO 362: I = J - K 363: * -> U(i,i) at A(i,j) 364: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 365: I = J - K + 1 366: * U(j,j) 367: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 368: DO I = J - K + 2, K - 1 369: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 370: END DO 371: END DO 372: END IF 373: END IF 374: ELSE 375: * n is even & k = n/2 376: IF( IFM.EQ.1 ) THEN 377: * A is n+1 by k 378: IF( ILU.EQ.1 ) THEN 379: * uplo ='L' 380: J = 0 381: * -> L(k,k) & j=1 -> L(0,0) 382: VALUE = MAX( VALUE, ABS( DBLE( A( J+J*LDA ) ) ) ) 383: VALUE = MAX( VALUE, ABS( DBLE( A( J+1+J*LDA ) ) ) ) 384: DO I = 2, N 385: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 386: END DO 387: DO J = 1, K - 1 388: DO I = 0, J - 1 389: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 390: END DO 391: I = J 392: * L(k+j,k+j) 393: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 394: I = J + 1 395: * -> L(j,j) 396: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 397: DO I = J + 2, N 398: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 399: END DO 400: END DO 401: ELSE 402: * uplo = 'U' 403: DO J = 0, K - 2 404: DO I = 0, K + J - 1 405: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 406: END DO 407: I = K + J 408: * -> U(i,i) 409: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 410: I = I + 1 411: * =k+j+1; i -> U(j,j) 412: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 413: DO I = K + J + 2, N 414: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 415: END DO 416: END DO 417: DO I = 0, N - 2 418: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 419: * j=k-1 420: END DO 421: * i=n-1 -> U(n-1,n-1) 422: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 423: I = N 424: * -> U(k-1,k-1) 425: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 426: END IF 427: ELSE 428: * xpose case; A is k by n+1 429: IF( ILU.EQ.1 ) THEN 430: * uplo ='L' 431: J = 0 432: * -> L(k,k) at A(0,0) 433: VALUE = MAX( VALUE, ABS( DBLE( A( J+J*LDA ) ) ) ) 434: DO I = 1, K - 1 435: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 436: END DO 437: DO J = 1, K - 1 438: DO I = 0, J - 2 439: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 440: END DO 441: I = J - 1 442: * L(i,i) 443: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 444: I = J 445: * L(j+k,j+k) 446: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 447: DO I = J + 1, K - 1 448: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 449: END DO 450: END DO 451: J = K 452: DO I = 0, K - 2 453: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 454: END DO 455: I = K - 1 456: * -> L(i,i) is at A(i,j) 457: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 458: DO J = K + 1, N 459: DO I = 0, K - 1 460: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 461: END DO 462: END DO 463: ELSE 464: * uplo = 'U' 465: DO J = 0, K - 1 466: DO I = 0, K - 1 467: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 468: END DO 469: END DO 470: J = K 471: * -> U(j,j) is at A(0,j) 472: VALUE = MAX( VALUE, ABS( DBLE( A( 0+J*LDA ) ) ) ) 473: DO I = 1, K - 1 474: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 475: END DO 476: DO J = K + 1, N - 1 477: DO I = 0, J - K - 2 478: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 479: END DO 480: I = J - K - 1 481: * -> U(i,i) at A(i,j) 482: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 483: I = J - K 484: * U(j,j) 485: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 486: DO I = J - K + 1, K - 1 487: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 488: END DO 489: END DO 490: J = N 491: DO I = 0, K - 2 492: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) ) 493: END DO 494: I = K - 1 495: * U(k,k) at A(i,j) 496: VALUE = MAX( VALUE, ABS( DBLE( A( I+J*LDA ) ) ) ) 497: END IF 498: END IF 499: END IF 500: ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR. 501: + ( NORM.EQ.'1' ) ) THEN 502: * 503: * Find normI(A) ( = norm1(A), since A is Hermitian). 504: * 505: IF( IFM.EQ.1 ) THEN 506: * A is 'N' 507: K = N / 2 508: IF( NOE.EQ.1 ) THEN 509: * n is odd & A is n by (n+1)/2 510: IF( ILU.EQ.0 ) THEN 511: * uplo = 'U' 512: DO I = 0, K - 1 513: WORK( I ) = ZERO 514: END DO 515: DO J = 0, K 516: S = ZERO 517: DO I = 0, K + J - 1 518: AA = ABS( A( I+J*LDA ) ) 519: * -> A(i,j+k) 520: S = S + AA 521: WORK( I ) = WORK( I ) + AA 522: END DO 523: AA = ABS( DBLE( A( I+J*LDA ) ) ) 524: * -> A(j+k,j+k) 525: WORK( J+K ) = S + AA 526: IF( I.EQ.K+K ) 527: + GO TO 10 528: I = I + 1 529: AA = ABS( DBLE( A( I+J*LDA ) ) ) 530: * -> A(j,j) 531: WORK( J ) = WORK( J ) + AA 532: S = ZERO 533: DO L = J + 1, K - 1 534: I = I + 1 535: AA = ABS( A( I+J*LDA ) ) 536: * -> A(l,j) 537: S = S + AA 538: WORK( L ) = WORK( L ) + AA 539: END DO 540: WORK( J ) = WORK( J ) + S 541: END DO 542: 10 CONTINUE 543: I = IDAMAX( N, WORK, 1 ) 544: VALUE = WORK( I-1 ) 545: ELSE 546: * ilu = 1 & uplo = 'L' 547: K = K + 1 548: * k=(n+1)/2 for n odd and ilu=1 549: DO I = K, N - 1 550: WORK( I ) = ZERO 551: END DO 552: DO J = K - 1, 0, -1 553: S = ZERO 554: DO I = 0, J - 2 555: AA = ABS( A( I+J*LDA ) ) 556: * -> A(j+k,i+k) 557: S = S + AA 558: WORK( I+K ) = WORK( I+K ) + AA 559: END DO 560: IF( J.GT.0 ) THEN 561: AA = ABS( DBLE( A( I+J*LDA ) ) ) 562: * -> A(j+k,j+k) 563: S = S + AA 564: WORK( I+K ) = WORK( I+K ) + S 565: * i=j 566: I = I + 1 567: END IF 568: AA = ABS( DBLE( A( I+J*LDA ) ) ) 569: * -> A(j,j) 570: WORK( J ) = AA 571: S = ZERO 572: DO L = J + 1, N - 1 573: I = I + 1 574: AA = ABS( A( I+J*LDA ) ) 575: * -> A(l,j) 576: S = S + AA 577: WORK( L ) = WORK( L ) + AA 578: END DO 579: WORK( J ) = WORK( J ) + S 580: END DO 581: I = IDAMAX( N, WORK, 1 ) 582: VALUE = WORK( I-1 ) 583: END IF 584: ELSE 585: * n is even & A is n+1 by k = n/2 586: IF( ILU.EQ.0 ) THEN 587: * uplo = 'U' 588: DO I = 0, K - 1 589: WORK( I ) = ZERO 590: END DO 591: DO J = 0, K - 1 592: S = ZERO 593: DO I = 0, K + J - 1 594: AA = ABS( A( I+J*LDA ) ) 595: * -> A(i,j+k) 596: S = S + AA 597: WORK( I ) = WORK( I ) + AA 598: END DO 599: AA = ABS( DBLE( A( I+J*LDA ) ) ) 600: * -> A(j+k,j+k) 601: WORK( J+K ) = S + AA 602: I = I + 1 603: AA = ABS( DBLE( A( I+J*LDA ) ) ) 604: * -> A(j,j) 605: WORK( J ) = WORK( J ) + AA 606: S = ZERO 607: DO L = J + 1, K - 1 608: I = I + 1 609: AA = ABS( A( I+J*LDA ) ) 610: * -> A(l,j) 611: S = S + AA 612: WORK( L ) = WORK( L ) + AA 613: END DO 614: WORK( J ) = WORK( J ) + S 615: END DO 616: I = IDAMAX( N, WORK, 1 ) 617: VALUE = WORK( I-1 ) 618: ELSE 619: * ilu = 1 & uplo = 'L' 620: DO I = K, N - 1 621: WORK( I ) = ZERO 622: END DO 623: DO J = K - 1, 0, -1 624: S = ZERO 625: DO I = 0, J - 1 626: AA = ABS( A( I+J*LDA ) ) 627: * -> A(j+k,i+k) 628: S = S + AA 629: WORK( I+K ) = WORK( I+K ) + AA 630: END DO 631: AA = ABS( DBLE( A( I+J*LDA ) ) ) 632: * -> A(j+k,j+k) 633: S = S + AA 634: WORK( I+K ) = WORK( I+K ) + S 635: * i=j 636: I = I + 1 637: AA = ABS( DBLE( A( I+J*LDA ) ) ) 638: * -> A(j,j) 639: WORK( J ) = AA 640: S = ZERO 641: DO L = J + 1, N - 1 642: I = I + 1 643: AA = ABS( A( I+J*LDA ) ) 644: * -> A(l,j) 645: S = S + AA 646: WORK( L ) = WORK( L ) + AA 647: END DO 648: WORK( J ) = WORK( J ) + S 649: END DO 650: I = IDAMAX( N, WORK, 1 ) 651: VALUE = WORK( I-1 ) 652: END IF 653: END IF 654: ELSE 655: * ifm=0 656: K = N / 2 657: IF( NOE.EQ.1 ) THEN 658: * n is odd & A is (n+1)/2 by n 659: IF( ILU.EQ.0 ) THEN 660: * uplo = 'U' 661: N1 = K 662: * n/2 663: K = K + 1 664: * k is the row size and lda 665: DO I = N1, N - 1 666: WORK( I ) = ZERO 667: END DO 668: DO J = 0, N1 - 1 669: S = ZERO 670: DO I = 0, K - 1 671: AA = ABS( A( I+J*LDA ) ) 672: * A(j,n1+i) 673: WORK( I+N1 ) = WORK( I+N1 ) + AA 674: S = S + AA 675: END DO 676: WORK( J ) = S 677: END DO 678: * j=n1=k-1 is special 679: S = ABS( DBLE( A( 0+J*LDA ) ) ) 680: * A(k-1,k-1) 681: DO I = 1, K - 1 682: AA = ABS( A( I+J*LDA ) ) 683: * A(k-1,i+n1) 684: WORK( I+N1 ) = WORK( I+N1 ) + AA 685: S = S + AA 686: END DO 687: WORK( J ) = WORK( J ) + S 688: DO J = K, N - 1 689: S = ZERO 690: DO I = 0, J - K - 1 691: AA = ABS( A( I+J*LDA ) ) 692: * A(i,j-k) 693: WORK( I ) = WORK( I ) + AA 694: S = S + AA 695: END DO 696: * i=j-k 697: AA = ABS( DBLE( A( I+J*LDA ) ) ) 698: * A(j-k,j-k) 699: S = S + AA 700: WORK( J-K ) = WORK( J-K ) + S 701: I = I + 1 702: S = ABS( DBLE( A( I+J*LDA ) ) ) 703: * A(j,j) 704: DO L = J + 1, N - 1 705: I = I + 1 706: AA = ABS( A( I+J*LDA ) ) 707: * A(j,l) 708: WORK( L ) = WORK( L ) + AA 709: S = S + AA 710: END DO 711: WORK( J ) = WORK( J ) + S 712: END DO 713: I = IDAMAX( N, WORK, 1 ) 714: VALUE = WORK( I-1 ) 715: ELSE 716: * ilu=1 & uplo = 'L' 717: K = K + 1 718: * k=(n+1)/2 for n odd and ilu=1 719: DO I = K, N - 1 720: WORK( I ) = ZERO 721: END DO 722: DO J = 0, K - 2 723: * process 724: S = ZERO 725: DO I = 0, J - 1 726: AA = ABS( A( I+J*LDA ) ) 727: * A(j,i) 728: WORK( I ) = WORK( I ) + AA 729: S = S + AA 730: END DO 731: AA = ABS( DBLE( A( I+J*LDA ) ) ) 732: * i=j so process of A(j,j) 733: S = S + AA 734: WORK( J ) = S 735: * is initialised here 736: I = I + 1 737: * i=j process A(j+k,j+k) 738: AA = ABS( DBLE( A( I+J*LDA ) ) ) 739: S = AA 740: DO L = K + J + 1, N - 1 741: I = I + 1 742: AA = ABS( A( I+J*LDA ) ) 743: * A(l,k+j) 744: S = S + AA 745: WORK( L ) = WORK( L ) + AA 746: END DO 747: WORK( K+J ) = WORK( K+J ) + S 748: END DO 749: * j=k-1 is special :process col A(k-1,0:k-1) 750: S = ZERO 751: DO I = 0, K - 2 752: AA = ABS( A( I+J*LDA ) ) 753: * A(k,i) 754: WORK( I ) = WORK( I ) + AA 755: S = S + AA 756: END DO 757: * i=k-1 758: AA = ABS( DBLE( A( I+J*LDA ) ) ) 759: * A(k-1,k-1) 760: S = S + AA 761: WORK( I ) = S 762: * done with col j=k+1 763: DO J = K, N - 1 764: * process col j of A = A(j,0:k-1) 765: S = ZERO 766: DO I = 0, K - 1 767: AA = ABS( A( I+J*LDA ) ) 768: * A(j,i) 769: WORK( I ) = WORK( I ) + AA 770: S = S + AA 771: END DO 772: WORK( J ) = WORK( J ) + S 773: END DO 774: I = IDAMAX( N, WORK, 1 ) 775: VALUE = WORK( I-1 ) 776: END IF 777: ELSE 778: * n is even & A is k=n/2 by n+1 779: IF( ILU.EQ.0 ) THEN 780: * uplo = 'U' 781: DO I = K, N - 1 782: WORK( I ) = ZERO 783: END DO 784: DO J = 0, K - 1 785: S = ZERO 786: DO I = 0, K - 1 787: AA = ABS( A( I+J*LDA ) ) 788: * A(j,i+k) 789: WORK( I+K ) = WORK( I+K ) + AA 790: S = S + AA 791: END DO 792: WORK( J ) = S 793: END DO 794: * j=k 795: AA = ABS( DBLE( A( 0+J*LDA ) ) ) 796: * A(k,k) 797: S = AA 798: DO I = 1, K - 1 799: AA = ABS( A( I+J*LDA ) ) 800: * A(k,k+i) 801: WORK( I+K ) = WORK( I+K ) + AA 802: S = S + AA 803: END DO 804: WORK( J ) = WORK( J ) + S 805: DO J = K + 1, N - 1 806: S = ZERO 807: DO I = 0, J - 2 - K 808: AA = ABS( A( I+J*LDA ) ) 809: * A(i,j-k-1) 810: WORK( I ) = WORK( I ) + AA 811: S = S + AA 812: END DO 813: * i=j-1-k 814: AA = ABS( DBLE( A( I+J*LDA ) ) ) 815: * A(j-k-1,j-k-1) 816: S = S + AA 817: WORK( J-K-1 ) = WORK( J-K-1 ) + S 818: I = I + 1 819: AA = ABS( DBLE( A( I+J*LDA ) ) ) 820: * A(j,j) 821: S = AA 822: DO L = J + 1, N - 1 823: I = I + 1 824: AA = ABS( A( I+J*LDA ) ) 825: * A(j,l) 826: WORK( L ) = WORK( L ) + AA 827: S = S + AA 828: END DO 829: WORK( J ) = WORK( J ) + S 830: END DO 831: * j=n 832: S = ZERO 833: DO I = 0, K - 2 834: AA = ABS( A( I+J*LDA ) ) 835: * A(i,k-1) 836: WORK( I ) = WORK( I ) + AA 837: S = S + AA 838: END DO 839: * i=k-1 840: AA = ABS( DBLE( A( I+J*LDA ) ) ) 841: * A(k-1,k-1) 842: S = S + AA 843: WORK( I ) = WORK( I ) + S 844: I = IDAMAX( N, WORK, 1 ) 845: VALUE = WORK( I-1 ) 846: ELSE 847: * ilu=1 & uplo = 'L' 848: DO I = K, N - 1 849: WORK( I ) = ZERO 850: END DO 851: * j=0 is special :process col A(k:n-1,k) 852: S = ABS( DBLE( A( 0 ) ) ) 853: * A(k,k) 854: DO I = 1, K - 1 855: AA = ABS( A( I ) ) 856: * A(k+i,k) 857: WORK( I+K ) = WORK( I+K ) + AA 858: S = S + AA 859: END DO 860: WORK( K ) = WORK( K ) + S 861: DO J = 1, K - 1 862: * process 863: S = ZERO 864: DO I = 0, J - 2 865: AA = ABS( A( I+J*LDA ) ) 866: * A(j-1,i) 867: WORK( I ) = WORK( I ) + AA 868: S = S + AA 869: END DO 870: AA = ABS( DBLE( A( I+J*LDA ) ) ) 871: * i=j-1 so process of A(j-1,j-1) 872: S = S + AA 873: WORK( J-1 ) = S 874: * is initialised here 875: I = I + 1 876: * i=j process A(j+k,j+k) 877: AA = ABS( DBLE( A( I+J*LDA ) ) ) 878: S = AA 879: DO L = K + J + 1, N - 1 880: I = I + 1 881: AA = ABS( A( I+J*LDA ) ) 882: * A(l,k+j) 883: S = S + AA 884: WORK( L ) = WORK( L ) + AA 885: END DO 886: WORK( K+J ) = WORK( K+J ) + S 887: END DO 888: * j=k is special :process col A(k,0:k-1) 889: S = ZERO 890: DO I = 0, K - 2 891: AA = ABS( A( I+J*LDA ) ) 892: * A(k,i) 893: WORK( I ) = WORK( I ) + AA 894: S = S + AA 895: END DO 896: * 897: * i=k-1 898: AA = ABS( DBLE( A( I+J*LDA ) ) ) 899: * A(k-1,k-1) 900: S = S + AA 901: WORK( I ) = S 902: * done with col j=k+1 903: DO J = K + 1, N 904: * 905: * process col j-1 of A = A(j-1,0:k-1) 906: S = ZERO 907: DO I = 0, K - 1 908: AA = ABS( A( I+J*LDA ) ) 909: * A(j-1,i) 910: WORK( I ) = WORK( I ) + AA 911: S = S + AA 912: END DO 913: WORK( J-1 ) = WORK( J-1 ) + S 914: END DO 915: I = IDAMAX( N, WORK, 1 ) 916: VALUE = WORK( I-1 ) 917: END IF 918: END IF 919: END IF 920: ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN 921: * 922: * Find normF(A). 923: * 924: K = ( N+1 ) / 2 925: SCALE = ZERO 926: S = ONE 927: IF( NOE.EQ.1 ) THEN 928: * n is odd 929: IF( IFM.EQ.1 ) THEN 930: * A is normal & A is n by k 931: IF( ILU.EQ.0 ) THEN 932: * A is upper 933: DO J = 0, K - 3 934: CALL ZLASSQ( K-J-2, A( K+J+1+J*LDA ), 1, SCALE, S ) 935: * L at A(k,0) 936: END DO 937: DO J = 0, K - 1 938: CALL ZLASSQ( K+J-1, A( 0+J*LDA ), 1, SCALE, S ) 939: * trap U at A(0,0) 940: END DO 941: S = S + S 942: * double s for the off diagonal elements 943: L = K - 1 944: * -> U(k,k) at A(k-1,0) 945: DO I = 0, K - 2 946: AA = DBLE( A( L ) ) 947: * U(k+i,k+i) 948: IF( AA.NE.ZERO ) THEN 949: IF( SCALE.LT.AA ) THEN 950: S = ONE + S*( SCALE / AA )**2 951: SCALE = AA 952: ELSE 953: S = S + ( AA / SCALE )**2 954: END IF 955: END IF 956: AA = DBLE( A( L+1 ) ) 957: * U(i,i) 958: IF( AA.NE.ZERO ) THEN 959: IF( SCALE.LT.AA ) THEN 960: S = ONE + S*( SCALE / AA )**2 961: SCALE = AA 962: ELSE 963: S = S + ( AA / SCALE )**2 964: END IF 965: END IF 966: L = L + LDA + 1 967: END DO 968: AA = DBLE( A( L ) ) 969: * U(n-1,n-1) 970: IF( AA.NE.ZERO ) THEN 971: IF( SCALE.LT.AA ) THEN 972: S = ONE + S*( SCALE / AA )**2 973: SCALE = AA 974: ELSE 975: S = S + ( AA / SCALE )**2 976: END IF 977: END IF 978: ELSE 979: * ilu=1 & A is lower 980: DO J = 0, K - 1 981: CALL ZLASSQ( N-J-1, A( J+1+J*LDA ), 1, SCALE, S ) 982: * trap L at A(0,0) 983: END DO 984: DO J = 1, K - 2 985: CALL ZLASSQ( J, A( 0+( 1+J )*LDA ), 1, SCALE, S ) 986: * U at A(0,1) 987: END DO 988: S = S + S 989: * double s for the off diagonal elements 990: AA = DBLE( A( 0 ) ) 991: * L(0,0) at A(0,0) 992: IF( AA.NE.ZERO ) THEN 993: IF( SCALE.LT.AA ) THEN 994: S = ONE + S*( SCALE / AA )**2 995: SCALE = AA 996: ELSE 997: S = S + ( AA / SCALE )**2 998: END IF 999: END IF 1000: L = LDA 1001: * -> L(k,k) at A(0,1) 1002: DO I = 1, K - 1 1003: AA = DBLE( A( L ) ) 1004: * L(k-1+i,k-1+i) 1005: IF( AA.NE.ZERO ) THEN 1006: IF( SCALE.LT.AA ) THEN 1007: S = ONE + S*( SCALE / AA )**2 1008: SCALE = AA 1009: ELSE 1010: S = S + ( AA / SCALE )**2 1011: END IF 1012: END IF 1013: AA = DBLE( A( L+1 ) ) 1014: * L(i,i) 1015: IF( AA.NE.ZERO ) THEN 1016: IF( SCALE.LT.AA ) THEN 1017: S = ONE + S*( SCALE / AA )**2 1018: SCALE = AA 1019: ELSE 1020: S = S + ( AA / SCALE )**2 1021: END IF 1022: END IF 1023: L = L + LDA + 1 1024: END DO 1025: END IF 1026: ELSE 1027: * A is xpose & A is k by n 1028: IF( ILU.EQ.0 ) THEN 1029: * A' is upper 1030: DO J = 1, K - 2 1031: CALL ZLASSQ( J, A( 0+( K+J )*LDA ), 1, SCALE, S ) 1032: * U at A(0,k) 1033: END DO 1034: DO J = 0, K - 2 1035: CALL ZLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 1036: * k by k-1 rect. at A(0,0) 1037: END DO 1038: DO J = 0, K - 2 1039: CALL ZLASSQ( K-J-1, A( J+1+( J+K-1 )*LDA ), 1, 1040: + SCALE, S ) 1041: * L at A(0,k-1) 1042: END DO 1043: S = S + S 1044: * double s for the off diagonal elements 1045: L = 0 + K*LDA - LDA 1046: * -> U(k-1,k-1) at A(0,k-1) 1047: AA = DBLE( A( L ) ) 1048: * U(k-1,k-1) 1049: IF( AA.NE.ZERO ) THEN 1050: IF( SCALE.LT.AA ) THEN 1051: S = ONE + S*( SCALE / AA )**2 1052: SCALE = AA 1053: ELSE 1054: S = S + ( AA / SCALE )**2 1055: END IF 1056: END IF 1057: L = L + LDA 1058: * -> U(0,0) at A(0,k) 1059: DO J = K, N - 1 1060: AA = DBLE( A( L ) ) 1061: * -> U(j-k,j-k) 1062: IF( AA.NE.ZERO ) THEN 1063: IF( SCALE.LT.AA ) THEN 1064: S = ONE + S*( SCALE / AA )**2 1065: SCALE = AA 1066: ELSE 1067: S = S + ( AA / SCALE )**2 1068: END IF 1069: END IF 1070: AA = DBLE( A( L+1 ) ) 1071: * -> U(j,j) 1072: IF( AA.NE.ZERO ) THEN 1073: IF( SCALE.LT.AA ) THEN 1074: S = ONE + S*( SCALE / AA )**2 1075: SCALE = AA 1076: ELSE 1077: S = S + ( AA / SCALE )**2 1078: END IF 1079: END IF 1080: L = L + LDA + 1 1081: END DO 1082: ELSE 1083: * A' is lower 1084: DO J = 1, K - 1 1085: CALL ZLASSQ( J, A( 0+J*LDA ), 1, SCALE, S ) 1086: * U at A(0,0) 1087: END DO 1088: DO J = K, N - 1 1089: CALL ZLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 1090: * k by k-1 rect. at A(0,k) 1091: END DO 1092: DO J = 0, K - 3 1093: CALL ZLASSQ( K-J-2, A( J+2+J*LDA ), 1, SCALE, S ) 1094: * L at A(1,0) 1095: END DO 1096: S = S + S 1097: * double s for the off diagonal elements 1098: L = 0 1099: * -> L(0,0) at A(0,0) 1100: DO I = 0, K - 2 1101: AA = DBLE( A( L ) ) 1102: * L(i,i) 1103: IF( AA.NE.ZERO ) THEN 1104: IF( SCALE.LT.AA ) THEN 1105: S = ONE + S*( SCALE / AA )**2 1106: SCALE = AA 1107: ELSE 1108: S = S + ( AA / SCALE )**2 1109: END IF 1110: END IF 1111: AA = DBLE( A( L+1 ) ) 1112: * L(k+i,k+i) 1113: IF( AA.NE.ZERO ) THEN 1114: IF( SCALE.LT.AA ) THEN 1115: S = ONE + S*( SCALE / AA )**2 1116: SCALE = AA 1117: ELSE 1118: S = S + ( AA / SCALE )**2 1119: END IF 1120: END IF 1121: L = L + LDA + 1 1122: END DO 1123: * L-> k-1 + (k-1)*lda or L(k-1,k-1) at A(k-1,k-1) 1124: AA = DBLE( A( L ) ) 1125: * L(k-1,k-1) at A(k-1,k-1) 1126: IF( AA.NE.ZERO ) THEN 1127: IF( SCALE.LT.AA ) THEN 1128: S = ONE + S*( SCALE / AA )**2 1129: SCALE = AA 1130: ELSE 1131: S = S + ( AA / SCALE )**2 1132: END IF 1133: END IF 1134: END IF 1135: END IF 1136: ELSE 1137: * n is even 1138: IF( IFM.EQ.1 ) THEN 1139: * A is normal 1140: IF( ILU.EQ.0 ) THEN 1141: * A is upper 1142: DO J = 0, K - 2 1143: CALL ZLASSQ( K-J-1, A( K+J+2+J*LDA ), 1, SCALE, S ) 1144: * L at A(k+1,0) 1145: END DO 1146: DO J = 0, K - 1 1147: CALL ZLASSQ( K+J, A( 0+J*LDA ), 1, SCALE, S ) 1148: * trap U at A(0,0) 1149: END DO 1150: S = S + S 1151: * double s for the off diagonal elements 1152: L = K 1153: * -> U(k,k) at A(k,0) 1154: DO I = 0, K - 1 1155: AA = DBLE( A( L ) ) 1156: * U(k+i,k+i) 1157: IF( AA.NE.ZERO ) THEN 1158: IF( SCALE.LT.AA ) THEN 1159: S = ONE + S*( SCALE / AA )**2 1160: SCALE = AA 1161: ELSE 1162: S = S + ( AA / SCALE )**2 1163: END IF 1164: END IF 1165: AA = DBLE( A( L+1 ) ) 1166: * U(i,i) 1167: IF( AA.NE.ZERO ) THEN 1168: IF( SCALE.LT.AA ) THEN 1169: S = ONE + S*( SCALE / AA )**2 1170: SCALE = AA 1171: ELSE 1172: S = S + ( AA / SCALE )**2 1173: END IF 1174: END IF 1175: L = L + LDA + 1 1176: END DO 1177: ELSE 1178: * ilu=1 & A is lower 1179: DO J = 0, K - 1 1180: CALL ZLASSQ( N-J-1, A( J+2+J*LDA ), 1, SCALE, S ) 1181: * trap L at A(1,0) 1182: END DO 1183: DO J = 1, K - 1 1184: CALL ZLASSQ( J, A( 0+J*LDA ), 1, SCALE, S ) 1185: * U at A(0,0) 1186: END DO 1187: S = S + S 1188: * double s for the off diagonal elements 1189: L = 0 1190: * -> L(k,k) at A(0,0) 1191: DO I = 0, K - 1 1192: AA = DBLE( A( L ) ) 1193: * L(k-1+i,k-1+i) 1194: IF( AA.NE.ZERO ) THEN 1195: IF( SCALE.LT.AA ) THEN 1196: S = ONE + S*( SCALE / AA )**2 1197: SCALE = AA 1198: ELSE 1199: S = S + ( AA / SCALE )**2 1200: END IF 1201: END IF 1202: AA = DBLE( A( L+1 ) ) 1203: * L(i,i) 1204: IF( AA.NE.ZERO ) THEN 1205: IF( SCALE.LT.AA ) THEN 1206: S = ONE + S*( SCALE / AA )**2 1207: SCALE = AA 1208: ELSE 1209: S = S + ( AA / SCALE )**2 1210: END IF 1211: END IF 1212: L = L + LDA + 1 1213: END DO 1214: END IF 1215: ELSE 1216: * A is xpose 1217: IF( ILU.EQ.0 ) THEN 1218: * A' is upper 1219: DO J = 1, K - 1 1220: CALL ZLASSQ( J, A( 0+( K+1+J )*LDA ), 1, SCALE, S ) 1221: * U at A(0,k+1) 1222: END DO 1223: DO J = 0, K - 1 1224: CALL ZLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 1225: * k by k rect. at A(0,0) 1226: END DO 1227: DO J = 0, K - 2 1228: CALL ZLASSQ( K-J-1, A( J+1+( J+K )*LDA ), 1, SCALE, 1229: + S ) 1230: * L at A(0,k) 1231: END DO 1232: S = S + S 1233: * double s for the off diagonal elements 1234: L = 0 + K*LDA 1235: * -> U(k,k) at A(0,k) 1236: AA = DBLE( A( L ) ) 1237: * U(k,k) 1238: IF( AA.NE.ZERO ) THEN 1239: IF( SCALE.LT.AA ) THEN 1240: S = ONE + S*( SCALE / AA )**2 1241: SCALE = AA 1242: ELSE 1243: S = S + ( AA / SCALE )**2 1244: END IF 1245: END IF 1246: L = L + LDA 1247: * -> U(0,0) at A(0,k+1) 1248: DO J = K + 1, N - 1 1249: AA = DBLE( A( L ) ) 1250: * -> U(j-k-1,j-k-1) 1251: IF( AA.NE.ZERO ) THEN 1252: IF( SCALE.LT.AA ) THEN 1253: S = ONE + S*( SCALE / AA )**2 1254: SCALE = AA 1255: ELSE 1256: S = S + ( AA / SCALE )**2 1257: END IF 1258: END IF 1259: AA = DBLE( A( L+1 ) ) 1260: * -> U(j,j) 1261: IF( AA.NE.ZERO ) THEN 1262: IF( SCALE.LT.AA ) THEN 1263: S = ONE + S*( SCALE / AA )**2 1264: SCALE = AA 1265: ELSE 1266: S = S + ( AA / SCALE )**2 1267: END IF 1268: END IF 1269: L = L + LDA + 1 1270: END DO 1271: * L=k-1+n*lda 1272: * -> U(k-1,k-1) at A(k-1,n) 1273: AA = DBLE( A( L ) ) 1274: * U(k,k) 1275: IF( AA.NE.ZERO ) THEN 1276: IF( SCALE.LT.AA ) THEN 1277: S = ONE + S*( SCALE / AA )**2 1278: SCALE = AA 1279: ELSE 1280: S = S + ( AA / SCALE )**2 1281: END IF 1282: END IF 1283: ELSE 1284: * A' is lower 1285: DO J = 1, K - 1 1286: CALL ZLASSQ( J, A( 0+( J+1 )*LDA ), 1, SCALE, S ) 1287: * U at A(0,1) 1288: END DO 1289: DO J = K + 1, N 1290: CALL ZLASSQ( K, A( 0+J*LDA ), 1, SCALE, S ) 1291: * k by k rect. at A(0,k+1) 1292: END DO 1293: DO J = 0, K - 2 1294: CALL ZLASSQ( K-J-1, A( J+1+J*LDA ), 1, SCALE, S ) 1295: * L at A(0,0) 1296: END DO 1297: S = S + S 1298: * double s for the off diagonal elements 1299: L = 0 1300: * -> L(k,k) at A(0,0) 1301: AA = DBLE( A( L ) ) 1302: * L(k,k) at A(0,0) 1303: IF( AA.NE.ZERO ) THEN 1304: IF( SCALE.LT.AA ) THEN 1305: S = ONE + S*( SCALE / AA )**2 1306: SCALE = AA 1307: ELSE 1308: S = S + ( AA / SCALE )**2 1309: END IF 1310: END IF 1311: L = LDA 1312: * -> L(0,0) at A(0,1) 1313: DO I = 0, K - 2 1314: AA = DBLE( A( L ) ) 1315: * L(i,i) 1316: IF( AA.NE.ZERO ) THEN 1317: IF( SCALE.LT.AA ) THEN 1318: S = ONE + S*( SCALE / AA )**2 1319: SCALE = AA 1320: ELSE 1321: S = S + ( AA / SCALE )**2 1322: END IF 1323: END IF 1324: AA = DBLE( A( L+1 ) ) 1325: * L(k+i+1,k+i+1) 1326: IF( AA.NE.ZERO ) THEN 1327: IF( SCALE.LT.AA ) THEN 1328: S = ONE + S*( SCALE / AA )**2 1329: SCALE = AA 1330: ELSE 1331: S = S + ( AA / SCALE )**2 1332: END IF 1333: END IF 1334: L = L + LDA + 1 1335: END DO 1336: * L-> k - 1 + k*lda or L(k-1,k-1) at A(k-1,k) 1337: AA = DBLE( A( L ) ) 1338: * L(k-1,k-1) at A(k-1,k) 1339: IF( AA.NE.ZERO ) THEN 1340: IF( SCALE.LT.AA ) THEN 1341: S = ONE + S*( SCALE / AA )**2 1342: SCALE = AA 1343: ELSE 1344: S = S + ( AA / SCALE )**2 1345: END IF 1346: END IF 1347: END IF 1348: END IF 1349: END IF 1350: VALUE = SCALE*SQRT( S ) 1351: END IF 1352: * 1353: ZLANHF = VALUE 1354: RETURN 1355: * 1356: * End of ZLANHF 1357: * 1358: END