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