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