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