LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
clanhf.f
Go to the documentation of this file.
1 *> \brief \b CLANHF
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CLANHF + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clanhf.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clanhf.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clanhf.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * REAL FUNCTION CLANHF( NORM, TRANSR, UPLO, N, A, WORK )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER NORM, TRANSR, UPLO
25 * INTEGER N
26 * ..
27 * .. Array Arguments ..
28 * REAL WORK( 0: * )
29 * COMPLEX A( 0: * )
30 * ..
31 *
32 *
33 *> \par Purpose:
34 * =============
35 *>
36 *> \verbatim
37 *>
38 *> CLANHF 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 CLANHF
44 *> \verbatim
45 *>
46 *> CLANHF = ( 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 CLANHF 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, CLANHF 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 REAL 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 complexOTHERcomputational
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  REAL FUNCTION clanhf( 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  REAL work( 0: * )
260  COMPLEX a( 0: * )
261 * ..
262 *
263 * =====================================================================
264 *
265 * .. Parameters ..
266  REAL one, zero
267  parameter( one = 1.0e+0, zero = 0.0e+0 )
268 * ..
269 * .. Local Scalars ..
270  INTEGER i, j, ifm, ilu, noe, n1, k, l, lda
271  REAL scale, s, value, aa
272 * ..
273 * .. External Functions ..
274  LOGICAL lsame
275  INTEGER isamax
276  EXTERNAL lsame, isamax
277 * ..
278 * .. External Subroutines ..
279  EXTERNAL classq
280 * ..
281 * .. Intrinsic Functions ..
282  INTRINSIC abs, REAL, max, sqrt
283 * ..
284 * .. Executable Statements ..
285 *
286  IF( n.EQ.0 ) THEN
287  clanhf = zero
288  RETURN
289  ELSE IF( n.EQ.1 ) THEN
290  clanhf = 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( REAL( 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( REAL( A( I+J*LDA ) ) ) )
353  i = j
354 * -> L(j,j)
355  value = max( value, abs( REAL( 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( REAL( A( I+J*LDA ) ) ) )
369  i = i + 1
370 * =k+j; i -> U(j,j)
371  value = max( value, abs( REAL( 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( REAL( 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( REAL( A( I+J*LDA ) ) ) )
394  i = j + 1
395 * L(j+k,j+k)
396  value = max( value, abs( REAL( 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( REAL( 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( REAL( 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( REAL( A( I+J*LDA ) ) ) )
433  i = j - k + 1
434 * U(j,j)
435  value = max( value, abs( REAL( 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( REAL( A( J+J*LDA ) ) ) )
451  value = max( value, abs( REAL( 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( REAL( A( I+J*LDA ) ) ) )
462  i = j + 1
463 * -> L(j,j)
464  value = max( value, abs( REAL( 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( REAL( A( I+J*LDA ) ) ) )
478  i = i + 1
479 * =k+j+1; i -> U(j,j)
480  value = max( value, abs( REAL( 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( REAL( A( I+J*LDA ) ) ) )
491  i = n
492 * -> U(k-1,k-1)
493  value = max( value, abs( REAL( 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( REAL( 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( REAL( A( I+J*LDA ) ) ) )
512  i = j
513 * L(j+k,j+k)
514  value = max( value, abs( REAL( 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( REAL( 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( REAL( 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( REAL( A( I+J*LDA ) ) ) )
551  i = j - k
552 * U(j,j)
553  value = max( value, abs( REAL( 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( REAL( 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( REAL( 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( REAL( 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 = isamax( 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( REAL( 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( REAL( 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 = isamax( 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( REAL( A( I+J*LDA ) ) )
668 * -> A(j+k,j+k)
669  work( j+k ) = s + aa
670  i = i + 1
671  aa = abs( REAL( 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 = isamax( 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( REAL( 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( REAL( 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 = isamax( 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( REAL( 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( REAL( 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( REAL( 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 = isamax( 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( REAL( 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( REAL( 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( REAL( 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 = isamax( 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( REAL( 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( REAL( 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( REAL( 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( REAL( A( I+J*LDA ) ) )
909 * A(k-1,k-1)
910  s = s + aa
911  work( i ) = work( i ) + s
912  i = isamax( 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( REAL( 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( REAL( 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( REAL( 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( REAL( 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 = isamax( 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 classq( 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 classq( 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 = REAL( 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 = REAL( 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 = REAL( 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 classq( 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 classq( 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 = REAL( 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 = REAL( 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 = REAL( 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 classq( 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 classq( 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 classq( 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 = REAL( 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 = REAL( 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 = REAL( 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 classq( 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 classq( 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 classq( 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 = REAL( 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 = REAL( 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 = REAL( 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 classq( 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 classq( 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 = REAL( 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 = REAL( 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 classq( 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 classq( 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 = REAL( 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 = REAL( 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 classq( 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 classq( 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 classq( 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 = REAL( 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 = REAL( 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 = REAL( 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 = REAL( 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 classq( 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 classq( 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 classq( 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 = REAL( 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 = REAL( 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 = REAL( 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 = REAL( 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  clanhf = value
1422  RETURN
1423 *
1424 * End of CLANHF
1425 *
1426  END