1: SUBROUTINE DPSTRF( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
2: *
3: * -- LAPACK 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: * DPSTRF 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 3 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: * LDA (input) INTEGER
57: * The leading dimension of the array A. LDA >= max(1,N).
58: *
59: * PIV (output) INTEGER array, dimension (N)
60: * PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
61: *
62: * RANK (output) INTEGER
63: * The rank of A given by the number of steps the algorithm
64: * completed.
65: *
66: * TOL (input) DOUBLE PRECISION
67: * User defined tolerance. If TOL < 0, then N*U*MAX( A(K,K) )
68: * will be used. The algorithm terminates at the (K-1)st step
69: * if the pivot <= TOL.
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, JB, K, NB, PVT
90: LOGICAL UPPER
91: * ..
92: * .. External Functions ..
93: DOUBLE PRECISION DLAMCH
94: INTEGER ILAENV
95: LOGICAL LSAME, DISNAN
96: EXTERNAL DLAMCH, ILAENV, LSAME, DISNAN
97: * ..
98: * .. External Subroutines ..
99: EXTERNAL DGEMV, DPSTF2, DSCAL, DSWAP, DSYRK, XERBLA
100: * ..
101: * .. Intrinsic Functions ..
102: INTRINSIC MAX, MIN, SQRT, MAXLOC
103: * ..
104: * .. Executable Statements ..
105: *
106: * Test the input parameters.
107: *
108: INFO = 0
109: UPPER = LSAME( UPLO, 'U' )
110: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
111: INFO = -1
112: ELSE IF( N.LT.0 ) THEN
113: INFO = -2
114: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
115: INFO = -4
116: END IF
117: IF( INFO.NE.0 ) THEN
118: CALL XERBLA( 'DPSTRF', -INFO )
119: RETURN
120: END IF
121: *
122: * Quick return if possible
123: *
124: IF( N.EQ.0 )
125: $ RETURN
126: *
127: * Get block size
128: *
129: NB = ILAENV( 1, 'DPOTRF', UPLO, N, -1, -1, -1 )
130: IF( NB.LE.1 .OR. NB.GE.N ) THEN
131: *
132: * Use unblocked code
133: *
134: CALL DPSTF2( UPLO, N, A( 1, 1 ), LDA, PIV, RANK, TOL, WORK,
135: $ INFO )
136: GO TO 200
137: *
138: ELSE
139: *
140: * Initialize PIV
141: *
142: DO 100 I = 1, N
143: PIV( I ) = I
144: 100 CONTINUE
145: *
146: * Compute stopping value
147: *
148: PVT = 1
149: AJJ = A( PVT, PVT )
150: DO I = 2, N
151: IF( A( I, I ).GT.AJJ ) THEN
152: PVT = I
153: AJJ = A( PVT, PVT )
154: END IF
155: END DO
156: IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
157: RANK = 0
158: INFO = 1
159: GO TO 200
160: END IF
161: *
162: * Compute stopping value if not supplied
163: *
164: IF( TOL.LT.ZERO ) THEN
165: DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ
166: ELSE
167: DSTOP = TOL
168: END IF
169: *
170: *
171: IF( UPPER ) THEN
172: *
173: * Compute the Cholesky factorization P' * A * P = U' * U
174: *
175: DO 140 K = 1, N, NB
176: *
177: * Account for last block not being NB wide
178: *
179: JB = MIN( NB, N-K+1 )
180: *
181: * Set relevant part of first half of WORK to zero,
182: * holds dot products
183: *
184: DO 110 I = K, N
185: WORK( I ) = 0
186: 110 CONTINUE
187: *
188: DO 130 J = K, K + JB - 1
189: *
190: * Find pivot, test for exit, else swap rows and columns
191: * Update dot products, compute possible pivots which are
192: * stored in the second half of WORK
193: *
194: DO 120 I = J, N
195: *
196: IF( J.GT.K ) THEN
197: WORK( I ) = WORK( I ) + A( J-1, I )**2
198: END IF
199: WORK( N+I ) = A( I, I ) - WORK( I )
200: *
201: 120 CONTINUE
202: *
203: IF( J.GT.1 ) THEN
204: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
205: PVT = ITEMP + J - 1
206: AJJ = WORK( N+PVT )
207: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
208: A( J, J ) = AJJ
209: GO TO 190
210: END IF
211: END IF
212: *
213: IF( J.NE.PVT ) THEN
214: *
215: * Pivot OK, so can now swap pivot rows and columns
216: *
217: A( PVT, PVT ) = A( J, J )
218: CALL DSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
219: IF( PVT.LT.N )
220: $ CALL DSWAP( N-PVT, A( J, PVT+1 ), LDA,
221: $ A( PVT, PVT+1 ), LDA )
222: CALL DSWAP( PVT-J-1, A( J, J+1 ), LDA,
223: $ A( J+1, PVT ), 1 )
224: *
225: * Swap dot products and PIV
226: *
227: DTEMP = WORK( J )
228: WORK( J ) = WORK( PVT )
229: WORK( PVT ) = DTEMP
230: ITEMP = PIV( PVT )
231: PIV( PVT ) = PIV( J )
232: PIV( J ) = ITEMP
233: END IF
234: *
235: AJJ = SQRT( AJJ )
236: A( J, J ) = AJJ
237: *
238: * Compute elements J+1:N of row J.
239: *
240: IF( J.LT.N ) THEN
241: CALL DGEMV( 'Trans', J-K, N-J, -ONE, A( K, J+1 ),
242: $ LDA, A( K, J ), 1, ONE, A( J, J+1 ),
243: $ LDA )
244: CALL DSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
245: END IF
246: *
247: 130 CONTINUE
248: *
249: * Update trailing matrix, J already incremented
250: *
251: IF( K+JB.LE.N ) THEN
252: CALL DSYRK( 'Upper', 'Trans', N-J+1, JB, -ONE,
253: $ A( K, J ), LDA, ONE, A( J, J ), LDA )
254: END IF
255: *
256: 140 CONTINUE
257: *
258: ELSE
259: *
260: * Compute the Cholesky factorization P' * A * P = L * L'
261: *
262: DO 180 K = 1, N, NB
263: *
264: * Account for last block not being NB wide
265: *
266: JB = MIN( NB, N-K+1 )
267: *
268: * Set relevant part of first half of WORK to zero,
269: * holds dot products
270: *
271: DO 150 I = K, N
272: WORK( I ) = 0
273: 150 CONTINUE
274: *
275: DO 170 J = K, K + JB - 1
276: *
277: * Find pivot, test for exit, else swap rows and columns
278: * Update dot products, compute possible pivots which are
279: * stored in the second half of WORK
280: *
281: DO 160 I = J, N
282: *
283: IF( J.GT.K ) THEN
284: WORK( I ) = WORK( I ) + A( I, J-1 )**2
285: END IF
286: WORK( N+I ) = A( I, I ) - WORK( I )
287: *
288: 160 CONTINUE
289: *
290: IF( J.GT.1 ) THEN
291: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
292: PVT = ITEMP + J - 1
293: AJJ = WORK( N+PVT )
294: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
295: A( J, J ) = AJJ
296: GO TO 190
297: END IF
298: END IF
299: *
300: IF( J.NE.PVT ) THEN
301: *
302: * Pivot OK, so can now swap pivot rows and columns
303: *
304: A( PVT, PVT ) = A( J, J )
305: CALL DSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
306: IF( PVT.LT.N )
307: $ CALL DSWAP( N-PVT, A( PVT+1, J ), 1,
308: $ A( PVT+1, PVT ), 1 )
309: CALL DSWAP( PVT-J-1, A( J+1, J ), 1, A( PVT, J+1 ),
310: $ LDA )
311: *
312: * Swap dot products and PIV
313: *
314: DTEMP = WORK( J )
315: WORK( J ) = WORK( PVT )
316: WORK( PVT ) = DTEMP
317: ITEMP = PIV( PVT )
318: PIV( PVT ) = PIV( J )
319: PIV( J ) = ITEMP
320: END IF
321: *
322: AJJ = SQRT( AJJ )
323: A( J, J ) = AJJ
324: *
325: * Compute elements J+1:N of column J.
326: *
327: IF( J.LT.N ) THEN
328: CALL DGEMV( 'No Trans', N-J, J-K, -ONE,
329: $ A( J+1, K ), LDA, A( J, K ), LDA, ONE,
330: $ A( J+1, J ), 1 )
331: CALL DSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
332: END IF
333: *
334: 170 CONTINUE
335: *
336: * Update trailing matrix, J already incremented
337: *
338: IF( K+JB.LE.N ) THEN
339: CALL DSYRK( 'Lower', 'No Trans', N-J+1, JB, -ONE,
340: $ A( J, K ), LDA, ONE, A( J, J ), LDA )
341: END IF
342: *
343: 180 CONTINUE
344: *
345: END IF
346: END IF
347: *
348: * Ran to completion, A has full rank
349: *
350: RANK = N
351: *
352: GO TO 200
353: 190 CONTINUE
354: *
355: * Rank is the number of steps completed. Set INFO = 1 to signal
356: * that the factorization cannot be used to solve a system.
357: *
358: RANK = J - 1
359: INFO = 1
360: *
361: 200 CONTINUE
362: RETURN
363: *
364: * End of DPSTRF
365: *
366: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dpstrf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack