1: *> \brief \b ZPSTF2 computes the Cholesky factorization with complete pivoting of a complex Hermitian positive semidefinite matrix.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZPSTF2 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zpstf2.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zpstf2.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zpstf2.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZPSTF2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
22: *
23: * .. Scalar Arguments ..
24: * DOUBLE PRECISION TOL
25: * INTEGER INFO, LDA, N, RANK
26: * CHARACTER UPLO
27: * ..
28: * .. Array Arguments ..
29: * COMPLEX*16 A( LDA, * )
30: * DOUBLE PRECISION WORK( 2*N )
31: * INTEGER PIV( N )
32: * ..
33: *
34: *
35: *> \par Purpose:
36: * =============
37: *>
38: *> \verbatim
39: *>
40: *> ZPSTF2 computes the Cholesky factorization with complete
41: *> pivoting of a complex Hermitian positive semidefinite matrix A.
42: *>
43: *> The factorization has the form
44: *> P**T * A * P = U**H * U , if UPLO = 'U',
45: *> P**T * A * P = L * L**H, if UPLO = 'L',
46: *> where U is an upper triangular matrix and L is lower triangular, and
47: *> P is stored as vector PIV.
48: *>
49: *> This algorithm does not attempt to check that A is positive
50: *> semidefinite. This version of the algorithm calls level 2 BLAS.
51: *> \endverbatim
52: *
53: * Arguments:
54: * ==========
55: *
56: *> \param[in] UPLO
57: *> \verbatim
58: *> UPLO is CHARACTER*1
59: *> Specifies whether the upper or lower triangular part of the
60: *> symmetric matrix A is stored.
61: *> = 'U': Upper triangular
62: *> = 'L': Lower triangular
63: *> \endverbatim
64: *>
65: *> \param[in] N
66: *> \verbatim
67: *> N is INTEGER
68: *> The order of the matrix A. N >= 0.
69: *> \endverbatim
70: *>
71: *> \param[in,out] A
72: *> \verbatim
73: *> A is COMPLEX*16 array, dimension (LDA,N)
74: *> On entry, the symmetric matrix A. If UPLO = 'U', the leading
75: *> n by n upper triangular part of A contains the upper
76: *> triangular part of the matrix A, and the strictly lower
77: *> triangular part of A is not referenced. If UPLO = 'L', the
78: *> leading n by n lower triangular part of A contains the lower
79: *> triangular part of the matrix A, and the strictly upper
80: *> triangular part of A is not referenced.
81: *>
82: *> On exit, if INFO = 0, the factor U or L from the Cholesky
83: *> factorization as above.
84: *> \endverbatim
85: *>
86: *> \param[out] PIV
87: *> \verbatim
88: *> PIV is INTEGER array, dimension (N)
89: *> PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
90: *> \endverbatim
91: *>
92: *> \param[out] RANK
93: *> \verbatim
94: *> RANK is INTEGER
95: *> The rank of A given by the number of steps the algorithm
96: *> completed.
97: *> \endverbatim
98: *>
99: *> \param[in] TOL
100: *> \verbatim
101: *> TOL is DOUBLE PRECISION
102: *> User defined tolerance. If TOL < 0, then N*U*MAX( A( K,K ) )
103: *> will be used. The algorithm terminates at the (K-1)st step
104: *> if the pivot <= TOL.
105: *> \endverbatim
106: *>
107: *> \param[in] LDA
108: *> \verbatim
109: *> LDA is INTEGER
110: *> The leading dimension of the array A. LDA >= max(1,N).
111: *> \endverbatim
112: *>
113: *> \param[out] WORK
114: *> \verbatim
115: *> WORK is DOUBLE PRECISION array, dimension (2*N)
116: *> Work space.
117: *> \endverbatim
118: *>
119: *> \param[out] INFO
120: *> \verbatim
121: *> INFO is INTEGER
122: *> < 0: If INFO = -K, the K-th argument had an illegal value,
123: *> = 0: algorithm completed successfully, and
124: *> > 0: the matrix A is either rank deficient with computed rank
125: *> as returned in RANK, or is not positive semidefinite. See
126: *> Section 7 of LAPACK Working Note #161 for further
127: *> information.
128: *> \endverbatim
129: *
130: * Authors:
131: * ========
132: *
133: *> \author Univ. of Tennessee
134: *> \author Univ. of California Berkeley
135: *> \author Univ. of Colorado Denver
136: *> \author NAG Ltd.
137: *
138: *> \ingroup complex16OTHERcomputational
139: *
140: * =====================================================================
141: SUBROUTINE ZPSTF2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
142: *
143: * -- LAPACK computational routine --
144: * -- LAPACK is a software package provided by Univ. of Tennessee, --
145: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
146: *
147: * .. Scalar Arguments ..
148: DOUBLE PRECISION TOL
149: INTEGER INFO, LDA, N, RANK
150: CHARACTER UPLO
151: * ..
152: * .. Array Arguments ..
153: COMPLEX*16 A( LDA, * )
154: DOUBLE PRECISION WORK( 2*N )
155: INTEGER PIV( N )
156: * ..
157: *
158: * =====================================================================
159: *
160: * .. Parameters ..
161: DOUBLE PRECISION ONE, ZERO
162: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
163: COMPLEX*16 CONE
164: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
165: * ..
166: * .. Local Scalars ..
167: COMPLEX*16 ZTEMP
168: DOUBLE PRECISION AJJ, DSTOP, DTEMP
169: INTEGER I, ITEMP, J, PVT
170: LOGICAL UPPER
171: * ..
172: * .. External Functions ..
173: DOUBLE PRECISION DLAMCH
174: LOGICAL LSAME, DISNAN
175: EXTERNAL DLAMCH, LSAME, DISNAN
176: * ..
177: * .. External Subroutines ..
178: EXTERNAL ZDSCAL, ZGEMV, ZLACGV, ZSWAP, XERBLA
179: * ..
180: * .. Intrinsic Functions ..
181: INTRINSIC DBLE, DCONJG, MAX, SQRT
182: * ..
183: * .. Executable Statements ..
184: *
185: * Test the input parameters
186: *
187: INFO = 0
188: UPPER = LSAME( UPLO, 'U' )
189: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
190: INFO = -1
191: ELSE IF( N.LT.0 ) THEN
192: INFO = -2
193: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
194: INFO = -4
195: END IF
196: IF( INFO.NE.0 ) THEN
197: CALL XERBLA( 'ZPSTF2', -INFO )
198: RETURN
199: END IF
200: *
201: * Quick return if possible
202: *
203: IF( N.EQ.0 )
204: $ RETURN
205: *
206: * Initialize PIV
207: *
208: DO 100 I = 1, N
209: PIV( I ) = I
210: 100 CONTINUE
211: *
212: * Compute stopping value
213: *
214: DO 110 I = 1, N
215: WORK( I ) = DBLE( A( I, I ) )
216: 110 CONTINUE
217: PVT = MAXLOC( WORK( 1:N ), 1 )
218: AJJ = DBLE( A( PVT, PVT ) )
219: IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN
220: RANK = 0
221: INFO = 1
222: GO TO 200
223: END IF
224: *
225: * Compute stopping value if not supplied
226: *
227: IF( TOL.LT.ZERO ) THEN
228: DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ
229: ELSE
230: DSTOP = TOL
231: END IF
232: *
233: * Set first half of WORK to zero, holds dot products
234: *
235: DO 120 I = 1, N
236: WORK( I ) = 0
237: 120 CONTINUE
238: *
239: IF( UPPER ) THEN
240: *
241: * Compute the Cholesky factorization P**T * A * P = U**H* U
242: *
243: DO 150 J = 1, N
244: *
245: * Find pivot, test for exit, else swap rows and columns
246: * Update dot products, compute possible pivots which are
247: * stored in the second half of WORK
248: *
249: DO 130 I = J, N
250: *
251: IF( J.GT.1 ) THEN
252: WORK( I ) = WORK( I ) +
253: $ DBLE( DCONJG( A( J-1, I ) )*
254: $ A( J-1, I ) )
255: END IF
256: WORK( N+I ) = DBLE( A( I, I ) ) - WORK( I )
257: *
258: 130 CONTINUE
259: *
260: IF( J.GT.1 ) THEN
261: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
262: PVT = ITEMP + J - 1
263: AJJ = WORK( N+PVT )
264: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
265: A( J, J ) = AJJ
266: GO TO 190
267: END IF
268: END IF
269: *
270: IF( J.NE.PVT ) THEN
271: *
272: * Pivot OK, so can now swap pivot rows and columns
273: *
274: A( PVT, PVT ) = A( J, J )
275: CALL ZSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
276: IF( PVT.LT.N )
277: $ CALL ZSWAP( N-PVT, A( J, PVT+1 ), LDA,
278: $ A( PVT, PVT+1 ), LDA )
279: DO 140 I = J + 1, PVT - 1
280: ZTEMP = DCONJG( A( J, I ) )
281: A( J, I ) = DCONJG( A( I, PVT ) )
282: A( I, PVT ) = ZTEMP
283: 140 CONTINUE
284: A( J, PVT ) = DCONJG( A( J, PVT ) )
285: *
286: * Swap dot products and PIV
287: *
288: DTEMP = WORK( J )
289: WORK( J ) = WORK( PVT )
290: WORK( PVT ) = DTEMP
291: ITEMP = PIV( PVT )
292: PIV( PVT ) = PIV( J )
293: PIV( J ) = ITEMP
294: END IF
295: *
296: AJJ = SQRT( AJJ )
297: A( J, J ) = AJJ
298: *
299: * Compute elements J+1:N of row J
300: *
301: IF( J.LT.N ) THEN
302: CALL ZLACGV( J-1, A( 1, J ), 1 )
303: CALL ZGEMV( 'Trans', J-1, N-J, -CONE, A( 1, J+1 ), LDA,
304: $ A( 1, J ), 1, CONE, A( J, J+1 ), LDA )
305: CALL ZLACGV( J-1, A( 1, J ), 1 )
306: CALL ZDSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
307: END IF
308: *
309: 150 CONTINUE
310: *
311: ELSE
312: *
313: * Compute the Cholesky factorization P**T * A * P = L * L**H
314: *
315: DO 180 J = 1, N
316: *
317: * Find pivot, test for exit, else swap rows and columns
318: * Update dot products, compute possible pivots which are
319: * stored in the second half of WORK
320: *
321: DO 160 I = J, N
322: *
323: IF( J.GT.1 ) THEN
324: WORK( I ) = WORK( I ) +
325: $ DBLE( DCONJG( A( I, J-1 ) )*
326: $ A( I, J-1 ) )
327: END IF
328: WORK( N+I ) = DBLE( A( I, I ) ) - WORK( I )
329: *
330: 160 CONTINUE
331: *
332: IF( J.GT.1 ) THEN
333: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
334: PVT = ITEMP + J - 1
335: AJJ = WORK( N+PVT )
336: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
337: A( J, J ) = AJJ
338: GO TO 190
339: END IF
340: END IF
341: *
342: IF( J.NE.PVT ) THEN
343: *
344: * Pivot OK, so can now swap pivot rows and columns
345: *
346: A( PVT, PVT ) = A( J, J )
347: CALL ZSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
348: IF( PVT.LT.N )
349: $ CALL ZSWAP( N-PVT, A( PVT+1, J ), 1, A( PVT+1, PVT ),
350: $ 1 )
351: DO 170 I = J + 1, PVT - 1
352: ZTEMP = DCONJG( A( I, J ) )
353: A( I, J ) = DCONJG( A( PVT, I ) )
354: A( PVT, I ) = ZTEMP
355: 170 CONTINUE
356: A( PVT, J ) = DCONJG( A( PVT, J ) )
357: *
358: * Swap dot products and PIV
359: *
360: DTEMP = WORK( J )
361: WORK( J ) = WORK( PVT )
362: WORK( PVT ) = DTEMP
363: ITEMP = PIV( PVT )
364: PIV( PVT ) = PIV( J )
365: PIV( J ) = ITEMP
366: END IF
367: *
368: AJJ = SQRT( AJJ )
369: A( J, J ) = AJJ
370: *
371: * Compute elements J+1:N of column J
372: *
373: IF( J.LT.N ) THEN
374: CALL ZLACGV( J-1, A( J, 1 ), LDA )
375: CALL ZGEMV( 'No Trans', N-J, J-1, -CONE, A( J+1, 1 ),
376: $ LDA, A( J, 1 ), LDA, CONE, A( J+1, J ), 1 )
377: CALL ZLACGV( J-1, A( J, 1 ), LDA )
378: CALL ZDSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
379: END IF
380: *
381: 180 CONTINUE
382: *
383: END IF
384: *
385: * Ran to completion, A has full rank
386: *
387: RANK = N
388: *
389: GO TO 200
390: 190 CONTINUE
391: *
392: * Rank is number of steps completed. Set INFO = 1 to signal
393: * that the factorization cannot be used to solve a system.
394: *
395: RANK = J - 1
396: INFO = 1
397: *
398: 200 CONTINUE
399: RETURN
400: *
401: * End of ZPSTF2
402: *
403: 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