LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlansf.f
Go to the documentation of this file.
1*> \brief \b DLANSF returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric 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 DLANSF + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlansf.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlansf.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlansf.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* DOUBLE PRECISION FUNCTION DLANSF( NORM, TRANSR, UPLO, N, A, WORK )
20*
21* .. Scalar Arguments ..
22* CHARACTER NORM, TRANSR, UPLO
23* INTEGER N
24* ..
25* .. Array Arguments ..
26* DOUBLE PRECISION A( 0: * ), WORK( 0: * )
27* ..
28*
29*
30*> \par Purpose:
31* =============
32*>
33*> \verbatim
34*>
35*> DLANSF returns the value of the one norm, or the Frobenius norm, or
36*> the infinity norm, or the element of largest absolute value of a
37*> real symmetric matrix A in RFP format.
38*> \endverbatim
39*>
40*> \return DLANSF
41*> \verbatim
42*>
43*> DLANSF = ( max(abs(A(i,j))), NORM = 'M' or 'm'
44*> (
45*> ( norm1(A), NORM = '1', 'O' or 'o'
46*> (
47*> ( normI(A), NORM = 'I' or 'i'
48*> (
49*> ( normF(A), NORM = 'F', 'f', 'E' or 'e'
50*>
51*> where norm1 denotes the one norm of a matrix (maximum column sum),
52*> normI denotes the infinity norm of a matrix (maximum row sum) and
53*> normF denotes the Frobenius norm of a matrix (square root of sum of
54*> squares). Note that max(abs(A(i,j))) is not a matrix norm.
55*> \endverbatim
56*
57* Arguments:
58* ==========
59*
60*> \param[in] NORM
61*> \verbatim
62*> NORM is CHARACTER*1
63*> Specifies the value to be returned in DLANSF as described
64*> above.
65*> \endverbatim
66*>
67*> \param[in] TRANSR
68*> \verbatim
69*> TRANSR is CHARACTER*1
70*> Specifies whether the RFP format of A is normal or
71*> transposed format.
72*> = 'N': RFP format is Normal;
73*> = 'T': RFP format is Transpose.
74*> \endverbatim
75*>
76*> \param[in] UPLO
77*> \verbatim
78*> UPLO is CHARACTER*1
79*> On entry, UPLO specifies whether the RFP matrix A came from
80*> an upper or lower triangular matrix as follows:
81*> = 'U': RFP A came from an upper triangular matrix;
82*> = 'L': RFP A came from a lower triangular matrix.
83*> \endverbatim
84*>
85*> \param[in] N
86*> \verbatim
87*> N is INTEGER
88*> The order of the matrix A. N >= 0. When N = 0, DLANSF is
89*> set to zero.
90*> \endverbatim
91*>
92*> \param[in] A
93*> \verbatim
94*> A is DOUBLE PRECISION array, dimension ( N*(N+1)/2 );
95*> On entry, the upper (if UPLO = 'U') or lower (if UPLO = 'L')
96*> part of the symmetric matrix A stored in RFP format. See the
97*> "Notes" below for more details.
98*> Unchanged on exit.
99*> \endverbatim
100*>
101*> \param[out] WORK
102*> \verbatim
103*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
104*> where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
105*> WORK is not referenced.
106*> \endverbatim
107*
108* Authors:
109* ========
110*
111*> \author Univ. of Tennessee
112*> \author Univ. of California Berkeley
113*> \author Univ. of Colorado Denver
114*> \author NAG Ltd.
115*
116*> \ingroup lanhf
117*
118*> \par Further Details:
119* =====================
120*>
121*> \verbatim
122*>
123*> We first consider Rectangular Full Packed (RFP) Format when N is
124*> even. We give an example where N = 6.
125*>
126*> AP is Upper AP is Lower
127*>
128*> 00 01 02 03 04 05 00
129*> 11 12 13 14 15 10 11
130*> 22 23 24 25 20 21 22
131*> 33 34 35 30 31 32 33
132*> 44 45 40 41 42 43 44
133*> 55 50 51 52 53 54 55
134*>
135*>
136*> Let TRANSR = 'N'. RFP holds AP as follows:
137*> For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
138*> three columns of AP upper. The lower triangle A(4:6,0:2) consists of
139*> the transpose of the first three columns of AP upper.
140*> For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
141*> three columns of AP lower. The upper triangle A(0:2,0:2) consists of
142*> the transpose of the last three columns of AP lower.
143*> This covers the case N even and TRANSR = 'N'.
144*>
145*> RFP A RFP A
146*>
147*> 03 04 05 33 43 53
148*> 13 14 15 00 44 54
149*> 23 24 25 10 11 55
150*> 33 34 35 20 21 22
151*> 00 44 45 30 31 32
152*> 01 11 55 40 41 42
153*> 02 12 22 50 51 52
154*>
155*> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
156*> transpose of RFP A above. One therefore gets:
157*>
158*>
159*> RFP A RFP A
160*>
161*> 03 13 23 33 00 01 02 33 00 10 20 30 40 50
162*> 04 14 24 34 44 11 12 43 44 11 21 31 41 51
163*> 05 15 25 35 45 55 22 53 54 55 22 32 42 52
164*>
165*>
166*> We then consider Rectangular Full Packed (RFP) Format when N is
167*> odd. We give an example where N = 5.
168*>
169*> AP is Upper AP is Lower
170*>
171*> 00 01 02 03 04 00
172*> 11 12 13 14 10 11
173*> 22 23 24 20 21 22
174*> 33 34 30 31 32 33
175*> 44 40 41 42 43 44
176*>
177*>
178*> Let TRANSR = 'N'. RFP holds AP as follows:
179*> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
180*> three columns of AP upper. The lower triangle A(3:4,0:1) consists of
181*> the transpose of the first two columns of AP upper.
182*> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
183*> three columns of AP lower. The upper triangle A(0:1,1:2) consists of
184*> the transpose of the last two columns of AP lower.
185*> This covers the case N odd and TRANSR = 'N'.
186*>
187*> RFP A RFP A
188*>
189*> 02 03 04 00 33 43
190*> 12 13 14 10 11 44
191*> 22 23 24 20 21 22
192*> 00 33 34 30 31 32
193*> 01 11 44 40 41 42
194*>
195*> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
196*> transpose of RFP A above. One therefore gets:
197*>
198*> RFP A RFP A
199*>
200*> 02 12 22 00 01 00 10 20 30 40 50
201*> 03 13 23 33 11 33 11 21 31 41 51
202*> 04 14 24 34 44 43 44 22 32 42 52
203*> \endverbatim
204*
205* =====================================================================
206 DOUBLE PRECISION FUNCTION dlansf( NORM, TRANSR, UPLO, N, A,
207 $ WORK )
208*
209* -- LAPACK computational routine --
210* -- LAPACK is a software package provided by Univ. of Tennessee, --
211* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
212*
213* .. Scalar Arguments ..
214 CHARACTER norm, transr, uplo
215 INTEGER n
216* ..
217* .. Array Arguments ..
218 DOUBLE PRECISION a( 0: * ), work( 0: * )
219* ..
220*
221* =====================================================================
222*
223* .. Parameters ..
224 DOUBLE PRECISION one, zero
225 parameter( one = 1.0d+0, zero = 0.0d+0 )
226* ..
227* .. Local Scalars ..
228 INTEGER i, j, ifm, ilu, noe, n1, k, l, lda
229 DOUBLE PRECISION scale, s, VALUE, aa, temp
230* ..
231* .. External Functions ..
232 LOGICAL lsame, disnan
233 EXTERNAL lsame, disnan
234* ..
235* .. External Subroutines ..
236 EXTERNAL dlassq
237* ..
238* .. Intrinsic Functions ..
239 INTRINSIC abs, max, sqrt
240* ..
241* .. Executable Statements ..
242*
243 IF( n.EQ.0 ) THEN
244 dlansf = zero
245 RETURN
246 ELSE IF( n.EQ.1 ) THEN
247 dlansf = abs( a(0) )
248 RETURN
249 END IF
250*
251* set noe = 1 if n is odd. if n is even set noe=0
252*
253 noe = 1
254 IF( mod( n, 2 ).EQ.0 )
255 $ noe = 0
256*
257* set ifm = 0 when form='T or 't' and 1 otherwise
258*
259 ifm = 1
260 IF( lsame( transr, 'T' ) )
261 $ ifm = 0
262*
263* set ilu = 0 when uplo='U or 'u' and 1 otherwise
264*
265 ilu = 1
266 IF( lsame( uplo, 'U' ) )
267 $ ilu = 0
268*
269* set lda = (n+1)/2 when ifm = 0
270* set lda = n when ifm = 1 and noe = 1
271* set lda = n+1 when ifm = 1 and noe = 0
272*
273 IF( ifm.EQ.1 ) THEN
274 IF( noe.EQ.1 ) THEN
275 lda = n
276 ELSE
277* noe=0
278 lda = n + 1
279 END IF
280 ELSE
281* ifm=0
282 lda = ( n+1 ) / 2
283 END IF
284*
285 IF( lsame( norm, 'M' ) ) THEN
286*
287* Find max(abs(A(i,j))).
288*
289 k = ( n+1 ) / 2
290 VALUE = zero
291 IF( noe.EQ.1 ) THEN
292* n is odd
293 IF( ifm.EQ.1 ) THEN
294* A is n by k
295 DO j = 0, k - 1
296 DO i = 0, n - 1
297 temp = abs( a( i+j*lda ) )
298 IF( VALUE .LT. temp .OR. disnan( temp ) )
299 $ VALUE = temp
300 END DO
301 END DO
302 ELSE
303* xpose case; A is k by n
304 DO j = 0, n - 1
305 DO i = 0, k - 1
306 temp = abs( a( i+j*lda ) )
307 IF( VALUE .LT. temp .OR. disnan( temp ) )
308 $ VALUE = temp
309 END DO
310 END DO
311 END IF
312 ELSE
313* n is even
314 IF( ifm.EQ.1 ) THEN
315* A is n+1 by k
316 DO j = 0, k - 1
317 DO i = 0, n
318 temp = abs( a( i+j*lda ) )
319 IF( VALUE .LT. temp .OR. disnan( temp ) )
320 $ VALUE = temp
321 END DO
322 END DO
323 ELSE
324* xpose case; A is k by n+1
325 DO j = 0, n
326 DO i = 0, k - 1
327 temp = abs( a( i+j*lda ) )
328 IF( VALUE .LT. temp .OR. disnan( temp ) )
329 $ VALUE = temp
330 END DO
331 END DO
332 END IF
333 END IF
334 ELSE IF( ( lsame( norm, 'I' ) ) .OR.
335 $ ( lsame( norm, 'O' ) ) .OR.
336 $ ( norm.EQ.'1' ) ) THEN
337*
338* Find normI(A) ( = norm1(A), since A is symmetric).
339*
340 IF( ifm.EQ.1 ) THEN
341 k = n / 2
342 IF( noe.EQ.1 ) THEN
343* n is odd
344 IF( ilu.EQ.0 ) THEN
345 DO i = 0, k - 1
346 work( i ) = zero
347 END DO
348 DO j = 0, k
349 s = zero
350 DO i = 0, k + j - 1
351 aa = abs( a( i+j*lda ) )
352* -> A(i,j+k)
353 s = s + aa
354 work( i ) = work( i ) + aa
355 END DO
356 aa = abs( a( i+j*lda ) )
357* -> A(j+k,j+k)
358 work( j+k ) = s + aa
359 IF( i.EQ.k+k )
360 $ GO TO 10
361 i = i + 1
362 aa = abs( a( i+j*lda ) )
363* -> A(j,j)
364 work( j ) = work( j ) + aa
365 s = zero
366 DO l = j + 1, k - 1
367 i = i + 1
368 aa = abs( a( i+j*lda ) )
369* -> A(l,j)
370 s = s + aa
371 work( l ) = work( l ) + aa
372 END DO
373 work( j ) = work( j ) + s
374 END DO
375 10 CONTINUE
376 VALUE = work( 0 )
377 DO i = 1, n-1
378 temp = work( i )
379 IF( VALUE .LT. temp .OR. disnan( temp ) )
380 $ VALUE = temp
381 END DO
382 ELSE
383* ilu = 1
384 k = k + 1
385* k=(n+1)/2 for n odd and ilu=1
386 DO i = k, n - 1
387 work( i ) = zero
388 END DO
389 DO j = k - 1, 0, -1
390 s = zero
391 DO i = 0, j - 2
392 aa = abs( a( i+j*lda ) )
393* -> A(j+k,i+k)
394 s = s + aa
395 work( i+k ) = work( i+k ) + aa
396 END DO
397 IF( j.GT.0 ) THEN
398 aa = abs( a( i+j*lda ) )
399* -> A(j+k,j+k)
400 s = s + aa
401 work( i+k ) = work( i+k ) + s
402* i=j
403 i = i + 1
404 END IF
405 aa = abs( a( i+j*lda ) )
406* -> A(j,j)
407 work( j ) = aa
408 s = zero
409 DO l = j + 1, n - 1
410 i = i + 1
411 aa = abs( a( i+j*lda ) )
412* -> A(l,j)
413 s = s + aa
414 work( l ) = work( l ) + aa
415 END DO
416 work( j ) = work( j ) + s
417 END DO
418 VALUE = work( 0 )
419 DO i = 1, n-1
420 temp = work( i )
421 IF( VALUE .LT. temp .OR. disnan( temp ) )
422 $ VALUE = temp
423 END DO
424 END IF
425 ELSE
426* n is even
427 IF( ilu.EQ.0 ) THEN
428 DO i = 0, k - 1
429 work( i ) = zero
430 END DO
431 DO j = 0, k - 1
432 s = zero
433 DO i = 0, k + j - 1
434 aa = abs( a( i+j*lda ) )
435* -> A(i,j+k)
436 s = s + aa
437 work( i ) = work( i ) + aa
438 END DO
439 aa = abs( a( i+j*lda ) )
440* -> A(j+k,j+k)
441 work( j+k ) = s + aa
442 i = i + 1
443 aa = abs( a( i+j*lda ) )
444* -> A(j,j)
445 work( j ) = work( j ) + aa
446 s = zero
447 DO l = j + 1, k - 1
448 i = i + 1
449 aa = abs( a( i+j*lda ) )
450* -> A(l,j)
451 s = s + aa
452 work( l ) = work( l ) + aa
453 END DO
454 work( j ) = work( j ) + s
455 END DO
456 VALUE = work( 0 )
457 DO i = 1, n-1
458 temp = work( i )
459 IF( VALUE .LT. temp .OR. disnan( temp ) )
460 $ VALUE = temp
461 END DO
462 ELSE
463* ilu = 1
464 DO i = k, n - 1
465 work( i ) = zero
466 END DO
467 DO j = k - 1, 0, -1
468 s = zero
469 DO i = 0, j - 1
470 aa = abs( a( i+j*lda ) )
471* -> A(j+k,i+k)
472 s = s + aa
473 work( i+k ) = work( i+k ) + aa
474 END DO
475 aa = abs( a( i+j*lda ) )
476* -> A(j+k,j+k)
477 s = s + aa
478 work( i+k ) = work( i+k ) + s
479* i=j
480 i = i + 1
481 aa = abs( a( i+j*lda ) )
482* -> A(j,j)
483 work( j ) = aa
484 s = zero
485 DO l = j + 1, n - 1
486 i = i + 1
487 aa = abs( a( i+j*lda ) )
488* -> A(l,j)
489 s = s + aa
490 work( l ) = work( l ) + aa
491 END DO
492 work( j ) = work( j ) + s
493 END DO
494 VALUE = work( 0 )
495 DO i = 1, n-1
496 temp = work( i )
497 IF( VALUE .LT. temp .OR. disnan( temp ) )
498 $ VALUE = temp
499 END DO
500 END IF
501 END IF
502 ELSE
503* ifm=0
504 k = n / 2
505 IF( noe.EQ.1 ) THEN
506* n is odd
507 IF( ilu.EQ.0 ) THEN
508 n1 = k
509* n/2
510 k = k + 1
511* k is the row size and lda
512 DO i = n1, n - 1
513 work( i ) = zero
514 END DO
515 DO j = 0, n1 - 1
516 s = zero
517 DO i = 0, k - 1
518 aa = abs( a( i+j*lda ) )
519* A(j,n1+i)
520 work( i+n1 ) = work( i+n1 ) + aa
521 s = s + aa
522 END DO
523 work( j ) = s
524 END DO
525* j=n1=k-1 is special
526 s = abs( a( 0+j*lda ) )
527* A(k-1,k-1)
528 DO i = 1, k - 1
529 aa = abs( a( i+j*lda ) )
530* A(k-1,i+n1)
531 work( i+n1 ) = work( i+n1 ) + aa
532 s = s + aa
533 END DO
534 work( j ) = work( j ) + s
535 DO j = k, n - 1
536 s = zero
537 DO i = 0, j - k - 1
538 aa = abs( a( i+j*lda ) )
539* A(i,j-k)
540 work( i ) = work( i ) + aa
541 s = s + aa
542 END DO
543* i=j-k
544 aa = abs( a( i+j*lda ) )
545* A(j-k,j-k)
546 s = s + aa
547 work( j-k ) = work( j-k ) + s
548 i = i + 1
549 s = abs( a( i+j*lda ) )
550* A(j,j)
551 DO l = j + 1, n - 1
552 i = i + 1
553 aa = abs( a( i+j*lda ) )
554* A(j,l)
555 work( l ) = work( l ) + aa
556 s = s + aa
557 END DO
558 work( j ) = work( j ) + s
559 END DO
560 VALUE = work( 0 )
561 DO i = 1, n-1
562 temp = work( i )
563 IF( VALUE .LT. temp .OR. disnan( temp ) )
564 $ VALUE = temp
565 END DO
566 ELSE
567* ilu=1
568 k = k + 1
569* k=(n+1)/2 for n odd and ilu=1
570 DO i = k, n - 1
571 work( i ) = zero
572 END DO
573 DO j = 0, k - 2
574* process
575 s = zero
576 DO i = 0, j - 1
577 aa = abs( a( i+j*lda ) )
578* A(j,i)
579 work( i ) = work( i ) + aa
580 s = s + aa
581 END DO
582 aa = abs( a( i+j*lda ) )
583* i=j so process of A(j,j)
584 s = s + aa
585 work( j ) = s
586* is initialised here
587 i = i + 1
588* i=j process A(j+k,j+k)
589 aa = abs( a( i+j*lda ) )
590 s = aa
591 DO l = k + j + 1, n - 1
592 i = i + 1
593 aa = abs( a( i+j*lda ) )
594* A(l,k+j)
595 s = s + aa
596 work( l ) = work( l ) + aa
597 END DO
598 work( k+j ) = work( k+j ) + s
599 END DO
600* j=k-1 is special :process col A(k-1,0:k-1)
601 s = zero
602 DO i = 0, k - 2
603 aa = abs( a( i+j*lda ) )
604* A(k,i)
605 work( i ) = work( i ) + aa
606 s = s + aa
607 END DO
608* i=k-1
609 aa = abs( a( i+j*lda ) )
610* A(k-1,k-1)
611 s = s + aa
612 work( i ) = s
613* done with col j=k+1
614 DO j = k, n - 1
615* process col j of A = A(j,0:k-1)
616 s = zero
617 DO i = 0, k - 1
618 aa = abs( a( i+j*lda ) )
619* A(j,i)
620 work( i ) = work( i ) + aa
621 s = s + aa
622 END DO
623 work( j ) = work( j ) + s
624 END DO
625 VALUE = work( 0 )
626 DO i = 1, n-1
627 temp = work( i )
628 IF( VALUE .LT. temp .OR. disnan( temp ) )
629 $ VALUE = temp
630 END DO
631 END IF
632 ELSE
633* n is even
634 IF( ilu.EQ.0 ) THEN
635 DO i = k, n - 1
636 work( i ) = zero
637 END DO
638 DO j = 0, k - 1
639 s = zero
640 DO i = 0, k - 1
641 aa = abs( a( i+j*lda ) )
642* A(j,i+k)
643 work( i+k ) = work( i+k ) + aa
644 s = s + aa
645 END DO
646 work( j ) = s
647 END DO
648* j=k
649 aa = abs( a( 0+j*lda ) )
650* A(k,k)
651 s = aa
652 DO i = 1, k - 1
653 aa = abs( a( i+j*lda ) )
654* A(k,k+i)
655 work( i+k ) = work( i+k ) + aa
656 s = s + aa
657 END DO
658 work( j ) = work( j ) + s
659 DO j = k + 1, n - 1
660 s = zero
661 DO i = 0, j - 2 - k
662 aa = abs( a( i+j*lda ) )
663* A(i,j-k-1)
664 work( i ) = work( i ) + aa
665 s = s + aa
666 END DO
667* i=j-1-k
668 aa = abs( a( i+j*lda ) )
669* A(j-k-1,j-k-1)
670 s = s + aa
671 work( j-k-1 ) = work( j-k-1 ) + s
672 i = i + 1
673 aa = abs( a( i+j*lda ) )
674* A(j,j)
675 s = aa
676 DO l = j + 1, n - 1
677 i = i + 1
678 aa = abs( a( i+j*lda ) )
679* A(j,l)
680 work( l ) = work( l ) + aa
681 s = s + aa
682 END DO
683 work( j ) = work( j ) + s
684 END DO
685* j=n
686 s = zero
687 DO i = 0, k - 2
688 aa = abs( a( i+j*lda ) )
689* A(i,k-1)
690 work( i ) = work( i ) + aa
691 s = s + aa
692 END DO
693* i=k-1
694 aa = abs( a( i+j*lda ) )
695* A(k-1,k-1)
696 s = s + aa
697 work( i ) = work( i ) + s
698 VALUE = work( 0 )
699 DO i = 1, n-1
700 temp = work( i )
701 IF( VALUE .LT. temp .OR. disnan( temp ) )
702 $ VALUE = temp
703 END DO
704 ELSE
705* ilu=1
706 DO i = k, n - 1
707 work( i ) = zero
708 END DO
709* j=0 is special :process col A(k:n-1,k)
710 s = abs( a( 0 ) )
711* A(k,k)
712 DO i = 1, k - 1
713 aa = abs( a( i ) )
714* A(k+i,k)
715 work( i+k ) = work( i+k ) + aa
716 s = s + aa
717 END DO
718 work( k ) = work( k ) + s
719 DO j = 1, k - 1
720* process
721 s = zero
722 DO i = 0, j - 2
723 aa = abs( a( i+j*lda ) )
724* A(j-1,i)
725 work( i ) = work( i ) + aa
726 s = s + aa
727 END DO
728 aa = abs( a( i+j*lda ) )
729* i=j-1 so process of A(j-1,j-1)
730 s = s + aa
731 work( j-1 ) = s
732* is initialised here
733 i = i + 1
734* i=j process A(j+k,j+k)
735 aa = abs( a( i+j*lda ) )
736 s = aa
737 DO l = k + j + 1, n - 1
738 i = i + 1
739 aa = abs( a( i+j*lda ) )
740* A(l,k+j)
741 s = s + aa
742 work( l ) = work( l ) + aa
743 END DO
744 work( k+j ) = work( k+j ) + s
745 END DO
746* j=k is special :process col A(k,0:k-1)
747 s = zero
748 DO i = 0, k - 2
749 aa = abs( a( i+j*lda ) )
750* A(k,i)
751 work( i ) = work( i ) + aa
752 s = s + aa
753 END DO
754* i=k-1
755 aa = abs( a( i+j*lda ) )
756* A(k-1,k-1)
757 s = s + aa
758 work( i ) = s
759* done with col j=k+1
760 DO j = k + 1, n
761* process col j-1 of A = A(j-1,0:k-1)
762 s = zero
763 DO i = 0, k - 1
764 aa = abs( a( i+j*lda ) )
765* A(j-1,i)
766 work( i ) = work( i ) + aa
767 s = s + aa
768 END DO
769 work( j-1 ) = work( j-1 ) + s
770 END DO
771 VALUE = work( 0 )
772 DO i = 1, n-1
773 temp = work( i )
774 IF( VALUE .LT. temp .OR. disnan( temp ) )
775 $ VALUE = temp
776 END DO
777 END IF
778 END IF
779 END IF
780 ELSE IF( ( lsame( norm, 'F' ) ) .OR.
781 $ ( lsame( norm, 'E' ) ) ) THEN
782*
783* Find normF(A).
784*
785 k = ( n+1 ) / 2
786 scale = zero
787 s = one
788 IF( noe.EQ.1 ) THEN
789* n is odd
790 IF( ifm.EQ.1 ) THEN
791* A is normal
792 IF( ilu.EQ.0 ) THEN
793* A is upper
794 DO j = 0, k - 3
795 CALL dlassq( k-j-2, a( k+j+1+j*lda ), 1, scale,
796 $ s )
797* L at A(k,0)
798 END DO
799 DO j = 0, k - 1
800 CALL dlassq( k+j-1, a( 0+j*lda ), 1, scale, s )
801* trap U at A(0,0)
802 END DO
803 s = s + s
804* double s for the off diagonal elements
805 CALL dlassq( k-1, a( k ), lda+1, scale, s )
806* tri L at A(k,0)
807 CALL dlassq( k, a( k-1 ), lda+1, scale, s )
808* tri U at A(k-1,0)
809 ELSE
810* ilu=1 & A is lower
811 DO j = 0, k - 1
812 CALL dlassq( n-j-1, a( j+1+j*lda ), 1, scale,
813 $ s )
814* trap L at A(0,0)
815 END DO
816 DO j = 0, k - 2
817 CALL dlassq( j, a( 0+( 1+j )*lda ), 1, scale,
818 $ s )
819* U at A(0,1)
820 END DO
821 s = s + s
822* double s for the off diagonal elements
823 CALL dlassq( k, a( 0 ), lda+1, scale, s )
824* tri L at A(0,0)
825 CALL dlassq( k-1, a( 0+lda ), lda+1, scale, s )
826* tri U at A(0,1)
827 END IF
828 ELSE
829* A is xpose
830 IF( ilu.EQ.0 ) THEN
831* A**T is upper
832 DO j = 1, k - 2
833 CALL dlassq( j, a( 0+( k+j )*lda ), 1, scale,
834 $ s )
835* U at A(0,k)
836 END DO
837 DO j = 0, k - 2
838 CALL dlassq( k, a( 0+j*lda ), 1, scale, s )
839* k by k-1 rect. at A(0,0)
840 END DO
841 DO j = 0, k - 2
842 CALL dlassq( k-j-1, a( j+1+( j+k-1 )*lda ), 1,
843 $ scale, s )
844* L at A(0,k-1)
845 END DO
846 s = s + s
847* double s for the off diagonal elements
848 CALL dlassq( k-1, a( 0+k*lda ), lda+1, scale, s )
849* tri U at A(0,k)
850 CALL dlassq( k, a( 0+( k-1 )*lda ), lda+1, scale,
851 $ s )
852* tri L at A(0,k-1)
853 ELSE
854* A**T is lower
855 DO j = 1, k - 1
856 CALL dlassq( j, a( 0+j*lda ), 1, scale, s )
857* U at A(0,0)
858 END DO
859 DO j = k, n - 1
860 CALL dlassq( k, a( 0+j*lda ), 1, scale, s )
861* k by k-1 rect. at A(0,k)
862 END DO
863 DO j = 0, k - 3
864 CALL dlassq( k-j-2, a( j+2+j*lda ), 1, scale,
865 $ s )
866* L at A(1,0)
867 END DO
868 s = s + s
869* double s for the off diagonal elements
870 CALL dlassq( k, a( 0 ), lda+1, scale, s )
871* tri U at A(0,0)
872 CALL dlassq( k-1, a( 1 ), lda+1, scale, s )
873* tri L at A(1,0)
874 END IF
875 END IF
876 ELSE
877* n is even
878 IF( ifm.EQ.1 ) THEN
879* A is normal
880 IF( ilu.EQ.0 ) THEN
881* A is upper
882 DO j = 0, k - 2
883 CALL dlassq( k-j-1, a( k+j+2+j*lda ), 1, scale,
884 $ s )
885* L at A(k+1,0)
886 END DO
887 DO j = 0, k - 1
888 CALL dlassq( k+j, a( 0+j*lda ), 1, scale, s )
889* trap U at A(0,0)
890 END DO
891 s = s + s
892* double s for the off diagonal elements
893 CALL dlassq( k, a( k+1 ), lda+1, scale, s )
894* tri L at A(k+1,0)
895 CALL dlassq( k, a( k ), lda+1, scale, s )
896* tri U at A(k,0)
897 ELSE
898* ilu=1 & A is lower
899 DO j = 0, k - 1
900 CALL dlassq( n-j-1, a( j+2+j*lda ), 1, scale,
901 $ s )
902* trap L at A(1,0)
903 END DO
904 DO j = 1, k - 1
905 CALL dlassq( j, a( 0+j*lda ), 1, scale, s )
906* U at A(0,0)
907 END DO
908 s = s + s
909* double s for the off diagonal elements
910 CALL dlassq( k, a( 1 ), lda+1, scale, s )
911* tri L at A(1,0)
912 CALL dlassq( k, a( 0 ), lda+1, scale, s )
913* tri U at A(0,0)
914 END IF
915 ELSE
916* A is xpose
917 IF( ilu.EQ.0 ) THEN
918* A**T is upper
919 DO j = 1, k - 1
920 CALL dlassq( j, a( 0+( k+1+j )*lda ), 1, scale,
921 $ s )
922* U at A(0,k+1)
923 END DO
924 DO j = 0, k - 1
925 CALL dlassq( k, a( 0+j*lda ), 1, scale, s )
926* k by k rect. at A(0,0)
927 END DO
928 DO j = 0, k - 2
929 CALL dlassq( k-j-1, a( j+1+( j+k )*lda ), 1,
930 $ scale,
931 $ s )
932* L at A(0,k)
933 END DO
934 s = s + s
935* double s for the off diagonal elements
936 CALL dlassq( k, a( 0+( k+1 )*lda ), lda+1, scale,
937 $ s )
938* tri U at A(0,k+1)
939 CALL dlassq( k, a( 0+k*lda ), lda+1, scale, s )
940* tri L at A(0,k)
941 ELSE
942* A**T is lower
943 DO j = 1, k - 1
944 CALL dlassq( j, a( 0+( j+1 )*lda ), 1, scale,
945 $ s )
946* U at A(0,1)
947 END DO
948 DO j = k + 1, n
949 CALL dlassq( k, a( 0+j*lda ), 1, scale, s )
950* k by k rect. at A(0,k+1)
951 END DO
952 DO j = 0, k - 2
953 CALL dlassq( k-j-1, a( j+1+j*lda ), 1, scale,
954 $ s )
955* L at A(0,0)
956 END DO
957 s = s + s
958* double s for the off diagonal elements
959 CALL dlassq( k, a( lda ), lda+1, scale, s )
960* tri L at A(0,1)
961 CALL dlassq( k, a( 0 ), lda+1, scale, s )
962* tri U at A(0,0)
963 END IF
964 END IF
965 END IF
966 VALUE = scale*sqrt( s )
967 END IF
968*
969 dlansf = VALUE
970 RETURN
971*
972* End of DLANSF
973*
974 END
disnan
logical function disnan(din)
DISNAN tests input for NaN.
Definition disnan.f:57
dlansf
double precision function dlansf(norm, transr, uplo, n, a, work)
DLANSF returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
Definition dlansf.f:208
dlassq
subroutine dlassq(n, x, incx, scale, sumsq)
DLASSQ updates a sum of squares represented in scaled form.
Definition dlassq.f90:122
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48