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