LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 *> \htmlonly
9 *> Download DLANSF + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlansf.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlansf.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlansf.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * DOUBLE PRECISION FUNCTION DLANSF( NORM, TRANSR, UPLO, N, A, WORK )
22 *
23 * .. Scalar Arguments ..
24 * CHARACTER NORM, TRANSR, UPLO
25 * INTEGER N
26 * ..
27 * .. Array Arguments ..
28 * DOUBLE PRECISION A( 0: * ), WORK( 0: * )
29 * ..
30 *
31 *
32 *> \par Purpose:
33 * =============
34 *>
35 *> \verbatim
36 *>
37 *> DLANSF returns the value of the one norm, or the Frobenius norm, or
38 *> the infinity norm, or the element of largest absolute value of a
39 *> real symmetric matrix A in RFP format.
40 *> \endverbatim
41 *>
42 *> \return DLANSF
43 *> \verbatim
44 *>
45 *> DLANSF = ( max(abs(A(i,j))), NORM = 'M' or 'm'
46 *> (
47 *> ( norm1(A), NORM = '1', 'O' or 'o'
48 *> (
49 *> ( normI(A), NORM = 'I' or 'i'
50 *> (
51 *> ( normF(A), NORM = 'F', 'f', 'E' or 'e'
52 *>
53 *> where norm1 denotes the one norm of a matrix (maximum column sum),
54 *> normI denotes the infinity norm of a matrix (maximum row sum) and
55 *> normF denotes the Frobenius norm of a matrix (square root of sum of
56 *> squares). Note that max(abs(A(i,j))) is not a matrix norm.
57 *> \endverbatim
58 *
59 * Arguments:
60 * ==========
61 *
62 *> \param[in] NORM
63 *> \verbatim
64 *> NORM is CHARACTER*1
65 *> Specifies the value to be returned in DLANSF as described
66 *> above.
67 *> \endverbatim
68 *>
69 *> \param[in] TRANSR
70 *> \verbatim
71 *> TRANSR is CHARACTER*1
72 *> Specifies whether the RFP format of A is normal or
73 *> transposed format.
74 *> = 'N': RFP format is Normal;
75 *> = 'T': RFP format is Transpose.
76 *> \endverbatim
77 *>
78 *> \param[in] UPLO
79 *> \verbatim
80 *> UPLO is CHARACTER*1
81 *> On entry, UPLO specifies whether the RFP matrix A came from
82 *> an upper or lower triangular matrix as follows:
83 *> = 'U': RFP A came from an upper triangular matrix;
84 *> = 'L': RFP A came from a lower triangular matrix.
85 *> \endverbatim
86 *>
87 *> \param[in] N
88 *> \verbatim
89 *> N is INTEGER
90 *> The order of the matrix A. N >= 0. When N = 0, DLANSF is
91 *> set to zero.
92 *> \endverbatim
93 *>
94 *> \param[in] A
95 *> \verbatim
96 *> A is DOUBLE PRECISION array, dimension ( N*(N+1)/2 );
97 *> On entry, the upper (if UPLO = 'U') or lower (if UPLO = 'L')
98 *> part of the symmetric matrix A stored in RFP format. See the
99 *> "Notes" below for more details.
100 *> Unchanged on exit.
101 *> \endverbatim
102 *>
103 *> \param[out] WORK
104 *> \verbatim
105 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
106 *> where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
107 *> WORK is not referenced.
108 *> \endverbatim
109 *
110 * Authors:
111 * ========
112 *
113 *> \author Univ. of Tennessee
114 *> \author Univ. of California Berkeley
115 *> \author Univ. of Colorado Denver
116 *> \author NAG Ltd.
117 *
118 *> \ingroup doubleOTHERcomputational
119 *
120 *> \par Further Details:
121 * =====================
122 *>
123 *> \verbatim
124 *>
125 *> We first consider Rectangular Full Packed (RFP) Format when N is
126 *> even. We give an example where N = 6.
127 *>
128 *> AP is Upper AP is Lower
129 *>
130 *> 00 01 02 03 04 05 00
131 *> 11 12 13 14 15 10 11
132 *> 22 23 24 25 20 21 22
133 *> 33 34 35 30 31 32 33
134 *> 44 45 40 41 42 43 44
135 *> 55 50 51 52 53 54 55
136 *>
137 *>
138 *> Let TRANSR = 'N'. RFP holds AP as follows:
139 *> For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
140 *> three columns of AP upper. The lower triangle A(4:6,0:2) consists of
141 *> the transpose of the first three columns of AP upper.
142 *> For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
143 *> three columns of AP lower. The upper triangle A(0:2,0:2) consists of
144 *> the transpose of the last three columns of AP lower.
145 *> This covers the case N even and TRANSR = 'N'.
146 *>
147 *> RFP A RFP A
148 *>
149 *> 03 04 05 33 43 53
150 *> 13 14 15 00 44 54
151 *> 23 24 25 10 11 55
152 *> 33 34 35 20 21 22
153 *> 00 44 45 30 31 32
154 *> 01 11 55 40 41 42
155 *> 02 12 22 50 51 52
156 *>
157 *> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
158 *> transpose of RFP A above. One therefore gets:
159 *>
160 *>
161 *> RFP A RFP A
162 *>
163 *> 03 13 23 33 00 01 02 33 00 10 20 30 40 50
164 *> 04 14 24 34 44 11 12 43 44 11 21 31 41 51
165 *> 05 15 25 35 45 55 22 53 54 55 22 32 42 52
166 *>
167 *>
168 *> We then consider Rectangular Full Packed (RFP) Format when N is
169 *> odd. We give an example where N = 5.
170 *>
171 *> AP is Upper AP is Lower
172 *>
173 *> 00 01 02 03 04 00
174 *> 11 12 13 14 10 11
175 *> 22 23 24 20 21 22
176 *> 33 34 30 31 32 33
177 *> 44 40 41 42 43 44
178 *>
179 *>
180 *> Let TRANSR = 'N'. RFP holds AP as follows:
181 *> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
182 *> three columns of AP upper. The lower triangle A(3:4,0:1) consists of
183 *> the transpose of the first two columns of AP upper.
184 *> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
185 *> three columns of AP lower. The upper triangle A(0:1,1:2) consists of
186 *> the transpose of the last two columns of AP lower.
187 *> This covers the case N odd and TRANSR = 'N'.
188 *>
189 *> RFP A RFP A
190 *>
191 *> 02 03 04 00 33 43
192 *> 12 13 14 10 11 44
193 *> 22 23 24 20 21 22
194 *> 00 33 34 30 31 32
195 *> 01 11 44 40 41 42
196 *>
197 *> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
198 *> transpose of RFP A above. One therefore gets:
199 *>
200 *> RFP A RFP A
201 *>
202 *> 02 12 22 00 01 00 10 20 30 40 50
203 *> 03 13 23 33 11 33 11 21 31 41 51
204 *> 04 14 24 34 44 43 44 22 32 42 52
205 *> \endverbatim
206 *
207 * =====================================================================
208  DOUBLE PRECISION FUNCTION dlansf( NORM, TRANSR, UPLO, N, A, WORK )
209 *
210 * -- LAPACK computational routine --
211 * -- LAPACK is a software package provided by Univ. of Tennessee, --
212 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
213 *
214 * .. Scalar Arguments ..
215  CHARACTER norm, transr, uplo
216  INTEGER n
217 * ..
218 * .. Array Arguments ..
219  DOUBLE PRECISION a( 0: * ), work( 0: * )
220 * ..
221 *
222 * =====================================================================
223 *
224 * .. Parameters ..
225  DOUBLE PRECISION one, zero
226  parameter( one = 1.0d+0, zero = 0.0d+0 )
227 * ..
228 * .. Local Scalars ..
229  INTEGER i, j, ifm, ilu, noe, n1, k, l, lda
230  DOUBLE PRECISION scale, s, VALUE, aa, temp
231 * ..
232 * .. External Functions ..
233  LOGICAL lsame, disnan
234  EXTERNAL lsame, disnan
235 * ..
236 * .. External Subroutines ..
237  EXTERNAL dlassq
238 * ..
239 * .. Intrinsic Functions ..
240  INTRINSIC abs, max, sqrt
241 * ..
242 * .. Executable Statements ..
243 *
244  IF( n.EQ.0 ) THEN
245  dlansf = zero
246  RETURN
247  ELSE IF( n.EQ.1 ) THEN
248  dlansf = abs( a(0) )
249  RETURN
250  END IF
251 *
252 * set noe = 1 if n is odd. if n is even set noe=0
253 *
254  noe = 1
255  IF( mod( n, 2 ).EQ.0 )
256  $ noe = 0
257 *
258 * set ifm = 0 when form='T or 't' and 1 otherwise
259 *
260  ifm = 1
261  IF( lsame( transr, 'T' ) )
262  $ ifm = 0
263 *
264 * set ilu = 0 when uplo='U or 'u' and 1 otherwise
265 *
266  ilu = 1
267  IF( lsame( uplo, 'U' ) )
268  $ ilu = 0
269 *
270 * set lda = (n+1)/2 when ifm = 0
271 * set lda = n when ifm = 1 and noe = 1
272 * set lda = n+1 when ifm = 1 and noe = 0
273 *
274  IF( ifm.EQ.1 ) THEN
275  IF( noe.EQ.1 ) THEN
276  lda = n
277  ELSE
278 * noe=0
279  lda = n + 1
280  END IF
281  ELSE
282 * ifm=0
283  lda = ( n+1 ) / 2
284  END IF
285 *
286  IF( lsame( norm, 'M' ) ) THEN
287 *
288 * Find max(abs(A(i,j))).
289 *
290  k = ( n+1 ) / 2
291  VALUE = zero
292  IF( noe.EQ.1 ) THEN
293 * n is odd
294  IF( ifm.EQ.1 ) THEN
295 * A is n by k
296  DO j = 0, k - 1
297  DO i = 0, n - 1
298  temp = abs( a( i+j*lda ) )
299  IF( VALUE .LT. temp .OR. disnan( temp ) )
300  $ VALUE = temp
301  END DO
302  END DO
303  ELSE
304 * xpose case; A is k by n
305  DO j = 0, n - 1
306  DO i = 0, k - 1
307  temp = abs( a( i+j*lda ) )
308  IF( VALUE .LT. temp .OR. disnan( temp ) )
309  $ VALUE = temp
310  END DO
311  END DO
312  END IF
313  ELSE
314 * n is even
315  IF( ifm.EQ.1 ) THEN
316 * A is n+1 by k
317  DO j = 0, k - 1
318  DO i = 0, n
319  temp = abs( a( i+j*lda ) )
320  IF( VALUE .LT. temp .OR. disnan( temp ) )
321  $ VALUE = temp
322  END DO
323  END DO
324  ELSE
325 * xpose case; A is k by n+1
326  DO j = 0, n
327  DO i = 0, k - 1
328  temp = abs( a( i+j*lda ) )
329  IF( VALUE .LT. temp .OR. disnan( temp ) )
330  $ VALUE = temp
331  END DO
332  END DO
333  END IF
334  END IF
335  ELSE IF( ( lsame( norm, 'I' ) ) .OR. ( 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. ( lsame( norm, 'E' ) ) ) THEN
781 *
782 * Find normF(A).
783 *
784  k = ( n+1 ) / 2
785  scale = zero
786  s = one
787  IF( noe.EQ.1 ) THEN
788 * n is odd
789  IF( ifm.EQ.1 ) THEN
790 * A is normal
791  IF( ilu.EQ.0 ) THEN
792 * A is upper
793  DO j = 0, k - 3
794  CALL dlassq( k-j-2, a( k+j+1+j*lda ), 1, scale, s )
795 * L at A(k,0)
796  END DO
797  DO j = 0, k - 1
798  CALL dlassq( k+j-1, a( 0+j*lda ), 1, scale, s )
799 * trap U at A(0,0)
800  END DO
801  s = s + s
802 * double s for the off diagonal elements
803  CALL dlassq( k-1, a( k ), lda+1, scale, s )
804 * tri L at A(k,0)
805  CALL dlassq( k, a( k-1 ), lda+1, scale, s )
806 * tri U at A(k-1,0)
807  ELSE
808 * ilu=1 & A is lower
809  DO j = 0, k - 1
810  CALL dlassq( n-j-1, a( j+1+j*lda ), 1, scale, s )
811 * trap L at A(0,0)
812  END DO
813  DO j = 0, k - 2
814  CALL dlassq( j, a( 0+( 1+j )*lda ), 1, scale, s )
815 * U at A(0,1)
816  END DO
817  s = s + s
818 * double s for the off diagonal elements
819  CALL dlassq( k, a( 0 ), lda+1, scale, s )
820 * tri L at A(0,0)
821  CALL dlassq( k-1, a( 0+lda ), lda+1, scale, s )
822 * tri U at A(0,1)
823  END IF
824  ELSE
825 * A is xpose
826  IF( ilu.EQ.0 ) THEN
827 * A**T is upper
828  DO j = 1, k - 2
829  CALL dlassq( j, a( 0+( k+j )*lda ), 1, scale, s )
830 * U at A(0,k)
831  END DO
832  DO j = 0, k - 2
833  CALL dlassq( k, a( 0+j*lda ), 1, scale, s )
834 * k by k-1 rect. at A(0,0)
835  END DO
836  DO j = 0, k - 2
837  CALL dlassq( k-j-1, a( j+1+( j+k-1 )*lda ), 1,
838  $ scale, s )
839 * L at A(0,k-1)
840  END DO
841  s = s + s
842 * double s for the off diagonal elements
843  CALL dlassq( k-1, a( 0+k*lda ), lda+1, scale, s )
844 * tri U at A(0,k)
845  CALL dlassq( k, a( 0+( k-1 )*lda ), lda+1, scale, s )
846 * tri L at A(0,k-1)
847  ELSE
848 * A**T is lower
849  DO j = 1, k - 1
850  CALL dlassq( j, a( 0+j*lda ), 1, scale, s )
851 * U at A(0,0)
852  END DO
853  DO j = k, n - 1
854  CALL dlassq( k, a( 0+j*lda ), 1, scale, s )
855 * k by k-1 rect. at A(0,k)
856  END DO
857  DO j = 0, k - 3
858  CALL dlassq( k-j-2, a( j+2+j*lda ), 1, scale, s )
859 * L at A(1,0)
860  END DO
861  s = s + s
862 * double s for the off diagonal elements
863  CALL dlassq( k, a( 0 ), lda+1, scale, s )
864 * tri U at A(0,0)
865  CALL dlassq( k-1, a( 1 ), lda+1, scale, s )
866 * tri L at A(1,0)
867  END IF
868  END IF
869  ELSE
870 * n is even
871  IF( ifm.EQ.1 ) THEN
872 * A is normal
873  IF( ilu.EQ.0 ) THEN
874 * A is upper
875  DO j = 0, k - 2
876  CALL dlassq( k-j-1, a( k+j+2+j*lda ), 1, scale, s )
877 * L at A(k+1,0)
878  END DO
879  DO j = 0, k - 1
880  CALL dlassq( k+j, a( 0+j*lda ), 1, scale, s )
881 * trap U at A(0,0)
882  END DO
883  s = s + s
884 * double s for the off diagonal elements
885  CALL dlassq( k, a( k+1 ), lda+1, scale, s )
886 * tri L at A(k+1,0)
887  CALL dlassq( k, a( k ), lda+1, scale, s )
888 * tri U at A(k,0)
889  ELSE
890 * ilu=1 & A is lower
891  DO j = 0, k - 1
892  CALL dlassq( n-j-1, a( j+2+j*lda ), 1, scale, s )
893 * trap L at A(1,0)
894  END DO
895  DO j = 1, k - 1
896  CALL dlassq( j, a( 0+j*lda ), 1, scale, s )
897 * U at A(0,0)
898  END DO
899  s = s + s
900 * double s for the off diagonal elements
901  CALL dlassq( k, a( 1 ), lda+1, scale, s )
902 * tri L at A(1,0)
903  CALL dlassq( k, a( 0 ), lda+1, scale, s )
904 * tri U at A(0,0)
905  END IF
906  ELSE
907 * A is xpose
908  IF( ilu.EQ.0 ) THEN
909 * A**T is upper
910  DO j = 1, k - 1
911  CALL dlassq( j, a( 0+( k+1+j )*lda ), 1, scale, s )
912 * U at A(0,k+1)
913  END DO
914  DO j = 0, k - 1
915  CALL dlassq( k, a( 0+j*lda ), 1, scale, s )
916 * k by k rect. at A(0,0)
917  END DO
918  DO j = 0, k - 2
919  CALL dlassq( k-j-1, a( j+1+( j+k )*lda ), 1, scale,
920  $ s )
921 * L at A(0,k)
922  END DO
923  s = s + s
924 * double s for the off diagonal elements
925  CALL dlassq( k, a( 0+( k+1 )*lda ), lda+1, scale, s )
926 * tri U at A(0,k+1)
927  CALL dlassq( k, a( 0+k*lda ), lda+1, scale, s )
928 * tri L at A(0,k)
929  ELSE
930 * A**T is lower
931  DO j = 1, k - 1
932  CALL dlassq( j, a( 0+( j+1 )*lda ), 1, scale, s )
933 * U at A(0,1)
934  END DO
935  DO j = k + 1, n
936  CALL dlassq( k, a( 0+j*lda ), 1, scale, s )
937 * k by k rect. at A(0,k+1)
938  END DO
939  DO j = 0, k - 2
940  CALL dlassq( k-j-1, a( j+1+j*lda ), 1, scale, s )
941 * L at A(0,0)
942  END DO
943  s = s + s
944 * double s for the off diagonal elements
945  CALL dlassq( k, a( lda ), lda+1, scale, s )
946 * tri L at A(0,1)
947  CALL dlassq( k, a( 0 ), lda+1, scale, s )
948 * tri U at A(0,0)
949  END IF
950  END IF
951  END IF
952  VALUE = scale*sqrt( s )
953  END IF
954 *
955  dlansf = VALUE
956  RETURN
957 *
958 * End of DLANSF
959 *
960  END
logical function disnan(DIN)
DISNAN tests input for NaN.
Definition: disnan.f:59
subroutine dlassq(n, x, incx, scl, sumsq)
DLASSQ updates a sum of squares represented in scaled form.
Definition: dlassq.f90:126
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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:209