Annotation of rpl/lapack/lapack/zpstf2.f, revision 1.12
1.12 ! bertrand 1: *> \brief \b ZPSTF2 computes the Cholesky factorization with complete pivoting of a complex Hermitian positive semidefinite matrix.
1.6 bertrand 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
1.12 ! bertrand 125: *> as returned in RANK, or is not positive semidefinite. See
! 126: *> Section 7 of LAPACK Working Note #161 for further
! 127: *> information.
1.6 bertrand 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: *
1.12 ! bertrand 138: *> \date November 2015
1.6 bertrand 139: *
140: *> \ingroup complex16OTHERcomputational
141: *
142: * =====================================================================
1.1 bertrand 143: SUBROUTINE ZPSTF2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
144: *
1.12 ! bertrand 145: * -- LAPACK computational routine (version 3.6.0) --
1.6 bertrand 146: * -- LAPACK is a software package provided by Univ. of Tennessee, --
147: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.12 ! bertrand 148: * November 2015
1.1 bertrand 149: *
150: * .. Scalar Arguments ..
151: DOUBLE PRECISION TOL
152: INTEGER INFO, LDA, N, RANK
153: CHARACTER UPLO
154: * ..
155: * .. Array Arguments ..
156: COMPLEX*16 A( LDA, * )
157: DOUBLE PRECISION WORK( 2*N )
158: INTEGER PIV( N )
159: * ..
160: *
161: * =====================================================================
162: *
163: * .. Parameters ..
164: DOUBLE PRECISION ONE, ZERO
165: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
166: COMPLEX*16 CONE
167: PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
168: * ..
169: * .. Local Scalars ..
170: COMPLEX*16 ZTEMP
171: DOUBLE PRECISION AJJ, DSTOP, DTEMP
172: INTEGER I, ITEMP, J, PVT
173: LOGICAL UPPER
174: * ..
175: * .. External Functions ..
176: DOUBLE PRECISION DLAMCH
177: LOGICAL LSAME, DISNAN
178: EXTERNAL DLAMCH, LSAME, DISNAN
179: * ..
180: * .. External Subroutines ..
181: EXTERNAL ZDSCAL, ZGEMV, ZLACGV, ZSWAP, XERBLA
182: * ..
183: * .. Intrinsic Functions ..
184: INTRINSIC DBLE, DCONJG, MAX, SQRT
185: * ..
186: * .. Executable Statements ..
187: *
188: * Test the input parameters
189: *
190: INFO = 0
191: UPPER = LSAME( UPLO, 'U' )
192: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
193: INFO = -1
194: ELSE IF( N.LT.0 ) THEN
195: INFO = -2
196: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
197: INFO = -4
198: END IF
199: IF( INFO.NE.0 ) THEN
200: CALL XERBLA( 'ZPSTF2', -INFO )
201: RETURN
202: END IF
203: *
204: * Quick return if possible
205: *
206: IF( N.EQ.0 )
207: $ RETURN
208: *
209: * Initialize PIV
210: *
211: DO 100 I = 1, N
212: PIV( I ) = I
213: 100 CONTINUE
214: *
215: * Compute stopping value
216: *
217: DO 110 I = 1, N
218: WORK( I ) = DBLE( A( I, I ) )
219: 110 CONTINUE
220: PVT = MAXLOC( WORK( 1:N ), 1 )
221: AJJ = DBLE( A( PVT, PVT ) )
1.12 ! bertrand 222: IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN
1.1 bertrand 223: RANK = 0
224: INFO = 1
225: GO TO 200
226: END IF
227: *
228: * Compute stopping value if not supplied
229: *
230: IF( TOL.LT.ZERO ) THEN
231: DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ
232: ELSE
233: DSTOP = TOL
234: END IF
235: *
236: * Set first half of WORK to zero, holds dot products
237: *
238: DO 120 I = 1, N
239: WORK( I ) = 0
240: 120 CONTINUE
241: *
242: IF( UPPER ) THEN
243: *
1.5 bertrand 244: * Compute the Cholesky factorization P**T * A * P = U**H* U
1.1 bertrand 245: *
246: DO 150 J = 1, N
247: *
248: * Find pivot, test for exit, else swap rows and columns
249: * Update dot products, compute possible pivots which are
250: * stored in the second half of WORK
251: *
252: DO 130 I = J, N
253: *
254: IF( J.GT.1 ) THEN
255: WORK( I ) = WORK( I ) +
256: $ DBLE( DCONJG( A( J-1, I ) )*
257: $ A( J-1, I ) )
258: END IF
259: WORK( N+I ) = DBLE( A( I, I ) ) - WORK( I )
260: *
261: 130 CONTINUE
262: *
263: IF( J.GT.1 ) THEN
264: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
265: PVT = ITEMP + J - 1
266: AJJ = WORK( N+PVT )
267: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
268: A( J, J ) = AJJ
269: GO TO 190
270: END IF
271: END IF
272: *
273: IF( J.NE.PVT ) THEN
274: *
275: * Pivot OK, so can now swap pivot rows and columns
276: *
277: A( PVT, PVT ) = A( J, J )
278: CALL ZSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
279: IF( PVT.LT.N )
280: $ CALL ZSWAP( N-PVT, A( J, PVT+1 ), LDA,
281: $ A( PVT, PVT+1 ), LDA )
282: DO 140 I = J + 1, PVT - 1
283: ZTEMP = DCONJG( A( J, I ) )
284: A( J, I ) = DCONJG( A( I, PVT ) )
285: A( I, PVT ) = ZTEMP
286: 140 CONTINUE
287: A( J, PVT ) = DCONJG( A( J, PVT ) )
288: *
289: * Swap dot products and PIV
290: *
291: DTEMP = WORK( J )
292: WORK( J ) = WORK( PVT )
293: WORK( PVT ) = DTEMP
294: ITEMP = PIV( PVT )
295: PIV( PVT ) = PIV( J )
296: PIV( J ) = ITEMP
297: END IF
298: *
299: AJJ = SQRT( AJJ )
300: A( J, J ) = AJJ
301: *
302: * Compute elements J+1:N of row J
303: *
304: IF( J.LT.N ) THEN
305: CALL ZLACGV( J-1, A( 1, J ), 1 )
306: CALL ZGEMV( 'Trans', J-1, N-J, -CONE, A( 1, J+1 ), LDA,
307: $ A( 1, J ), 1, CONE, A( J, J+1 ), LDA )
308: CALL ZLACGV( J-1, A( 1, J ), 1 )
309: CALL ZDSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
310: END IF
311: *
312: 150 CONTINUE
313: *
314: ELSE
315: *
1.5 bertrand 316: * Compute the Cholesky factorization P**T * A * P = L * L**H
1.1 bertrand 317: *
318: DO 180 J = 1, N
319: *
320: * Find pivot, test for exit, else swap rows and columns
321: * Update dot products, compute possible pivots which are
322: * stored in the second half of WORK
323: *
324: DO 160 I = J, N
325: *
326: IF( J.GT.1 ) THEN
327: WORK( I ) = WORK( I ) +
328: $ DBLE( DCONJG( A( I, J-1 ) )*
329: $ A( I, J-1 ) )
330: END IF
331: WORK( N+I ) = DBLE( A( I, I ) ) - WORK( I )
332: *
333: 160 CONTINUE
334: *
335: IF( J.GT.1 ) THEN
336: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
337: PVT = ITEMP + J - 1
338: AJJ = WORK( N+PVT )
339: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
340: A( J, J ) = AJJ
341: GO TO 190
342: END IF
343: END IF
344: *
345: IF( J.NE.PVT ) THEN
346: *
347: * Pivot OK, so can now swap pivot rows and columns
348: *
349: A( PVT, PVT ) = A( J, J )
350: CALL ZSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
351: IF( PVT.LT.N )
352: $ CALL ZSWAP( N-PVT, A( PVT+1, J ), 1, A( PVT+1, PVT ),
353: $ 1 )
354: DO 170 I = J + 1, PVT - 1
355: ZTEMP = DCONJG( A( I, J ) )
356: A( I, J ) = DCONJG( A( PVT, I ) )
357: A( PVT, I ) = ZTEMP
358: 170 CONTINUE
359: A( PVT, J ) = DCONJG( A( PVT, J ) )
360: *
361: * Swap dot products and PIV
362: *
363: DTEMP = WORK( J )
364: WORK( J ) = WORK( PVT )
365: WORK( PVT ) = DTEMP
366: ITEMP = PIV( PVT )
367: PIV( PVT ) = PIV( J )
368: PIV( J ) = ITEMP
369: END IF
370: *
371: AJJ = SQRT( AJJ )
372: A( J, J ) = AJJ
373: *
374: * Compute elements J+1:N of column J
375: *
376: IF( J.LT.N ) THEN
377: CALL ZLACGV( J-1, A( J, 1 ), LDA )
378: CALL ZGEMV( 'No Trans', N-J, J-1, -CONE, A( J+1, 1 ),
379: $ LDA, A( J, 1 ), LDA, CONE, A( J+1, J ), 1 )
380: CALL ZLACGV( J-1, A( J, 1 ), LDA )
381: CALL ZDSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
382: END IF
383: *
384: 180 CONTINUE
385: *
386: END IF
387: *
388: * Ran to completion, A has full rank
389: *
390: RANK = N
391: *
392: GO TO 200
393: 190 CONTINUE
394: *
395: * Rank is number of steps completed. Set INFO = 1 to signal
396: * that the factorization cannot be used to solve a system.
397: *
398: RANK = J - 1
399: INFO = 1
400: *
401: 200 CONTINUE
402: RETURN
403: *
404: * End of ZPSTF2
405: *
406: 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