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