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