1: SUBROUTINE DPSTF2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
2: *
3: * -- LAPACK PROTOTYPE routine (version 3.2.2) --
4: * Craig Lucas, University of Manchester / NAG Ltd.
5: * October, 2008
6: *
7: * .. Scalar Arguments ..
8: DOUBLE PRECISION TOL
9: INTEGER INFO, LDA, N, RANK
10: CHARACTER UPLO
11: * ..
12: * .. Array Arguments ..
13: DOUBLE PRECISION A( LDA, * ), WORK( 2*N )
14: INTEGER PIV( N )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * DPSTF2 computes the Cholesky factorization with complete
21: * pivoting of a real symmetric positive semidefinite matrix A.
22: *
23: * The factorization has the form
24: * P' * A * P = U' * U , if UPLO = 'U',
25: * P' * A * P = L * L', if UPLO = 'L',
26: * where U is an upper triangular matrix and L is lower triangular, and
27: * P is stored as vector PIV.
28: *
29: * This algorithm does not attempt to check that A is positive
30: * semidefinite. This version of the algorithm calls level 2 BLAS.
31: *
32: * Arguments
33: * =========
34: *
35: * UPLO (input) CHARACTER*1
36: * Specifies whether the upper or lower triangular part of the
37: * symmetric matrix A is stored.
38: * = 'U': Upper triangular
39: * = 'L': Lower triangular
40: *
41: * N (input) INTEGER
42: * The order of the matrix A. N >= 0.
43: *
44: * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
45: * On entry, the symmetric matrix A. If UPLO = 'U', the leading
46: * n by n upper triangular part of A contains the upper
47: * triangular part of the matrix A, and the strictly lower
48: * triangular part of A is not referenced. If UPLO = 'L', the
49: * leading n by n lower triangular part of A contains the lower
50: * triangular part of the matrix A, and the strictly upper
51: * triangular part of A is not referenced.
52: *
53: * On exit, if INFO = 0, the factor U or L from the Cholesky
54: * factorization as above.
55: *
56: * PIV (output) INTEGER array, dimension (N)
57: * PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
58: *
59: * RANK (output) INTEGER
60: * The rank of A given by the number of steps the algorithm
61: * completed.
62: *
63: * TOL (input) DOUBLE PRECISION
64: * User defined tolerance. If TOL < 0, then N*U*MAX( A( K,K ) )
65: * will be used. The algorithm terminates at the (K-1)st step
66: * if the pivot <= TOL.
67: *
68: * LDA (input) INTEGER
69: * The leading dimension of the array A. LDA >= max(1,N).
70: *
71: * WORK (workspace) DOUBLE PRECISION array, dimension (2*N)
72: * Work space.
73: *
74: * INFO (output) INTEGER
75: * < 0: If INFO = -K, the K-th argument had an illegal value,
76: * = 0: algorithm completed successfully, and
77: * > 0: the matrix A is either rank deficient with computed rank
78: * as returned in RANK, or is indefinite. See Section 7 of
79: * LAPACK Working Note #161 for further information.
80: *
81: * =====================================================================
82: *
83: * .. Parameters ..
84: DOUBLE PRECISION ONE, ZERO
85: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
86: * ..
87: * .. Local Scalars ..
88: DOUBLE PRECISION AJJ, DSTOP, DTEMP
89: INTEGER I, ITEMP, J, PVT
90: LOGICAL UPPER
91: * ..
92: * .. External Functions ..
93: DOUBLE PRECISION DLAMCH
94: LOGICAL LSAME, DISNAN
95: EXTERNAL DLAMCH, LSAME, DISNAN
96: * ..
97: * .. External Subroutines ..
98: EXTERNAL DGEMV, DSCAL, DSWAP, XERBLA
99: * ..
100: * .. Intrinsic Functions ..
101: INTRINSIC MAX, SQRT, MAXLOC
102: * ..
103: * .. Executable Statements ..
104: *
105: * Test the input parameters
106: *
107: INFO = 0
108: UPPER = LSAME( UPLO, 'U' )
109: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
110: INFO = -1
111: ELSE IF( N.LT.0 ) THEN
112: INFO = -2
113: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
114: INFO = -4
115: END IF
116: IF( INFO.NE.0 ) THEN
117: CALL XERBLA( 'DPSTF2', -INFO )
118: RETURN
119: END IF
120: *
121: * Quick return if possible
122: *
123: IF( N.EQ.0 )
124: $ RETURN
125: *
126: * Initialize PIV
127: *
128: DO 100 I = 1, N
129: PIV( I ) = I
130: 100 CONTINUE
131: *
132: * Compute stopping value
133: *
134: PVT = 1
135: AJJ = A( PVT, PVT )
136: DO I = 2, N
137: IF( A( I, I ).GT.AJJ ) THEN
138: PVT = I
139: AJJ = A( PVT, PVT )
140: END IF
141: END DO
142: IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
143: RANK = 0
144: INFO = 1
145: GO TO 170
146: END IF
147: *
148: * Compute stopping value if not supplied
149: *
150: IF( TOL.LT.ZERO ) THEN
151: DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ
152: ELSE
153: DSTOP = TOL
154: END IF
155: *
156: * Set first half of WORK to zero, holds dot products
157: *
158: DO 110 I = 1, N
159: WORK( I ) = 0
160: 110 CONTINUE
161: *
162: IF( UPPER ) THEN
163: *
164: * Compute the Cholesky factorization P' * A * P = U' * U
165: *
166: DO 130 J = 1, N
167: *
168: * Find pivot, test for exit, else swap rows and columns
169: * Update dot products, compute possible pivots which are
170: * stored in the second half of WORK
171: *
172: DO 120 I = J, N
173: *
174: IF( J.GT.1 ) THEN
175: WORK( I ) = WORK( I ) + A( J-1, I )**2
176: END IF
177: WORK( N+I ) = A( I, I ) - WORK( I )
178: *
179: 120 CONTINUE
180: *
181: IF( J.GT.1 ) THEN
182: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
183: PVT = ITEMP + J - 1
184: AJJ = WORK( N+PVT )
185: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
186: A( J, J ) = AJJ
187: GO TO 160
188: END IF
189: END IF
190: *
191: IF( J.NE.PVT ) THEN
192: *
193: * Pivot OK, so can now swap pivot rows and columns
194: *
195: A( PVT, PVT ) = A( J, J )
196: CALL DSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
197: IF( PVT.LT.N )
198: $ CALL DSWAP( N-PVT, A( J, PVT+1 ), LDA,
199: $ A( PVT, PVT+1 ), LDA )
200: CALL DSWAP( PVT-J-1, A( J, J+1 ), LDA, A( J+1, PVT ), 1 )
201: *
202: * Swap dot products and PIV
203: *
204: DTEMP = WORK( J )
205: WORK( J ) = WORK( PVT )
206: WORK( PVT ) = DTEMP
207: ITEMP = PIV( PVT )
208: PIV( PVT ) = PIV( J )
209: PIV( J ) = ITEMP
210: END IF
211: *
212: AJJ = SQRT( AJJ )
213: A( J, J ) = AJJ
214: *
215: * Compute elements J+1:N of row J
216: *
217: IF( J.LT.N ) THEN
218: CALL DGEMV( 'Trans', J-1, N-J, -ONE, A( 1, J+1 ), LDA,
219: $ A( 1, J ), 1, ONE, A( J, J+1 ), LDA )
220: CALL DSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
221: END IF
222: *
223: 130 CONTINUE
224: *
225: ELSE
226: *
227: * Compute the Cholesky factorization P' * A * P = L * L'
228: *
229: DO 150 J = 1, N
230: *
231: * Find pivot, test for exit, else swap rows and columns
232: * Update dot products, compute possible pivots which are
233: * stored in the second half of WORK
234: *
235: DO 140 I = J, N
236: *
237: IF( J.GT.1 ) THEN
238: WORK( I ) = WORK( I ) + A( I, J-1 )**2
239: END IF
240: WORK( N+I ) = A( I, I ) - WORK( I )
241: *
242: 140 CONTINUE
243: *
244: IF( J.GT.1 ) THEN
245: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
246: PVT = ITEMP + J - 1
247: AJJ = WORK( N+PVT )
248: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
249: A( J, J ) = AJJ
250: GO TO 160
251: END IF
252: END IF
253: *
254: IF( J.NE.PVT ) THEN
255: *
256: * Pivot OK, so can now swap pivot rows and columns
257: *
258: A( PVT, PVT ) = A( J, J )
259: CALL DSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
260: IF( PVT.LT.N )
261: $ CALL DSWAP( N-PVT, A( PVT+1, J ), 1, A( PVT+1, PVT ),
262: $ 1 )
263: CALL DSWAP( PVT-J-1, A( J+1, J ), 1, A( PVT, J+1 ), LDA )
264: *
265: * Swap dot products and PIV
266: *
267: DTEMP = WORK( J )
268: WORK( J ) = WORK( PVT )
269: WORK( PVT ) = DTEMP
270: ITEMP = PIV( PVT )
271: PIV( PVT ) = PIV( J )
272: PIV( J ) = ITEMP
273: END IF
274: *
275: AJJ = SQRT( AJJ )
276: A( J, J ) = AJJ
277: *
278: * Compute elements J+1:N of column J
279: *
280: IF( J.LT.N ) THEN
281: CALL DGEMV( 'No Trans', N-J, J-1, -ONE, A( J+1, 1 ), LDA,
282: $ A( J, 1 ), LDA, ONE, A( J+1, J ), 1 )
283: CALL DSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
284: END IF
285: *
286: 150 CONTINUE
287: *
288: END IF
289: *
290: * Ran to completion, A has full rank
291: *
292: RANK = N
293: *
294: GO TO 170
295: 160 CONTINUE
296: *
297: * Rank is number of steps completed. Set INFO = 1 to signal
298: * that the factorization cannot be used to solve a system.
299: *
300: RANK = J - 1
301: INFO = 1
302: *
303: 170 CONTINUE
304: RETURN
305: *
306: * End of DPSTF2
307: *
308: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dpstf2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack