File:
[local] /
rpl /
lapack /
lapack /
dlansf.f
Revision
1.13:
download - view:
text,
annotated -
select for diffs -
revision graph
Mon Jan 27 09:28:20 2014 UTC (10 years, 3 months ago) by
bertrand
Branches:
MAIN
CVS tags:
rpl-4_1_24,
rpl-4_1_23,
rpl-4_1_22,
rpl-4_1_21,
rpl-4_1_20,
rpl-4_1_19,
rpl-4_1_18,
rpl-4_1_17,
HEAD
Cohérence.
1: *> \brief \b DLANSF returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric matrix in RFP format.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DLANSF + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlansf.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlansf.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlansf.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * DOUBLE PRECISION FUNCTION DLANSF( NORM, TRANSR, UPLO, N, A, WORK )
22: *
23: * .. Scalar Arguments ..
24: * CHARACTER NORM, TRANSR, UPLO
25: * INTEGER N
26: * ..
27: * .. Array Arguments ..
28: * DOUBLE PRECISION A( 0: * ), WORK( 0: * )
29: * ..
30: *
31: *
32: *> \par Purpose:
33: * =============
34: *>
35: *> \verbatim
36: *>
37: *> DLANSF returns the value of the one norm, or the Frobenius norm, or
38: *> the infinity norm, or the element of largest absolute value of a
39: *> real symmetric matrix A in RFP format.
40: *> \endverbatim
41: *>
42: *> \return DLANSF
43: *> \verbatim
44: *>
45: *> DLANSF = ( max(abs(A(i,j))), NORM = 'M' or 'm'
46: *> (
47: *> ( norm1(A), NORM = '1', 'O' or 'o'
48: *> (
49: *> ( normI(A), NORM = 'I' or 'i'
50: *> (
51: *> ( normF(A), NORM = 'F', 'f', 'E' or 'e'
52: *>
53: *> where norm1 denotes the one norm of a matrix (maximum column sum),
54: *> normI denotes the infinity norm of a matrix (maximum row sum) and
55: *> normF denotes the Frobenius norm of a matrix (square root of sum of
56: *> squares). Note that max(abs(A(i,j))) is not a matrix norm.
57: *> \endverbatim
58: *
59: * Arguments:
60: * ==========
61: *
62: *> \param[in] NORM
63: *> \verbatim
64: *> NORM is CHARACTER*1
65: *> Specifies the value to be returned in DLANSF as described
66: *> above.
67: *> \endverbatim
68: *>
69: *> \param[in] TRANSR
70: *> \verbatim
71: *> TRANSR is CHARACTER*1
72: *> Specifies whether the RFP format of A is normal or
73: *> transposed format.
74: *> = 'N': RFP format is Normal;
75: *> = 'T': RFP format is Transpose.
76: *> \endverbatim
77: *>
78: *> \param[in] UPLO
79: *> \verbatim
80: *> UPLO is CHARACTER*1
81: *> On entry, UPLO specifies whether the RFP matrix A came from
82: *> an upper or lower triangular matrix as follows:
83: *> = 'U': RFP A came from an upper triangular matrix;
84: *> = 'L': RFP A came from a lower triangular matrix.
85: *> \endverbatim
86: *>
87: *> \param[in] N
88: *> \verbatim
89: *> N is INTEGER
90: *> The order of the matrix A. N >= 0. When N = 0, DLANSF is
91: *> set to zero.
92: *> \endverbatim
93: *>
94: *> \param[in] A
95: *> \verbatim
96: *> A is DOUBLE PRECISION array, dimension ( N*(N+1)/2 );
97: *> On entry, the upper (if UPLO = 'U') or lower (if UPLO = 'L')
98: *> part of the symmetric matrix A stored in RFP format. See the
99: *> "Notes" below for more details.
100: *> Unchanged on exit.
101: *> \endverbatim
102: *>
103: *> \param[out] WORK
104: *> \verbatim
105: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
106: *> where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
107: *> WORK is not referenced.
108: *> \endverbatim
109: *
110: * Authors:
111: * ========
112: *
113: *> \author Univ. of Tennessee
114: *> \author Univ. of California Berkeley
115: *> \author Univ. of Colorado Denver
116: *> \author NAG Ltd.
117: *
118: *> \date September 2012
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.2) --
213: * -- LAPACK is a software package provided by Univ. of Tennessee, --
214: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
215: * September 2012
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
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