![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlansf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to dlansf.f CVS log | Up to [local] / rpl / lapack / lapack |
1.11 bertrand 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.
1.7 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.15 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.7 bertrand 7: *
8: *> \htmlonly
1.15 bertrand 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">
1.7 bertrand 15: *> [TXT]</a>
1.15 bertrand 16: *> \endhtmlonly
1.7 bertrand 17: *
18: * Definition:
19: * ===========
1.1 bertrand 20: *
1.7 bertrand 21: * DOUBLE PRECISION FUNCTION DLANSF( NORM, TRANSR, UPLO, N, A, WORK )
1.15 bertrand 22: *
1.7 bertrand 23: * .. Scalar Arguments ..
24: * CHARACTER NORM, TRANSR, UPLO
25: * INTEGER N
26: * ..
27: * .. Array Arguments ..
28: * DOUBLE PRECISION A( 0: * ), WORK( 0: * )
29: * ..
1.15 bertrand 30: *
1.7 bertrand 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: *
1.15 bertrand 113: *> \author Univ. of Tennessee
114: *> \author Univ. of California Berkeley
115: *> \author Univ. of Colorado Denver
116: *> \author NAG Ltd.
1.7 bertrand 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
1.1 bertrand 206: *
1.7 bertrand 207: * =====================================================================
208: DOUBLE PRECISION FUNCTION DLANSF( NORM, TRANSR, UPLO, N, A, WORK )
1.1 bertrand 209: *
1.18 ! bertrand 210: * -- LAPACK computational routine --
1.1 bertrand 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
1.11 bertrand 230: DOUBLE PRECISION SCALE, S, VALUE, AA, TEMP
1.1 bertrand 231: * ..
232: * .. External Functions ..
1.11 bertrand 233: LOGICAL LSAME, DISNAN
234: EXTERNAL LSAME, DISNAN
1.1 bertrand 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
1.9 bertrand 247: ELSE IF( N.EQ.1 ) THEN
248: DLANSF = ABS( A(0) )
249: RETURN
1.1 bertrand 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 )
1.6 bertrand 256: $ NOE = 0
1.1 bertrand 257: *
258: * set ifm = 0 when form='T or 't' and 1 otherwise
259: *
260: IFM = 1
261: IF( LSAME( TRANSR, 'T' ) )
1.6 bertrand 262: $ IFM = 0
1.1 bertrand 263: *
264: * set ilu = 0 when uplo='U or 'u' and 1 otherwise
265: *
266: ILU = 1
267: IF( LSAME( UPLO, 'U' ) )
1.6 bertrand 268: $ ILU = 0
1.1 bertrand 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
1.11 bertrand 298: TEMP = ABS( A( I+J*LDA ) )
1.15 bertrand 299: IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
1.11 bertrand 300: $ VALUE = TEMP
1.1 bertrand 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
1.11 bertrand 307: TEMP = ABS( A( I+J*LDA ) )
1.15 bertrand 308: IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
1.11 bertrand 309: $ VALUE = TEMP
1.1 bertrand 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
1.11 bertrand 319: TEMP = ABS( A( I+J*LDA ) )
1.15 bertrand 320: IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
1.11 bertrand 321: $ VALUE = TEMP
1.1 bertrand 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
1.11 bertrand 328: TEMP = ABS( A( I+J*LDA ) )
1.15 bertrand 329: IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
1.11 bertrand 330: $ VALUE = TEMP
1.1 bertrand 331: END DO
332: END DO
333: END IF
334: END IF
335: ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR.
1.6 bertrand 336: $ ( NORM.EQ.'1' ) ) THEN
1.1 bertrand 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 )
1.6 bertrand 360: $ GO TO 10
1.1 bertrand 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
1.11 bertrand 376: VALUE = WORK( 0 )
377: DO I = 1, N-1
378: TEMP = WORK( I )
1.15 bertrand 379: IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
1.11 bertrand 380: $ VALUE = TEMP
381: END DO
1.1 bertrand 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
1.11 bertrand 418: VALUE = WORK( 0 )
419: DO I = 1, N-1
420: TEMP = WORK( I )
1.15 bertrand 421: IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
1.11 bertrand 422: $ VALUE = TEMP
423: END DO
1.1 bertrand 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
1.11 bertrand 456: VALUE = WORK( 0 )
457: DO I = 1, N-1
458: TEMP = WORK( I )
1.15 bertrand 459: IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
1.11 bertrand 460: $ VALUE = TEMP
461: END DO
1.1 bertrand 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
1.11 bertrand 494: VALUE = WORK( 0 )
495: DO I = 1, N-1
496: TEMP = WORK( I )
1.15 bertrand 497: IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
1.11 bertrand 498: $ VALUE = TEMP
499: END DO
1.1 bertrand 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
1.11 bertrand 560: VALUE = WORK( 0 )
561: DO I = 1, N-1
562: TEMP = WORK( I )
1.15 bertrand 563: IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
1.11 bertrand 564: $ VALUE = TEMP
565: END DO
1.1 bertrand 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
1.11 bertrand 625: VALUE = WORK( 0 )
626: DO I = 1, N-1
627: TEMP = WORK( I )
1.15 bertrand 628: IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
1.11 bertrand 629: $ VALUE = TEMP
630: END DO
1.1 bertrand 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
1.11 bertrand 698: VALUE = WORK( 0 )
699: DO I = 1, N-1
700: TEMP = WORK( I )
1.15 bertrand 701: IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
1.11 bertrand 702: $ VALUE = TEMP
703: END DO
1.1 bertrand 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
1.11 bertrand 771: VALUE = WORK( 0 )
772: DO I = 1, N-1
773: TEMP = WORK( I )
1.15 bertrand 774: IF( VALUE .LT. TEMP .OR. DISNAN( TEMP ) )
1.11 bertrand 775: $ VALUE = TEMP
776: END DO
1.1 bertrand 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
1.6 bertrand 827: * A**T is upper
1.1 bertrand 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,
1.6 bertrand 838: $ SCALE, S )
1.1 bertrand 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
1.6 bertrand 848: * A**T is lower
1.1 bertrand 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
1.6 bertrand 909: * A**T is upper
1.1 bertrand 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,
1.6 bertrand 920: $ S )
1.1 bertrand 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
1.6 bertrand 930: * A**T is lower
1.1 bertrand 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