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