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