![[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.1 bertrand 1: DOUBLE PRECISION FUNCTION DLANSF( NORM, TRANSR, UPLO, N, A, WORK )
2: *
1.6 ! bertrand 3: * -- LAPACK routine (version 3.3.1) --
1.1 bertrand 4: *
5: * -- Contributed by Fred Gustavson of the IBM Watson Research Center --
1.6 ! bertrand 6: * -- April 2011 --
1.1 bertrand 7: *
8: * -- LAPACK is a software package provided by Univ. of Tennessee, --
9: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
10: *
11: * .. Scalar Arguments ..
12: CHARACTER NORM, TRANSR, UPLO
13: INTEGER N
14: * ..
15: * .. Array Arguments ..
16: DOUBLE PRECISION A( 0: * ), WORK( 0: * )
17: * ..
18: *
19: * Purpose
20: * =======
21: *
22: * DLANSF returns the value of the one norm, or the Frobenius norm, or
23: * the infinity norm, or the element of largest absolute value of a
24: * real symmetric matrix A in RFP format.
25: *
26: * Description
27: * ===========
28: *
29: * DLANSF returns the value
30: *
31: * DLANSF = ( max(abs(A(i,j))), NORM = 'M' or 'm'
32: * (
33: * ( norm1(A), NORM = '1', 'O' or 'o'
34: * (
35: * ( normI(A), NORM = 'I' or 'i'
36: * (
37: * ( normF(A), NORM = 'F', 'f', 'E' or 'e'
38: *
39: * where norm1 denotes the one norm of a matrix (maximum column sum),
40: * normI denotes the infinity norm of a matrix (maximum row sum) and
41: * normF denotes the Frobenius norm of a matrix (square root of sum of
42: * squares). Note that max(abs(A(i,j))) is not a matrix norm.
43: *
44: * Arguments
45: * =========
46: *
1.4 bertrand 47: * NORM (input) CHARACTER*1
1.1 bertrand 48: * Specifies the value to be returned in DLANSF as described
49: * above.
50: *
1.4 bertrand 51: * TRANSR (input) CHARACTER*1
1.1 bertrand 52: * Specifies whether the RFP format of A is normal or
53: * transposed format.
54: * = 'N': RFP format is Normal;
55: * = 'T': RFP format is Transpose.
56: *
1.4 bertrand 57: * UPLO (input) CHARACTER*1
1.1 bertrand 58: * On entry, UPLO specifies whether the RFP matrix A came from
59: * an upper or lower triangular matrix as follows:
60: * = 'U': RFP A came from an upper triangular matrix;
61: * = 'L': RFP A came from a lower triangular matrix.
62: *
63: * N (input) INTEGER
64: * The order of the matrix A. N >= 0. When N = 0, DLANSF is
65: * set to zero.
66: *
67: * A (input) DOUBLE PRECISION array, dimension ( N*(N+1)/2 );
68: * On entry, the upper (if UPLO = 'U') or lower (if UPLO = 'L')
69: * part of the symmetric matrix A stored in RFP format. See the
70: * "Notes" below for more details.
71: * Unchanged on exit.
72: *
73: * WORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
74: * where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
75: * WORK is not referenced.
76: *
77: * Further Details
78: * ===============
79: *
80: * We first consider Rectangular Full Packed (RFP) Format when N is
81: * even. We give an example where N = 6.
82: *
83: * AP is Upper AP is Lower
84: *
85: * 00 01 02 03 04 05 00
86: * 11 12 13 14 15 10 11
87: * 22 23 24 25 20 21 22
88: * 33 34 35 30 31 32 33
89: * 44 45 40 41 42 43 44
90: * 55 50 51 52 53 54 55
91: *
92: *
93: * Let TRANSR = 'N'. RFP holds AP as follows:
94: * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
95: * three columns of AP upper. The lower triangle A(4:6,0:2) consists of
96: * the transpose of the first three columns of AP upper.
97: * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
98: * three columns of AP lower. The upper triangle A(0:2,0:2) consists of
99: * the transpose of the last three columns of AP lower.
100: * This covers the case N even and TRANSR = 'N'.
101: *
102: * RFP A RFP A
103: *
104: * 03 04 05 33 43 53
105: * 13 14 15 00 44 54
106: * 23 24 25 10 11 55
107: * 33 34 35 20 21 22
108: * 00 44 45 30 31 32
109: * 01 11 55 40 41 42
110: * 02 12 22 50 51 52
111: *
112: * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
113: * transpose of RFP A above. One therefore gets:
114: *
115: *
116: * RFP A RFP A
117: *
118: * 03 13 23 33 00 01 02 33 00 10 20 30 40 50
119: * 04 14 24 34 44 11 12 43 44 11 21 31 41 51
120: * 05 15 25 35 45 55 22 53 54 55 22 32 42 52
121: *
122: *
123: * We then consider Rectangular Full Packed (RFP) Format when N is
124: * odd. We give an example where N = 5.
125: *
126: * AP is Upper AP is Lower
127: *
128: * 00 01 02 03 04 00
129: * 11 12 13 14 10 11
130: * 22 23 24 20 21 22
131: * 33 34 30 31 32 33
132: * 44 40 41 42 43 44
133: *
134: *
135: * Let TRANSR = 'N'. RFP holds AP as follows:
136: * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
137: * three columns of AP upper. The lower triangle A(3:4,0:1) consists of
138: * the transpose of the first two columns of AP upper.
139: * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
140: * three columns of AP lower. The upper triangle A(0:1,1:2) consists of
141: * the transpose of the last two columns of AP lower.
142: * This covers the case N odd and TRANSR = 'N'.
143: *
144: * RFP A RFP A
145: *
146: * 02 03 04 00 33 43
147: * 12 13 14 10 11 44
148: * 22 23 24 20 21 22
149: * 00 33 34 30 31 32
150: * 01 11 44 40 41 42
151: *
152: * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
153: * transpose of RFP A above. One therefore gets:
154: *
155: * RFP A RFP A
156: *
157: * 02 12 22 00 01 00 10 20 30 40 50
158: * 03 13 23 33 11 33 11 21 31 41 51
159: * 04 14 24 34 44 43 44 22 32 42 52
160: *
161: * Reference
162: * =========
163: *
164: * =====================================================================
165: *
166: * .. Parameters ..
167: DOUBLE PRECISION ONE, ZERO
168: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
169: * ..
170: * .. Local Scalars ..
171: INTEGER I, J, IFM, ILU, NOE, N1, K, L, LDA
172: DOUBLE PRECISION SCALE, S, VALUE, AA
173: * ..
174: * .. External Functions ..
175: LOGICAL LSAME
176: INTEGER IDAMAX
177: EXTERNAL LSAME, IDAMAX
178: * ..
179: * .. External Subroutines ..
180: EXTERNAL DLASSQ
181: * ..
182: * .. Intrinsic Functions ..
183: INTRINSIC ABS, MAX, SQRT
184: * ..
185: * .. Executable Statements ..
186: *
187: IF( N.EQ.0 ) THEN
188: DLANSF = ZERO
189: RETURN
190: END IF
191: *
192: * set noe = 1 if n is odd. if n is even set noe=0
193: *
194: NOE = 1
195: IF( MOD( N, 2 ).EQ.0 )
1.6 ! bertrand 196: $ NOE = 0
1.1 bertrand 197: *
198: * set ifm = 0 when form='T or 't' and 1 otherwise
199: *
200: IFM = 1
201: IF( LSAME( TRANSR, 'T' ) )
1.6 ! bertrand 202: $ IFM = 0
1.1 bertrand 203: *
204: * set ilu = 0 when uplo='U or 'u' and 1 otherwise
205: *
206: ILU = 1
207: IF( LSAME( UPLO, 'U' ) )
1.6 ! bertrand 208: $ ILU = 0
1.1 bertrand 209: *
210: * set lda = (n+1)/2 when ifm = 0
211: * set lda = n when ifm = 1 and noe = 1
212: * set lda = n+1 when ifm = 1 and noe = 0
213: *
214: IF( IFM.EQ.1 ) THEN
215: IF( NOE.EQ.1 ) THEN
216: LDA = N
217: ELSE
218: * noe=0
219: LDA = N + 1
220: END IF
221: ELSE
222: * ifm=0
223: LDA = ( N+1 ) / 2
224: END IF
225: *
226: IF( LSAME( NORM, 'M' ) ) THEN
227: *
228: * Find max(abs(A(i,j))).
229: *
230: K = ( N+1 ) / 2
231: VALUE = ZERO
232: IF( NOE.EQ.1 ) THEN
233: * n is odd
234: IF( IFM.EQ.1 ) THEN
235: * A is n by k
236: DO J = 0, K - 1
237: DO I = 0, N - 1
238: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) )
239: END DO
240: END DO
241: ELSE
242: * xpose case; A is k by n
243: DO J = 0, N - 1
244: DO I = 0, K - 1
245: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) )
246: END DO
247: END DO
248: END IF
249: ELSE
250: * n is even
251: IF( IFM.EQ.1 ) THEN
252: * A is n+1 by k
253: DO J = 0, K - 1
254: DO I = 0, N
255: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) )
256: END DO
257: END DO
258: ELSE
259: * xpose case; A is k by n+1
260: DO J = 0, N
261: DO I = 0, K - 1
262: VALUE = MAX( VALUE, ABS( A( I+J*LDA ) ) )
263: END DO
264: END DO
265: END IF
266: END IF
267: ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR.
1.6 ! bertrand 268: $ ( NORM.EQ.'1' ) ) THEN
1.1 bertrand 269: *
270: * Find normI(A) ( = norm1(A), since A is symmetric).
271: *
272: IF( IFM.EQ.1 ) THEN
273: K = N / 2
274: IF( NOE.EQ.1 ) THEN
275: * n is odd
276: IF( ILU.EQ.0 ) THEN
277: DO I = 0, K - 1
278: WORK( I ) = ZERO
279: END DO
280: DO J = 0, K
281: S = ZERO
282: DO I = 0, K + J - 1
283: AA = ABS( A( I+J*LDA ) )
284: * -> A(i,j+k)
285: S = S + AA
286: WORK( I ) = WORK( I ) + AA
287: END DO
288: AA = ABS( A( I+J*LDA ) )
289: * -> A(j+k,j+k)
290: WORK( J+K ) = S + AA
291: IF( I.EQ.K+K )
1.6 ! bertrand 292: $ GO TO 10
1.1 bertrand 293: I = I + 1
294: AA = ABS( A( I+J*LDA ) )
295: * -> A(j,j)
296: WORK( J ) = WORK( J ) + AA
297: S = ZERO
298: DO L = J + 1, K - 1
299: I = I + 1
300: AA = ABS( A( I+J*LDA ) )
301: * -> A(l,j)
302: S = S + AA
303: WORK( L ) = WORK( L ) + AA
304: END DO
305: WORK( J ) = WORK( J ) + S
306: END DO
307: 10 CONTINUE
308: I = IDAMAX( N, WORK, 1 )
309: VALUE = WORK( I-1 )
310: ELSE
311: * ilu = 1
312: K = K + 1
313: * k=(n+1)/2 for n odd and ilu=1
314: DO I = K, N - 1
315: WORK( I ) = ZERO
316: END DO
317: DO J = K - 1, 0, -1
318: S = ZERO
319: DO I = 0, J - 2
320: AA = ABS( A( I+J*LDA ) )
321: * -> A(j+k,i+k)
322: S = S + AA
323: WORK( I+K ) = WORK( I+K ) + AA
324: END DO
325: IF( J.GT.0 ) THEN
326: AA = ABS( A( I+J*LDA ) )
327: * -> A(j+k,j+k)
328: S = S + AA
329: WORK( I+K ) = WORK( I+K ) + S
330: * i=j
331: I = I + 1
332: END IF
333: AA = ABS( A( I+J*LDA ) )
334: * -> A(j,j)
335: WORK( J ) = AA
336: S = ZERO
337: DO L = J + 1, N - 1
338: I = I + 1
339: AA = ABS( A( I+J*LDA ) )
340: * -> A(l,j)
341: S = S + AA
342: WORK( L ) = WORK( L ) + AA
343: END DO
344: WORK( J ) = WORK( J ) + S
345: END DO
346: I = IDAMAX( N, WORK, 1 )
347: VALUE = WORK( I-1 )
348: END IF
349: ELSE
350: * n is even
351: IF( ILU.EQ.0 ) THEN
352: DO I = 0, K - 1
353: WORK( I ) = ZERO
354: END DO
355: DO J = 0, K - 1
356: S = ZERO
357: DO I = 0, K + J - 1
358: AA = ABS( A( I+J*LDA ) )
359: * -> A(i,j+k)
360: S = S + AA
361: WORK( I ) = WORK( I ) + AA
362: END DO
363: AA = ABS( A( I+J*LDA ) )
364: * -> A(j+k,j+k)
365: WORK( J+K ) = S + AA
366: I = I + 1
367: AA = ABS( A( I+J*LDA ) )
368: * -> A(j,j)
369: WORK( J ) = WORK( J ) + AA
370: S = ZERO
371: DO L = J + 1, K - 1
372: I = I + 1
373: AA = ABS( A( I+J*LDA ) )
374: * -> A(l,j)
375: S = S + AA
376: WORK( L ) = WORK( L ) + AA
377: END DO
378: WORK( J ) = WORK( J ) + S
379: END DO
380: I = IDAMAX( N, WORK, 1 )
381: VALUE = WORK( I-1 )
382: ELSE
383: * ilu = 1
384: DO I = K, N - 1
385: WORK( I ) = ZERO
386: END DO
387: DO J = K - 1, 0, -1
388: S = ZERO
389: DO I = 0, J - 1
390: AA = ABS( A( I+J*LDA ) )
391: * -> A(j+k,i+k)
392: S = S + AA
393: WORK( I+K ) = WORK( I+K ) + AA
394: END DO
395: AA = ABS( A( I+J*LDA ) )
396: * -> A(j+k,j+k)
397: S = S + AA
398: WORK( I+K ) = WORK( I+K ) + S
399: * i=j
400: I = I + 1
401: AA = ABS( A( I+J*LDA ) )
402: * -> A(j,j)
403: WORK( J ) = AA
404: S = ZERO
405: DO L = J + 1, N - 1
406: I = I + 1
407: AA = ABS( A( I+J*LDA ) )
408: * -> A(l,j)
409: S = S + AA
410: WORK( L ) = WORK( L ) + AA
411: END DO
412: WORK( J ) = WORK( J ) + S
413: END DO
414: I = IDAMAX( N, WORK, 1 )
415: VALUE = WORK( I-1 )
416: END IF
417: END IF
418: ELSE
419: * ifm=0
420: K = N / 2
421: IF( NOE.EQ.1 ) THEN
422: * n is odd
423: IF( ILU.EQ.0 ) THEN
424: N1 = K
425: * n/2
426: K = K + 1
427: * k is the row size and lda
428: DO I = N1, N - 1
429: WORK( I ) = ZERO
430: END DO
431: DO J = 0, N1 - 1
432: S = ZERO
433: DO I = 0, K - 1
434: AA = ABS( A( I+J*LDA ) )
435: * A(j,n1+i)
436: WORK( I+N1 ) = WORK( I+N1 ) + AA
437: S = S + AA
438: END DO
439: WORK( J ) = S
440: END DO
441: * j=n1=k-1 is special
442: S = ABS( A( 0+J*LDA ) )
443: * A(k-1,k-1)
444: DO I = 1, K - 1
445: AA = ABS( A( I+J*LDA ) )
446: * A(k-1,i+n1)
447: WORK( I+N1 ) = WORK( I+N1 ) + AA
448: S = S + AA
449: END DO
450: WORK( J ) = WORK( J ) + S
451: DO J = K, N - 1
452: S = ZERO
453: DO I = 0, J - K - 1
454: AA = ABS( A( I+J*LDA ) )
455: * A(i,j-k)
456: WORK( I ) = WORK( I ) + AA
457: S = S + AA
458: END DO
459: * i=j-k
460: AA = ABS( A( I+J*LDA ) )
461: * A(j-k,j-k)
462: S = S + AA
463: WORK( J-K ) = WORK( J-K ) + S
464: I = I + 1
465: S = ABS( A( I+J*LDA ) )
466: * A(j,j)
467: DO L = J + 1, N - 1
468: I = I + 1
469: AA = ABS( A( I+J*LDA ) )
470: * A(j,l)
471: WORK( L ) = WORK( L ) + AA
472: S = S + AA
473: END DO
474: WORK( J ) = WORK( J ) + S
475: END DO
476: I = IDAMAX( N, WORK, 1 )
477: VALUE = WORK( I-1 )
478: ELSE
479: * ilu=1
480: K = K + 1
481: * k=(n+1)/2 for n odd and ilu=1
482: DO I = K, N - 1
483: WORK( I ) = ZERO
484: END DO
485: DO J = 0, K - 2
486: * process
487: S = ZERO
488: DO I = 0, J - 1
489: AA = ABS( A( I+J*LDA ) )
490: * A(j,i)
491: WORK( I ) = WORK( I ) + AA
492: S = S + AA
493: END DO
494: AA = ABS( A( I+J*LDA ) )
495: * i=j so process of A(j,j)
496: S = S + AA
497: WORK( J ) = S
498: * is initialised here
499: I = I + 1
500: * i=j process A(j+k,j+k)
501: AA = ABS( A( I+J*LDA ) )
502: S = AA
503: DO L = K + J + 1, N - 1
504: I = I + 1
505: AA = ABS( A( I+J*LDA ) )
506: * A(l,k+j)
507: S = S + AA
508: WORK( L ) = WORK( L ) + AA
509: END DO
510: WORK( K+J ) = WORK( K+J ) + S
511: END DO
512: * j=k-1 is special :process col A(k-1,0:k-1)
513: S = ZERO
514: DO I = 0, K - 2
515: AA = ABS( A( I+J*LDA ) )
516: * A(k,i)
517: WORK( I ) = WORK( I ) + AA
518: S = S + AA
519: END DO
520: * i=k-1
521: AA = ABS( A( I+J*LDA ) )
522: * A(k-1,k-1)
523: S = S + AA
524: WORK( I ) = S
525: * done with col j=k+1
526: DO J = K, N - 1
527: * process col j of A = A(j,0:k-1)
528: S = ZERO
529: DO I = 0, K - 1
530: AA = ABS( A( I+J*LDA ) )
531: * A(j,i)
532: WORK( I ) = WORK( I ) + 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: END IF
540: ELSE
541: * n is even
542: IF( ILU.EQ.0 ) THEN
543: DO I = K, N - 1
544: WORK( I ) = ZERO
545: END DO
546: DO J = 0, K - 1
547: S = ZERO
548: DO I = 0, K - 1
549: AA = ABS( A( I+J*LDA ) )
550: * A(j,i+k)
551: WORK( I+K ) = WORK( I+K ) + AA
552: S = S + AA
553: END DO
554: WORK( J ) = S
555: END DO
556: * j=k
557: AA = ABS( A( 0+J*LDA ) )
558: * A(k,k)
559: S = AA
560: DO I = 1, K - 1
561: AA = ABS( A( I+J*LDA ) )
562: * A(k,k+i)
563: WORK( I+K ) = WORK( I+K ) + AA
564: S = S + AA
565: END DO
566: WORK( J ) = WORK( J ) + S
567: DO J = K + 1, N - 1
568: S = ZERO
569: DO I = 0, J - 2 - K
570: AA = ABS( A( I+J*LDA ) )
571: * A(i,j-k-1)
572: WORK( I ) = WORK( I ) + AA
573: S = S + AA
574: END DO
575: * i=j-1-k
576: AA = ABS( A( I+J*LDA ) )
577: * A(j-k-1,j-k-1)
578: S = S + AA
579: WORK( J-K-1 ) = WORK( J-K-1 ) + S
580: I = I + 1
581: AA = ABS( A( I+J*LDA ) )
582: * A(j,j)
583: S = AA
584: DO L = J + 1, N - 1
585: I = I + 1
586: AA = ABS( A( I+J*LDA ) )
587: * A(j,l)
588: WORK( L ) = WORK( L ) + AA
589: S = S + AA
590: END DO
591: WORK( J ) = WORK( J ) + S
592: END DO
593: * j=n
594: S = ZERO
595: DO I = 0, K - 2
596: AA = ABS( A( I+J*LDA ) )
597: * A(i,k-1)
598: WORK( I ) = WORK( I ) + AA
599: S = S + AA
600: END DO
601: * i=k-1
602: AA = ABS( A( I+J*LDA ) )
603: * A(k-1,k-1)
604: S = S + AA
605: WORK( I ) = WORK( I ) + S
606: I = IDAMAX( N, WORK, 1 )
607: VALUE = WORK( I-1 )
608: ELSE
609: * ilu=1
610: DO I = K, N - 1
611: WORK( I ) = ZERO
612: END DO
613: * j=0 is special :process col A(k:n-1,k)
614: S = ABS( A( 0 ) )
615: * A(k,k)
616: DO I = 1, K - 1
617: AA = ABS( A( I ) )
618: * A(k+i,k)
619: WORK( I+K ) = WORK( I+K ) + AA
620: S = S + AA
621: END DO
622: WORK( K ) = WORK( K ) + S
623: DO J = 1, K - 1
624: * process
625: S = ZERO
626: DO I = 0, J - 2
627: AA = ABS( A( I+J*LDA ) )
628: * A(j-1,i)
629: WORK( I ) = WORK( I ) + AA
630: S = S + AA
631: END DO
632: AA = ABS( A( I+J*LDA ) )
633: * i=j-1 so process of A(j-1,j-1)
634: S = S + AA
635: WORK( J-1 ) = S
636: * is initialised here
637: I = I + 1
638: * i=j process A(j+k,j+k)
639: AA = ABS( A( I+J*LDA ) )
640: S = AA
641: DO L = K + J + 1, N - 1
642: I = I + 1
643: AA = ABS( A( I+J*LDA ) )
644: * A(l,k+j)
645: S = S + AA
646: WORK( L ) = WORK( L ) + AA
647: END DO
648: WORK( K+J ) = WORK( K+J ) + S
649: END DO
650: * j=k is special :process col A(k,0:k-1)
651: S = ZERO
652: DO I = 0, K - 2
653: AA = ABS( A( I+J*LDA ) )
654: * A(k,i)
655: WORK( I ) = WORK( I ) + AA
656: S = S + AA
657: END DO
658: * i=k-1
659: AA = ABS( A( I+J*LDA ) )
660: * A(k-1,k-1)
661: S = S + AA
662: WORK( I ) = S
663: * done with col j=k+1
664: DO J = K + 1, N
665: * process col j-1 of A = A(j-1,0:k-1)
666: S = ZERO
667: DO I = 0, K - 1
668: AA = ABS( A( I+J*LDA ) )
669: * A(j-1,i)
670: WORK( I ) = WORK( I ) + AA
671: S = S + AA
672: END DO
673: WORK( J-1 ) = WORK( J-1 ) + S
674: END DO
675: I = IDAMAX( N, WORK, 1 )
676: VALUE = WORK( I-1 )
677: END IF
678: END IF
679: END IF
680: ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
681: *
682: * Find normF(A).
683: *
684: K = ( N+1 ) / 2
685: SCALE = ZERO
686: S = ONE
687: IF( NOE.EQ.1 ) THEN
688: * n is odd
689: IF( IFM.EQ.1 ) THEN
690: * A is normal
691: IF( ILU.EQ.0 ) THEN
692: * A is upper
693: DO J = 0, K - 3
694: CALL DLASSQ( K-J-2, A( K+J+1+J*LDA ), 1, SCALE, S )
695: * L at A(k,0)
696: END DO
697: DO J = 0, K - 1
698: CALL DLASSQ( K+J-1, A( 0+J*LDA ), 1, SCALE, S )
699: * trap U at A(0,0)
700: END DO
701: S = S + S
702: * double s for the off diagonal elements
703: CALL DLASSQ( K-1, A( K ), LDA+1, SCALE, S )
704: * tri L at A(k,0)
705: CALL DLASSQ( K, A( K-1 ), LDA+1, SCALE, S )
706: * tri U at A(k-1,0)
707: ELSE
708: * ilu=1 & A is lower
709: DO J = 0, K - 1
710: CALL DLASSQ( N-J-1, A( J+1+J*LDA ), 1, SCALE, S )
711: * trap L at A(0,0)
712: END DO
713: DO J = 0, K - 2
714: CALL DLASSQ( J, A( 0+( 1+J )*LDA ), 1, SCALE, S )
715: * U at A(0,1)
716: END DO
717: S = S + S
718: * double s for the off diagonal elements
719: CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S )
720: * tri L at A(0,0)
721: CALL DLASSQ( K-1, A( 0+LDA ), LDA+1, SCALE, S )
722: * tri U at A(0,1)
723: END IF
724: ELSE
725: * A is xpose
726: IF( ILU.EQ.0 ) THEN
1.6 ! bertrand 727: * A**T is upper
1.1 bertrand 728: DO J = 1, K - 2
729: CALL DLASSQ( J, A( 0+( K+J )*LDA ), 1, SCALE, S )
730: * U at A(0,k)
731: END DO
732: DO J = 0, K - 2
733: CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
734: * k by k-1 rect. at A(0,0)
735: END DO
736: DO J = 0, K - 2
737: CALL DLASSQ( K-J-1, A( J+1+( J+K-1 )*LDA ), 1,
1.6 ! bertrand 738: $ SCALE, S )
1.1 bertrand 739: * L at A(0,k-1)
740: END DO
741: S = S + S
742: * double s for the off diagonal elements
743: CALL DLASSQ( K-1, A( 0+K*LDA ), LDA+1, SCALE, S )
744: * tri U at A(0,k)
745: CALL DLASSQ( K, A( 0+( K-1 )*LDA ), LDA+1, SCALE, S )
746: * tri L at A(0,k-1)
747: ELSE
1.6 ! bertrand 748: * A**T is lower
1.1 bertrand 749: DO J = 1, K - 1
750: CALL DLASSQ( J, A( 0+J*LDA ), 1, SCALE, S )
751: * U at A(0,0)
752: END DO
753: DO J = K, N - 1
754: CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
755: * k by k-1 rect. at A(0,k)
756: END DO
757: DO J = 0, K - 3
758: CALL DLASSQ( K-J-2, A( J+2+J*LDA ), 1, SCALE, S )
759: * L at A(1,0)
760: END DO
761: S = S + S
762: * double s for the off diagonal elements
763: CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S )
764: * tri U at A(0,0)
765: CALL DLASSQ( K-1, A( 1 ), LDA+1, SCALE, S )
766: * tri L at A(1,0)
767: END IF
768: END IF
769: ELSE
770: * n is even
771: IF( IFM.EQ.1 ) THEN
772: * A is normal
773: IF( ILU.EQ.0 ) THEN
774: * A is upper
775: DO J = 0, K - 2
776: CALL DLASSQ( K-J-1, A( K+J+2+J*LDA ), 1, SCALE, S )
777: * L at A(k+1,0)
778: END DO
779: DO J = 0, K - 1
780: CALL DLASSQ( K+J, A( 0+J*LDA ), 1, SCALE, S )
781: * trap U at A(0,0)
782: END DO
783: S = S + S
784: * double s for the off diagonal elements
785: CALL DLASSQ( K, A( K+1 ), LDA+1, SCALE, S )
786: * tri L at A(k+1,0)
787: CALL DLASSQ( K, A( K ), LDA+1, SCALE, S )
788: * tri U at A(k,0)
789: ELSE
790: * ilu=1 & A is lower
791: DO J = 0, K - 1
792: CALL DLASSQ( N-J-1, A( J+2+J*LDA ), 1, SCALE, S )
793: * trap L at A(1,0)
794: END DO
795: DO J = 1, K - 1
796: CALL DLASSQ( J, A( 0+J*LDA ), 1, SCALE, S )
797: * U at A(0,0)
798: END DO
799: S = S + S
800: * double s for the off diagonal elements
801: CALL DLASSQ( K, A( 1 ), LDA+1, SCALE, S )
802: * tri L at A(1,0)
803: CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S )
804: * tri U at A(0,0)
805: END IF
806: ELSE
807: * A is xpose
808: IF( ILU.EQ.0 ) THEN
1.6 ! bertrand 809: * A**T is upper
1.1 bertrand 810: DO J = 1, K - 1
811: CALL DLASSQ( J, A( 0+( K+1+J )*LDA ), 1, SCALE, S )
812: * U at A(0,k+1)
813: END DO
814: DO J = 0, K - 1
815: CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
816: * k by k rect. at A(0,0)
817: END DO
818: DO J = 0, K - 2
819: CALL DLASSQ( K-J-1, A( J+1+( J+K )*LDA ), 1, SCALE,
1.6 ! bertrand 820: $ S )
1.1 bertrand 821: * L at A(0,k)
822: END DO
823: S = S + S
824: * double s for the off diagonal elements
825: CALL DLASSQ( K, A( 0+( K+1 )*LDA ), LDA+1, SCALE, S )
826: * tri U at A(0,k+1)
827: CALL DLASSQ( K, A( 0+K*LDA ), LDA+1, SCALE, S )
828: * tri L at A(0,k)
829: ELSE
1.6 ! bertrand 830: * A**T is lower
1.1 bertrand 831: DO J = 1, K - 1
832: CALL DLASSQ( J, A( 0+( J+1 )*LDA ), 1, SCALE, S )
833: * U at A(0,1)
834: END DO
835: DO J = K + 1, N
836: CALL DLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
837: * k by k rect. at A(0,k+1)
838: END DO
839: DO J = 0, K - 2
840: CALL DLASSQ( K-J-1, A( J+1+J*LDA ), 1, SCALE, S )
841: * L at A(0,0)
842: END DO
843: S = S + S
844: * double s for the off diagonal elements
845: CALL DLASSQ( K, A( LDA ), LDA+1, SCALE, S )
846: * tri L at A(0,1)
847: CALL DLASSQ( K, A( 0 ), LDA+1, SCALE, S )
848: * tri U at A(0,0)
849: END IF
850: END IF
851: END IF
852: VALUE = SCALE*SQRT( S )
853: END IF
854: *
855: DLANSF = VALUE
856: RETURN
857: *
858: * End of DLANSF
859: *
860: END