1: SUBROUTINE DSYTF2( UPLO, N, A, LDA, IPIV, INFO )
2: *
3: * -- LAPACK routine (version 3.2) --
4: * -- LAPACK is a software package provided by Univ. of Tennessee, --
5: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6: * November 2006
7: *
8: * .. Scalar Arguments ..
9: CHARACTER UPLO
10: INTEGER INFO, LDA, N
11: * ..
12: * .. Array Arguments ..
13: INTEGER IPIV( * )
14: DOUBLE PRECISION A( LDA, * )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * DSYTF2 computes the factorization of a real symmetric matrix A using
21: * the Bunch-Kaufman diagonal pivoting method:
22: *
23: * A = U*D*U' or A = L*D*L'
24: *
25: * where U (or L) is a product of permutation and unit upper (lower)
26: * triangular matrices, U' is the transpose of U, and D is symmetric and
27: * block diagonal with 1-by-1 and 2-by-2 diagonal blocks.
28: *
29: * This is the unblocked version of the algorithm, calling Level 2 BLAS.
30: *
31: * Arguments
32: * =========
33: *
34: * UPLO (input) CHARACTER*1
35: * Specifies whether the upper or lower triangular part of the
36: * symmetric matrix A is stored:
37: * = 'U': Upper triangular
38: * = 'L': Lower triangular
39: *
40: * N (input) INTEGER
41: * The order of the matrix A. N >= 0.
42: *
43: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
44: * On entry, the symmetric matrix A. If UPLO = 'U', the leading
45: * n-by-n upper triangular part of A contains the upper
46: * triangular part of the matrix A, and the strictly lower
47: * triangular part of A is not referenced. If UPLO = 'L', the
48: * leading n-by-n lower triangular part of A contains the lower
49: * triangular part of the matrix A, and the strictly upper
50: * triangular part of A is not referenced.
51: *
52: * On exit, the block diagonal matrix D and the multipliers used
53: * to obtain the factor U or L (see below for further details).
54: *
55: * LDA (input) INTEGER
56: * The leading dimension of the array A. LDA >= max(1,N).
57: *
58: * IPIV (output) INTEGER array, dimension (N)
59: * Details of the interchanges and the block structure of D.
60: * If IPIV(k) > 0, then rows and columns k and IPIV(k) were
61: * interchanged and D(k,k) is a 1-by-1 diagonal block.
62: * If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
63: * columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
64: * is a 2-by-2 diagonal block. If UPLO = 'L' and IPIV(k) =
65: * IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
66: * interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
67: *
68: * INFO (output) INTEGER
69: * = 0: successful exit
70: * < 0: if INFO = -k, the k-th argument had an illegal value
71: * > 0: if INFO = k, D(k,k) is exactly zero. The factorization
72: * has been completed, but the block diagonal matrix D is
73: * exactly singular, and division by zero will occur if it
74: * is used to solve a system of equations.
75: *
76: * Further Details
77: * ===============
78: *
79: * 09-29-06 - patch from
80: * Bobby Cheng, MathWorks
81: *
82: * Replace l.204 and l.372
83: * IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
84: * by
85: * IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
86: *
87: * 01-01-96 - Based on modifications by
88: * J. Lewis, Boeing Computer Services Company
89: * A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
90: * 1-96 - Based on modifications by J. Lewis, Boeing Computer Services
91: * Company
92: *
93: * If UPLO = 'U', then A = U*D*U', where
94: * U = P(n)*U(n)* ... *P(k)U(k)* ...,
95: * i.e., U is a product of terms P(k)*U(k), where k decreases from n to
96: * 1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
97: * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
98: * defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
99: * that if the diagonal block D(k) is of order s (s = 1 or 2), then
100: *
101: * ( I v 0 ) k-s
102: * U(k) = ( 0 I 0 ) s
103: * ( 0 0 I ) n-k
104: * k-s s n-k
105: *
106: * If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
107: * If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
108: * and A(k,k), and v overwrites A(1:k-2,k-1:k).
109: *
110: * If UPLO = 'L', then A = L*D*L', where
111: * L = P(1)*L(1)* ... *P(k)*L(k)* ...,
112: * i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
113: * n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
114: * and 2-by-2 diagonal blocks D(k). P(k) is a permutation matrix as
115: * defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
116: * that if the diagonal block D(k) is of order s (s = 1 or 2), then
117: *
118: * ( I 0 0 ) k-1
119: * L(k) = ( 0 I 0 ) s
120: * ( 0 v I ) n-k-s+1
121: * k-1 s n-k-s+1
122: *
123: * If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
124: * If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
125: * and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
126: *
127: * =====================================================================
128: *
129: * .. Parameters ..
130: DOUBLE PRECISION ZERO, ONE
131: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
132: DOUBLE PRECISION EIGHT, SEVTEN
133: PARAMETER ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
134: * ..
135: * .. Local Scalars ..
136: LOGICAL UPPER
137: INTEGER I, IMAX, J, JMAX, K, KK, KP, KSTEP
138: DOUBLE PRECISION ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, R1,
139: $ ROWMAX, T, WK, WKM1, WKP1
140: * ..
141: * .. External Functions ..
142: LOGICAL LSAME, DISNAN
143: INTEGER IDAMAX
144: EXTERNAL LSAME, IDAMAX, DISNAN
145: * ..
146: * .. External Subroutines ..
147: EXTERNAL DSCAL, DSWAP, DSYR, XERBLA
148: * ..
149: * .. Intrinsic Functions ..
150: INTRINSIC ABS, MAX, SQRT
151: * ..
152: * .. Executable Statements ..
153: *
154: * Test the input parameters.
155: *
156: INFO = 0
157: UPPER = LSAME( UPLO, 'U' )
158: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
159: INFO = -1
160: ELSE IF( N.LT.0 ) THEN
161: INFO = -2
162: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
163: INFO = -4
164: END IF
165: IF( INFO.NE.0 ) THEN
166: CALL XERBLA( 'DSYTF2', -INFO )
167: RETURN
168: END IF
169: *
170: * Initialize ALPHA for use in choosing pivot block size.
171: *
172: ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
173: *
174: IF( UPPER ) THEN
175: *
176: * Factorize A as U*D*U' using the upper triangle of A
177: *
178: * K is the main loop index, decreasing from N to 1 in steps of
179: * 1 or 2
180: *
181: K = N
182: 10 CONTINUE
183: *
184: * If K < 1, exit from loop
185: *
186: IF( K.LT.1 )
187: $ GO TO 70
188: KSTEP = 1
189: *
190: * Determine rows and columns to be interchanged and whether
191: * a 1-by-1 or 2-by-2 pivot block will be used
192: *
193: ABSAKK = ABS( A( K, K ) )
194: *
195: * IMAX is the row-index of the largest off-diagonal element in
196: * column K, and COLMAX is its absolute value
197: *
198: IF( K.GT.1 ) THEN
199: IMAX = IDAMAX( K-1, A( 1, K ), 1 )
200: COLMAX = ABS( A( IMAX, K ) )
201: ELSE
202: COLMAX = ZERO
203: END IF
204: *
205: IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
206: *
207: * Column K is zero or contains a NaN: set INFO and continue
208: *
209: IF( INFO.EQ.0 )
210: $ INFO = K
211: KP = K
212: ELSE
213: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
214: *
215: * no interchange, use 1-by-1 pivot block
216: *
217: KP = K
218: ELSE
219: *
220: * JMAX is the column-index of the largest off-diagonal
221: * element in row IMAX, and ROWMAX is its absolute value
222: *
223: JMAX = IMAX + IDAMAX( K-IMAX, A( IMAX, IMAX+1 ), LDA )
224: ROWMAX = ABS( A( IMAX, JMAX ) )
225: IF( IMAX.GT.1 ) THEN
226: JMAX = IDAMAX( IMAX-1, A( 1, IMAX ), 1 )
227: ROWMAX = MAX( ROWMAX, ABS( A( JMAX, IMAX ) ) )
228: END IF
229: *
230: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
231: *
232: * no interchange, use 1-by-1 pivot block
233: *
234: KP = K
235: ELSE IF( ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN
236: *
237: * interchange rows and columns K and IMAX, use 1-by-1
238: * pivot block
239: *
240: KP = IMAX
241: ELSE
242: *
243: * interchange rows and columns K-1 and IMAX, use 2-by-2
244: * pivot block
245: *
246: KP = IMAX
247: KSTEP = 2
248: END IF
249: END IF
250: *
251: KK = K - KSTEP + 1
252: IF( KP.NE.KK ) THEN
253: *
254: * Interchange rows and columns KK and KP in the leading
255: * submatrix A(1:k,1:k)
256: *
257: CALL DSWAP( KP-1, A( 1, KK ), 1, A( 1, KP ), 1 )
258: CALL DSWAP( KK-KP-1, A( KP+1, KK ), 1, A( KP, KP+1 ),
259: $ LDA )
260: T = A( KK, KK )
261: A( KK, KK ) = A( KP, KP )
262: A( KP, KP ) = T
263: IF( KSTEP.EQ.2 ) THEN
264: T = A( K-1, K )
265: A( K-1, K ) = A( KP, K )
266: A( KP, K ) = T
267: END IF
268: END IF
269: *
270: * Update the leading submatrix
271: *
272: IF( KSTEP.EQ.1 ) THEN
273: *
274: * 1-by-1 pivot block D(k): column k now holds
275: *
276: * W(k) = U(k)*D(k)
277: *
278: * where U(k) is the k-th column of U
279: *
280: * Perform a rank-1 update of A(1:k-1,1:k-1) as
281: *
282: * A := A - U(k)*D(k)*U(k)' = A - W(k)*1/D(k)*W(k)'
283: *
284: R1 = ONE / A( K, K )
285: CALL DSYR( UPLO, K-1, -R1, A( 1, K ), 1, A, LDA )
286: *
287: * Store U(k) in column k
288: *
289: CALL DSCAL( K-1, R1, A( 1, K ), 1 )
290: ELSE
291: *
292: * 2-by-2 pivot block D(k): columns k and k-1 now hold
293: *
294: * ( W(k-1) W(k) ) = ( U(k-1) U(k) )*D(k)
295: *
296: * where U(k) and U(k-1) are the k-th and (k-1)-th columns
297: * of U
298: *
299: * Perform a rank-2 update of A(1:k-2,1:k-2) as
300: *
301: * A := A - ( U(k-1) U(k) )*D(k)*( U(k-1) U(k) )'
302: * = A - ( W(k-1) W(k) )*inv(D(k))*( W(k-1) W(k) )'
303: *
304: IF( K.GT.2 ) THEN
305: *
306: D12 = A( K-1, K )
307: D22 = A( K-1, K-1 ) / D12
308: D11 = A( K, K ) / D12
309: T = ONE / ( D11*D22-ONE )
310: D12 = T / D12
311: *
312: DO 30 J = K - 2, 1, -1
313: WKM1 = D12*( D11*A( J, K-1 )-A( J, K ) )
314: WK = D12*( D22*A( J, K )-A( J, K-1 ) )
315: DO 20 I = J, 1, -1
316: A( I, J ) = A( I, J ) - A( I, K )*WK -
317: $ A( I, K-1 )*WKM1
318: 20 CONTINUE
319: A( J, K ) = WK
320: A( J, K-1 ) = WKM1
321: 30 CONTINUE
322: *
323: END IF
324: *
325: END IF
326: END IF
327: *
328: * Store details of the interchanges in IPIV
329: *
330: IF( KSTEP.EQ.1 ) THEN
331: IPIV( K ) = KP
332: ELSE
333: IPIV( K ) = -KP
334: IPIV( K-1 ) = -KP
335: END IF
336: *
337: * Decrease K and return to the start of the main loop
338: *
339: K = K - KSTEP
340: GO TO 10
341: *
342: ELSE
343: *
344: * Factorize A as L*D*L' using the lower triangle of A
345: *
346: * K is the main loop index, increasing from 1 to N in steps of
347: * 1 or 2
348: *
349: K = 1
350: 40 CONTINUE
351: *
352: * If K > N, exit from loop
353: *
354: IF( K.GT.N )
355: $ GO TO 70
356: KSTEP = 1
357: *
358: * Determine rows and columns to be interchanged and whether
359: * a 1-by-1 or 2-by-2 pivot block will be used
360: *
361: ABSAKK = ABS( A( K, K ) )
362: *
363: * IMAX is the row-index of the largest off-diagonal element in
364: * column K, and COLMAX is its absolute value
365: *
366: IF( K.LT.N ) THEN
367: IMAX = K + IDAMAX( N-K, A( K+1, K ), 1 )
368: COLMAX = ABS( A( IMAX, K ) )
369: ELSE
370: COLMAX = ZERO
371: END IF
372: *
373: IF( (MAX( ABSAKK, COLMAX ).EQ.ZERO) .OR. DISNAN(ABSAKK) ) THEN
374: *
375: * Column K is zero or contains a NaN: set INFO and continue
376: *
377: IF( INFO.EQ.0 )
378: $ INFO = K
379: KP = K
380: ELSE
381: IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
382: *
383: * no interchange, use 1-by-1 pivot block
384: *
385: KP = K
386: ELSE
387: *
388: * JMAX is the column-index of the largest off-diagonal
389: * element in row IMAX, and ROWMAX is its absolute value
390: *
391: JMAX = K - 1 + IDAMAX( IMAX-K, A( IMAX, K ), LDA )
392: ROWMAX = ABS( A( IMAX, JMAX ) )
393: IF( IMAX.LT.N ) THEN
394: JMAX = IMAX + IDAMAX( N-IMAX, A( IMAX+1, IMAX ), 1 )
395: ROWMAX = MAX( ROWMAX, ABS( A( JMAX, IMAX ) ) )
396: END IF
397: *
398: IF( ABSAKK.GE.ALPHA*COLMAX*( COLMAX / ROWMAX ) ) THEN
399: *
400: * no interchange, use 1-by-1 pivot block
401: *
402: KP = K
403: ELSE IF( ABS( A( IMAX, IMAX ) ).GE.ALPHA*ROWMAX ) THEN
404: *
405: * interchange rows and columns K and IMAX, use 1-by-1
406: * pivot block
407: *
408: KP = IMAX
409: ELSE
410: *
411: * interchange rows and columns K+1 and IMAX, use 2-by-2
412: * pivot block
413: *
414: KP = IMAX
415: KSTEP = 2
416: END IF
417: END IF
418: *
419: KK = K + KSTEP - 1
420: IF( KP.NE.KK ) THEN
421: *
422: * Interchange rows and columns KK and KP in the trailing
423: * submatrix A(k:n,k:n)
424: *
425: IF( KP.LT.N )
426: $ CALL DSWAP( N-KP, A( KP+1, KK ), 1, A( KP+1, KP ), 1 )
427: CALL DSWAP( KP-KK-1, A( KK+1, KK ), 1, A( KP, KK+1 ),
428: $ LDA )
429: T = A( KK, KK )
430: A( KK, KK ) = A( KP, KP )
431: A( KP, KP ) = T
432: IF( KSTEP.EQ.2 ) THEN
433: T = A( K+1, K )
434: A( K+1, K ) = A( KP, K )
435: A( KP, K ) = T
436: END IF
437: END IF
438: *
439: * Update the trailing submatrix
440: *
441: IF( KSTEP.EQ.1 ) THEN
442: *
443: * 1-by-1 pivot block D(k): column k now holds
444: *
445: * W(k) = L(k)*D(k)
446: *
447: * where L(k) is the k-th column of L
448: *
449: IF( K.LT.N ) THEN
450: *
451: * Perform a rank-1 update of A(k+1:n,k+1:n) as
452: *
453: * A := A - L(k)*D(k)*L(k)' = A - W(k)*(1/D(k))*W(k)'
454: *
455: D11 = ONE / A( K, K )
456: CALL DSYR( UPLO, N-K, -D11, A( K+1, K ), 1,
457: $ A( K+1, K+1 ), LDA )
458: *
459: * Store L(k) in column K
460: *
461: CALL DSCAL( N-K, D11, A( K+1, K ), 1 )
462: END IF
463: ELSE
464: *
465: * 2-by-2 pivot block D(k)
466: *
467: IF( K.LT.N-1 ) THEN
468: *
469: * Perform a rank-2 update of A(k+2:n,k+2:n) as
470: *
471: * A := A - ( (A(k) A(k+1))*D(k)**(-1) ) * (A(k) A(k+1))'
472: *
473: * where L(k) and L(k+1) are the k-th and (k+1)-th
474: * columns of L
475: *
476: D21 = A( K+1, K )
477: D11 = A( K+1, K+1 ) / D21
478: D22 = A( K, K ) / D21
479: T = ONE / ( D11*D22-ONE )
480: D21 = T / D21
481: *
482: DO 60 J = K + 2, N
483: *
484: WK = D21*( D11*A( J, K )-A( J, K+1 ) )
485: WKP1 = D21*( D22*A( J, K+1 )-A( J, K ) )
486: *
487: DO 50 I = J, N
488: A( I, J ) = A( I, J ) - A( I, K )*WK -
489: $ A( I, K+1 )*WKP1
490: 50 CONTINUE
491: *
492: A( J, K ) = WK
493: A( J, K+1 ) = WKP1
494: *
495: 60 CONTINUE
496: END IF
497: END IF
498: END IF
499: *
500: * Store details of the interchanges in IPIV
501: *
502: IF( KSTEP.EQ.1 ) THEN
503: IPIV( K ) = KP
504: ELSE
505: IPIV( K ) = -KP
506: IPIV( K+1 ) = -KP
507: END IF
508: *
509: * Increase K and return to the start of the main loop
510: *
511: K = K + KSTEP
512: GO TO 40
513: *
514: END IF
515: *
516: 70 CONTINUE
517: *
518: RETURN
519: *
520: * End of DSYTF2
521: *
522: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dsytf2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack