![[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.7 bertrand 1: *> \brief \b DLANSF
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: * ===========
1.1 bertrand 20: *
1.7 bertrand 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 November 2011
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
1.1 bertrand 208: *
1.7 bertrand 209: * =====================================================================
210: DOUBLE PRECISION FUNCTION DLANSF( NORM, TRANSR, UPLO, N, A, WORK )
1.1 bertrand 211: *
1.7 bertrand 212: * -- LAPACK computational routine (version 3.4.0) --
1.1 bertrand 213: * -- LAPACK is a software package provided by Univ. of Tennessee, --
214: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.7 bertrand 215: * November 2011
1.1 bertrand 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
234: * ..
235: * .. External Functions ..
236: LOGICAL LSAME
237: INTEGER IDAMAX
238: EXTERNAL LSAME, IDAMAX
239: * ..
240: * .. External Subroutines ..
241: EXTERNAL DLASSQ
242: * ..
243: * .. Intrinsic Functions ..
244: INTRINSIC ABS, MAX, SQRT
245: * ..
246: * .. Executable Statements ..
247: *
248: IF( N.EQ.0 ) THEN
249: DLANSF = ZERO
250: RETURN
251: END IF
252: *
253: * set noe = 1 if n is odd. if n is even set noe=0
254: *
255: NOE = 1
256: IF( MOD( N, 2 ).EQ.0 )
1.6 bertrand 257: $ NOE = 0
1.1 bertrand 258: *
259: * set ifm = 0 when form='T or 't' and 1 otherwise
260: *
261: IFM = 1
262: IF( LSAME( TRANSR, 'T' ) )
1.6 bertrand 263: $ IFM = 0
1.1 bertrand 264: *
265: * set ilu = 0 when uplo='U or 'u' and 1 otherwise
266: *
267: ILU = 1
268: IF( LSAME( UPLO, 'U' ) )
1.6 bertrand 269: $ ILU = 0
1.1 bertrand 270: *
271: * set lda = (n+1)/2 when ifm = 0
272: * set lda = n when ifm = 1 and noe = 1
273: * set lda = n+1 when ifm = 1 and noe = 0
274: *
275: IF( IFM.EQ.1 ) THEN
276: IF( NOE.EQ.1 ) THEN
277: LDA = N
278: ELSE
279: * noe=0
280: LDA = N + 1
281: END IF
282: ELSE
283: * ifm=0
284: LDA = ( N+1 ) / 2
285: END IF
286: *
287: IF( LSAME( NORM, 'M' ) ) THEN
288: *
289: * Find max(abs(A(i,j))).
290: *
291: K = ( N+1 ) / 2
292: VALUE = ZERO
293: IF( NOE.EQ.1 ) THEN
294: * n is odd
295: IF( IFM.EQ.1 ) THEN
296: * A is n by k
297: DO J = 0, K - 1
298: DO I = 0, N - 1
299: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) )
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: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) )
307: END DO
308: END DO
309: END IF
310: ELSE
311: * n is even
312: IF( IFM.EQ.1 ) THEN
313: * A is n+1 by k
314: DO J = 0, K - 1
315: DO I = 0, N
316: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) )
317: END DO
318: END DO
319: ELSE
320: * xpose case; A is k by n+1
321: DO J = 0, N
322: DO I = 0, K - 1
323: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) )
324: END DO
325: END DO
326: END IF
327: END IF
328: ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR.
1.6 bertrand 329: $ ( NORM.EQ.'1' ) ) THEN
1.1 bertrand 330: *
331: * Find normI(A) ( = norm1(A), since A is symmetric).
332: *
333: IF( IFM.EQ.1 ) THEN
334: K = N / 2
335: IF( NOE.EQ.1 ) THEN
336: * n is odd
337: IF( ILU.EQ.0 ) THEN
338: DO I = 0, K - 1
339: WORK( I ) = ZERO
340: END DO
341: DO J = 0, K
342: S = ZERO
343: DO I = 0, K + J - 1
344: AA = ABS( A( I+J*LDA ) )
345: * -> A(i,j+k)
346: S = S + AA
347: WORK( I ) = WORK( I ) + AA
348: END DO
349: AA = ABS( A( I+J*LDA ) )
350: * -> A(j+k,j+k)
351: WORK( J+K ) = S + AA
352: IF( I.EQ.K+K )
1.6 bertrand 353: $ GO TO 10
1.1 bertrand 354: I = I + 1
355: AA = ABS( A( I+J*LDA ) )
356: * -> A(j,j)
357: WORK( J ) = WORK( J ) + AA
358: S = ZERO
359: DO L = J + 1, K - 1
360: I = I + 1
361: AA = ABS( A( I+J*LDA ) )
362: * -> A(l,j)
363: S = S + AA
364: WORK( L ) = WORK( L ) + AA
365: END DO
366: WORK( J ) = WORK( J ) + S
367: END DO
368: 10 CONTINUE
369: I = IDAMAX( N, WORK, 1 )
370: VALUE = WORK( I-1 )
371: ELSE
372: * ilu = 1
373: K = K + 1
374: * k=(n+1)/2 for n odd and ilu=1
375: DO I = K, N - 1
376: WORK( I ) = ZERO
377: END DO
378: DO J = K - 1, 0, -1
379: S = ZERO
380: DO I = 0, J - 2
381: AA = ABS( A( I+J*LDA ) )
382: * -> A(j+k,i+k)
383: S = S + AA
384: WORK( I+K ) = WORK( I+K ) + AA
385: END DO
386: IF( J.GT.0 ) THEN
387: AA = ABS( A( I+J*LDA ) )
388: * -> A(j+k,j+k)
389: S = S + AA
390: WORK( I+K ) = WORK( I+K ) + S
391: * i=j
392: I = I + 1
393: END IF
394: AA = ABS( A( I+J*LDA ) )
395: * -> A(j,j)
396: WORK( J ) = AA
397: S = ZERO
398: DO L = J + 1, N - 1
399: I = I + 1
400: AA = ABS( A( I+J*LDA ) )
401: * -> A(l,j)
402: S = S + AA
403: WORK( L ) = WORK( L ) + AA
404: END DO
405: WORK( J ) = WORK( J ) + S
406: END DO
407: I = IDAMAX( N, WORK, 1 )
408: VALUE = WORK( I-1 )
409: END IF
410: ELSE
411: * n is even
412: IF( ILU.EQ.0 ) THEN
413: DO I = 0, K - 1
414: WORK( I ) = ZERO
415: END DO
416: DO J = 0, K - 1
417: S = ZERO
418: DO I = 0, K + J - 1
419: AA = ABS( A( I+J*LDA ) )
420: * -> A(i,j+k)
421: S = S + AA
422: WORK( I ) = WORK( I ) + AA
423: END DO
424: AA = ABS( A( I+J*LDA ) )
425: * -> A(j+k,j+k)
426: WORK( J+K ) = S + AA
427: I = I + 1
428: AA = ABS( A( I+J*LDA ) )
429: * -> A(j,j)
430: WORK( J ) = WORK( J ) + AA
431: S = ZERO
432: DO L = J + 1, K - 1
433: I = I + 1
434: AA = ABS( A( I+J*LDA ) )
435: * -> A(l,j)
436: S = S + AA
437: WORK( L ) = WORK( L ) + AA
438: END DO
439: WORK( J ) = WORK( J ) + S
440: END DO
441: I = IDAMAX( N, WORK, 1 )
442: VALUE = WORK( I-1 )
443: ELSE
444: * ilu = 1
445: DO I = K, N - 1
446: WORK( I ) = ZERO
447: END DO
448: DO J = K - 1, 0, -1
449: S = ZERO
450: DO I = 0, J - 1
451: AA = ABS( A( I+J*LDA ) )
452: * -> A(j+k,i+k)
453: S = S + AA
454: WORK( I+K ) = WORK( I+K ) + AA
455: END DO
456: AA = ABS( A( I+J*LDA ) )
457: * -> A(j+k,j+k)
458: S = S + AA
459: WORK( I+K ) = WORK( I+K ) + S
460: * i=j
461: I = I + 1
462: AA = ABS( A( I+J*LDA ) )
463: * -> A(j,j)
464: WORK( J ) = AA
465: S = ZERO
466: DO L = J + 1, N - 1
467: I = I + 1
468: AA = ABS( A( I+J*LDA ) )
469: * -> A(l,j)
470: S = S + AA
471: WORK( L ) = WORK( L ) + AA
472: END DO
473: WORK( J ) = WORK( J ) + S
474: END DO
475: I = IDAMAX( N, WORK, 1 )
476: VALUE = WORK( I-1 )
477: END IF
478: END IF
479: ELSE
480: * ifm=0
481: K = N / 2
482: IF( NOE.EQ.1 ) THEN
483: * n is odd
484: IF( ILU.EQ.0 ) THEN
485: N1 = K
486: * n/2
487: K = K + 1
488: * k is the row size and lda
489: DO I = N1, N - 1
490: WORK( I ) = ZERO
491: END DO
492: DO J = 0, N1 - 1
493: S = ZERO
494: DO I = 0, K - 1
495: AA = ABS( A( I+J*LDA ) )
496: * A(j,n1+i)
497: WORK( I+N1 ) = WORK( I+N1 ) + AA
498: S = S + AA
499: END DO
500: WORK( J ) = S
501: END DO
502: * j=n1=k-1 is special
503: S = ABS( A( 0+J*LDA ) )
504: * A(k-1,k-1)
505: DO I = 1, K - 1
506: AA = ABS( A( I+J*LDA ) )
507: * A(k-1,i+n1)
508: WORK( I+N1 ) = WORK( I+N1 ) + AA
509: S = S + AA
510: END DO
511: WORK( J ) = WORK( J ) + S
512: DO J = K, N - 1
513: S = ZERO
514: DO I = 0, J - K - 1
515: AA = ABS( A( I+J*LDA ) )
516: * A(i,j-k)
517: WORK( I ) = WORK( I ) + AA
518: S = S + AA
519: END DO
520: * i=j-k
521: AA = ABS( A( I+J*LDA ) )
522: * A(j-k,j-k)
523: S = S + AA
524: WORK( J-K ) = WORK( J-K ) + S
525: I = I + 1
526: S = ABS( A( I+J*LDA ) )
527: * A(j,j)
528: DO L = J + 1, N - 1
529: I = I + 1
530: AA = ABS( A( I+J*LDA ) )
531: * A(j,l)
532: WORK( L ) = WORK( L ) + AA
533: S = S + AA
534: END DO
535: WORK( J ) = WORK( J ) + S
536: END DO
537: I = IDAMAX( N, WORK, 1 )
538: VALUE = WORK( I-1 )
539: ELSE
540: * ilu=1
541: K = K + 1
542: * k=(n+1)/2 for n odd and ilu=1
543: DO I = K, N - 1
544: WORK( I ) = ZERO
545: END DO
546: DO J = 0, K - 2
547: * process
548: S = ZERO
549: DO I = 0, J - 1
550: AA = ABS( A( I+J*LDA ) )
551: * A(j,i)
552: WORK( I ) = WORK( I ) + AA
553: S = S + AA
554: END DO
555: AA = ABS( A( I+J*LDA ) )
556: * i=j so process of A(j,j)
557: S = S + AA
558: WORK( J ) = S
559: * is initialised here
560: I = I + 1
561: * i=j process A(j+k,j+k)
562: AA = ABS( A( I+J*LDA ) )
563: S = AA
564: DO L = K + J + 1, N - 1
565: I = I + 1
566: AA = ABS( A( I+J*LDA ) )
567: * A(l,k+j)
568: S = S + AA
569: WORK( L ) = WORK( L ) + AA
570: END DO
571: WORK( K+J ) = WORK( K+J ) + S
572: END DO
573: * j=k-1 is special :process col A(k-1,0:k-1)
574: S = ZERO
575: DO I = 0, K - 2
576: AA = ABS( A( I+J*LDA ) )
577: * A(k,i)
578: WORK( I ) = WORK( I ) + AA
579: S = S + AA
580: END DO
581: * i=k-1
582: AA = ABS( A( I+J*LDA ) )
583: * A(k-1,k-1)
584: S = S + AA
585: WORK( I ) = S
586: * done with col j=k+1
587: DO J = K, N - 1
588: * process col j of A = A(j,0:k-1)
589: S = ZERO
590: DO I = 0, K - 1
591: AA = ABS( A( I+J*LDA ) )
592: * A(j,i)
593: WORK( I ) = WORK( I ) + AA
594: S = S + AA
595: END DO
596: WORK( J ) = WORK( J ) + S
597: END DO
598: I = IDAMAX( N, WORK, 1 )
599: VALUE = WORK( I-1 )
600: END IF
601: ELSE
602: * n is even
603: IF( ILU.EQ.0 ) THEN
604: DO I = K, N - 1
605: WORK( I ) = ZERO
606: END DO
607: DO J = 0, K - 1
608: S = ZERO
609: DO I = 0, K - 1
610: AA = ABS( A( I+J*LDA ) )
611: * A(j,i+k)
612: WORK( I+K ) = WORK( I+K ) + AA
613: S = S + AA
614: END DO
615: WORK( J ) = S
616: END DO
617: * j=k
618: AA = ABS( A( 0+J*LDA ) )
619: * A(k,k)
620: S = AA
621: DO I = 1, K - 1
622: AA = ABS( A( I+J*LDA ) )
623: * A(k,k+i)
624: WORK( I+K ) = WORK( I+K ) + AA
625: S = S + AA
626: END DO
627: WORK( J ) = WORK( J ) + S
628: DO J = K + 1, N - 1
629: S = ZERO
630: DO I = 0, J - 2 - K
631: AA = ABS( A( I+J*LDA ) )
632: * A(i,j-k-1)
633: WORK( I ) = WORK( I ) + AA
634: S = S + AA
635: END DO
636: * i=j-1-k
637: AA = ABS( A( I+J*LDA ) )
638: * A(j-k-1,j-k-1)
639: S = S + AA
640: WORK( J-K-1 ) = WORK( J-K-1 ) + S
641: I = I + 1
642: AA = ABS( A( I+J*LDA ) )
643: * A(j,j)
644: S = AA
645: DO L = J + 1, N - 1
646: I = I + 1
647: AA = ABS( A( I+J*LDA ) )
648: * A(j,l)
649: WORK( L ) = WORK( L ) + AA
650: S = S + AA
651: END DO
652: WORK( J ) = WORK( J ) + S
653: END DO
654: * j=n
655: S = ZERO
656: DO I = 0, K - 2
657: AA = ABS( A( I+J*LDA ) )
658: * A(i,k-1)
659: WORK( I ) = WORK( I ) + AA
660: S = S + AA
661: END DO
662: * i=k-1
663: AA = ABS( A( I+J*LDA ) )
664: * A(k-1,k-1)
665: S = S + AA
666: WORK( I ) = WORK( I ) + S
667: I = IDAMAX( N, WORK, 1 )
668: VALUE = WORK( I-1 )
669: ELSE
670: * ilu=1
671: DO I = K, N - 1
672: WORK( I ) = ZERO
673: END DO
674: * j=0 is special :process col A(k:n-1,k)
675: S = ABS( A( 0 ) )
676: * A(k,k)
677: DO I = 1, K - 1
678: AA = ABS( A( I ) )
679: * A(k+i,k)
680: WORK( I+K ) = WORK( I+K ) + AA
681: S = S + AA
682: END DO
683: WORK( K ) = WORK( K ) + S
684: DO J = 1, K - 1
685: * process
686: S = ZERO
687: DO I = 0, J - 2
688: AA = ABS( A( I+J*LDA ) )
689: * A(j-1,i)
690: WORK( I ) = WORK( I ) + AA
691: S = S + AA
692: END DO
693: AA = ABS( A( I+J*LDA ) )
694: * i=j-1 so process of A(j-1,j-1)
695: S = S + AA
696: WORK( J-1 ) = S
697: * is initialised here
698: I = I + 1
699: * i=j process A(j+k,j+k)
700: AA = ABS( A( I+J*LDA ) )
701: S = AA
702: DO L = K + J + 1, N - 1
703: I = I + 1
704: AA = ABS( A( I+J*LDA ) )
705: * A(l,k+j)
706: S = S + AA
707: WORK( L ) = WORK( L ) + AA
708: END DO
709: WORK( K+J ) = WORK( K+J ) + S
710: END DO
711: * j=k is special :process col A(k,0:k-1)
712: S = ZERO
713: DO I = 0, K - 2
714: AA = ABS( A( I+J*LDA ) )
715: * A(k,i)
716: WORK( I ) = WORK( I ) + AA
717: S = S + AA
718: END DO
719: * i=k-1
720: AA = ABS( A( I+J*LDA ) )
721: * A(k-1,k-1)
722: S = S + AA
723: WORK( I ) = S
724: * done with col j=k+1
725: DO J = K + 1, N
726: * process col j-1 of A = A(j-1,0:k-1)
727: S = ZERO
728: DO I = 0, K - 1
729: AA = ABS( A( I+J*LDA ) )
730: * A(j-1,i)
731: WORK( I ) = WORK( I ) + AA
732: S = S + AA
733: END DO
734: WORK( J-1 ) = WORK( J-1 ) + S
735: END DO
736: I = IDAMAX( N, WORK, 1 )
737: VALUE = WORK( I-1 )
738: END IF
739: END IF
740: END IF
741: ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
742: *
743: * Find normF(A).
744: *
745: K = ( N+1 ) / 2
746: SCALE = ZERO
747: S = ONE
748: IF( NOE.EQ.1 ) THEN
749: * n is odd
750: IF( IFM.EQ.1 ) THEN
751: * A is normal
752: IF( ILU.EQ.0 ) THEN
753: * A is upper
754: DO J = 0, K - 3
755: CALL DLASSQ( K-J-2, A( K+J+1+J*LDA ), 1, SCALE, S )
756: * L at A(k,0)
757: END DO
758: DO J = 0, K - 1
759: CALL DLASSQ( K+J-1, A( 0+J*LDA ), 1, SCALE, S )
760: * trap U at A(0,0)
761: END DO
762: S = S + S
763: * double s for the off diagonal elements
764: CALL DLASSQ( K-1, A( K ), LDA+1, SCALE, S )
765: * tri L at A(k,0)
766: CALL DLASSQ( K, A( K-1 ), LDA+1, SCALE, S )
767: * tri U at A(k-1,0)
768: ELSE
769: * ilu=1 & A is lower
770: DO J = 0, K - 1
771: CALL DLASSQ( N-J-1, A( J+1+J*LDA ), 1, SCALE, S )
772: * trap L at A(0,0)
773: END DO
774: DO J = 0, K - 2
775: CALL DLASSQ( J, A( 0+( 1+J )*LDA ), 1, SCALE, S )
776: * U at A(0,1)
777: END DO
778: S = S + S
779: * double s for the off diagonal elements
780: CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S )
781: * tri L at A(0,0)
782: CALL DLASSQ( K-1, A( 0+LDA ), LDA+1, SCALE, S )
783: * tri U at A(0,1)
784: END IF
785: ELSE
786: * A is xpose
787: IF( ILU.EQ.0 ) THEN
1.6 bertrand 788: * A**T is upper
1.1 bertrand 789: DO J = 1, K - 2
790: CALL DLASSQ( J, A( 0+( K+J )*LDA ), 1, SCALE, S )
791: * U at A(0,k)
792: END DO
793: DO J = 0, K - 2
794: CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
795: * k by k-1 rect. at A(0,0)
796: END DO
797: DO J = 0, K - 2
798: CALL DLASSQ( K-J-1, A( J+1+( J+K-1 )*LDA ), 1,
1.6 bertrand 799: $ SCALE, S )
1.1 bertrand 800: * L at A(0,k-1)
801: END DO
802: S = S + S
803: * double s for the off diagonal elements
804: CALL DLASSQ( K-1, A( 0+K*LDA ), LDA+1, SCALE, S )
805: * tri U at A(0,k)
806: CALL DLASSQ( K, A( 0+( K-1 )*LDA ), LDA+1, SCALE, S )
807: * tri L at A(0,k-1)
808: ELSE
1.6 bertrand 809: * A**T is lower
1.1 bertrand 810: DO J = 1, K - 1
811: CALL DLASSQ( J, A( 0+J*LDA ), 1, SCALE, S )
812: * U at A(0,0)
813: END DO
814: DO J = K, N - 1
815: CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
816: * k by k-1 rect. at A(0,k)
817: END DO
818: DO J = 0, K - 3
819: CALL DLASSQ( K-J-2, A( J+2+J*LDA ), 1, SCALE, S )
820: * L at A(1,0)
821: END DO
822: S = S + S
823: * double s for the off diagonal elements
824: CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S )
825: * tri U at A(0,0)
826: CALL DLASSQ( K-1, A( 1 ), LDA+1, SCALE, S )
827: * tri L at A(1,0)
828: END IF
829: END IF
830: ELSE
831: * n is even
832: IF( IFM.EQ.1 ) THEN
833: * A is normal
834: IF( ILU.EQ.0 ) THEN
835: * A is upper
836: DO J = 0, K - 2
837: CALL DLASSQ( K-J-1, A( K+J+2+J*LDA ), 1, SCALE, S )
838: * L at A(k+1,0)
839: END DO
840: DO J = 0, K - 1
841: CALL DLASSQ( K+J, A( 0+J*LDA ), 1, SCALE, S )
842: * trap U at A(0,0)
843: END DO
844: S = S + S
845: * double s for the off diagonal elements
846: CALL DLASSQ( K, A( K+1 ), LDA+1, SCALE, S )
847: * tri L at A(k+1,0)
848: CALL DLASSQ( K, A( K ), LDA+1, SCALE, S )
849: * tri U at A(k,0)
850: ELSE
851: * ilu=1 & A is lower
852: DO J = 0, K - 1
853: CALL DLASSQ( N-J-1, A( J+2+J*LDA ), 1, SCALE, S )
854: * trap L at A(1,0)
855: END DO
856: DO J = 1, K - 1
857: CALL DLASSQ( J, A( 0+J*LDA ), 1, SCALE, S )
858: * U at A(0,0)
859: END DO
860: S = S + S
861: * double s for the off diagonal elements
862: CALL DLASSQ( K, A( 1 ), LDA+1, SCALE, S )
863: * tri L at A(1,0)
864: CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S )
865: * tri U at A(0,0)
866: END IF
867: ELSE
868: * A is xpose
869: IF( ILU.EQ.0 ) THEN
1.6 bertrand 870: * A**T is upper
1.1 bertrand 871: DO J = 1, K - 1
872: CALL DLASSQ( J, A( 0+( K+1+J )*LDA ), 1, SCALE, S )
873: * U at A(0,k+1)
874: END DO
875: DO J = 0, K - 1
876: CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
877: * k by k rect. at A(0,0)
878: END DO
879: DO J = 0, K - 2
880: CALL DLASSQ( K-J-1, A( J+1+( J+K )*LDA ), 1, SCALE,
1.6 bertrand 881: $ S )
1.1 bertrand 882: * L at A(0,k)
883: END DO
884: S = S + S
885: * double s for the off diagonal elements
886: CALL DLASSQ( K, A( 0+( K+1 )*LDA ), LDA+1, SCALE, S )
887: * tri U at A(0,k+1)
888: CALL DLASSQ( K, A( 0+K*LDA ), LDA+1, SCALE, S )
889: * tri L at A(0,k)
890: ELSE
1.6 bertrand 891: * A**T is lower
1.1 bertrand 892: DO J = 1, K - 1
893: CALL DLASSQ( J, A( 0+( J+1 )*LDA ), 1, SCALE, S )
894: * U at A(0,1)
895: END DO
896: DO J = K + 1, N
897: CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
898: * k by k rect. at A(0,k+1)
899: END DO
900: DO J = 0, K - 2
901: CALL DLASSQ( K-J-1, A( J+1+J*LDA ), 1, SCALE, S )
902: * L at A(0,0)
903: END DO
904: S = S + S
905: * double s for the off diagonal elements
906: CALL DLASSQ( K, A( LDA ), LDA+1, SCALE, S )
907: * tri L at A(0,1)
908: CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S )
909: * tri U at A(0,0)
910: END IF
911: END IF
912: END IF
913: VALUE = SCALE*SQRT( S )
914: END IF
915: *
916: DLANSF = VALUE
917: RETURN
918: *
919: * End of DLANSF
920: *
921: END