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