1: *> \brief \b ZPSTRF 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 ZPSTRF + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zpstrf.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zpstrf.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zpstrf.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZPSTRF( 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: *> ZPSTRF 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 3 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[in] LDA
87: *> \verbatim
88: *> LDA is INTEGER
89: *> The leading dimension of the array A. LDA >= max(1,N).
90: *> \endverbatim
91: *>
92: *> \param[out] PIV
93: *> \verbatim
94: *> PIV is INTEGER array, dimension (N)
95: *> PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
96: *> \endverbatim
97: *>
98: *> \param[out] RANK
99: *> \verbatim
100: *> RANK is INTEGER
101: *> The rank of A given by the number of steps the algorithm
102: *> completed.
103: *> \endverbatim
104: *>
105: *> \param[in] TOL
106: *> \verbatim
107: *> TOL is DOUBLE PRECISION
108: *> User defined tolerance. If TOL < 0, then N*U*MAX( A(K,K) )
109: *> will be used. The algorithm terminates at the (K-1)st step
110: *> if the pivot <= TOL.
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 ZPSTRF( 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, JB, K, NB, PVT
170: LOGICAL UPPER
171: * ..
172: * .. External Functions ..
173: DOUBLE PRECISION DLAMCH
174: INTEGER ILAENV
175: LOGICAL LSAME, DISNAN
176: EXTERNAL DLAMCH, ILAENV, LSAME, DISNAN
177: * ..
178: * .. External Subroutines ..
179: EXTERNAL ZDSCAL, ZGEMV, ZHERK, ZLACGV, ZPSTF2, ZSWAP,
180: $ XERBLA
181: * ..
182: * .. Intrinsic Functions ..
183: INTRINSIC DBLE, DCONJG, MAX, MIN, SQRT, MAXLOC
184: * ..
185: * .. Executable Statements ..
186: *
187: * Test the input parameters.
188: *
189: INFO = 0
190: UPPER = LSAME( UPLO, 'U' )
191: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
192: INFO = -1
193: ELSE IF( N.LT.0 ) THEN
194: INFO = -2
195: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
196: INFO = -4
197: END IF
198: IF( INFO.NE.0 ) THEN
199: CALL XERBLA( 'ZPSTRF', -INFO )
200: RETURN
201: END IF
202: *
203: * Quick return if possible
204: *
205: IF( N.EQ.0 )
206: $ RETURN
207: *
208: * Get block size
209: *
210: NB = ILAENV( 1, 'ZPOTRF', UPLO, N, -1, -1, -1 )
211: IF( NB.LE.1 .OR. NB.GE.N ) THEN
212: *
213: * Use unblocked code
214: *
215: CALL ZPSTF2( UPLO, N, A( 1, 1 ), LDA, PIV, RANK, TOL, WORK,
216: $ INFO )
217: GO TO 230
218: *
219: ELSE
220: *
221: * Initialize PIV
222: *
223: DO 100 I = 1, N
224: PIV( I ) = I
225: 100 CONTINUE
226: *
227: * Compute stopping value
228: *
229: DO 110 I = 1, N
230: WORK( I ) = DBLE( A( I, I ) )
231: 110 CONTINUE
232: PVT = MAXLOC( WORK( 1:N ), 1 )
233: AJJ = DBLE( A( PVT, PVT ) )
234: IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN
235: RANK = 0
236: INFO = 1
237: GO TO 230
238: END IF
239: *
240: * Compute stopping value if not supplied
241: *
242: IF( TOL.LT.ZERO ) THEN
243: DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ
244: ELSE
245: DSTOP = TOL
246: END IF
247: *
248: *
249: IF( UPPER ) THEN
250: *
251: * Compute the Cholesky factorization P**T * A * P = U**H * U
252: *
253: DO 160 K = 1, N, NB
254: *
255: * Account for last block not being NB wide
256: *
257: JB = MIN( NB, N-K+1 )
258: *
259: * Set relevant part of first half of WORK to zero,
260: * holds dot products
261: *
262: DO 120 I = K, N
263: WORK( I ) = 0
264: 120 CONTINUE
265: *
266: DO 150 J = K, K + JB - 1
267: *
268: * Find pivot, test for exit, else swap rows and columns
269: * Update dot products, compute possible pivots which are
270: * stored in the second half of WORK
271: *
272: DO 130 I = J, N
273: *
274: IF( J.GT.K ) THEN
275: WORK( I ) = WORK( I ) +
276: $ DBLE( DCONJG( A( J-1, I ) )*
277: $ A( J-1, I ) )
278: END IF
279: WORK( N+I ) = DBLE( A( I, I ) ) - WORK( I )
280: *
281: 130 CONTINUE
282: *
283: IF( J.GT.1 ) THEN
284: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
285: PVT = ITEMP + J - 1
286: AJJ = WORK( N+PVT )
287: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
288: A( J, J ) = AJJ
289: GO TO 220
290: END IF
291: END IF
292: *
293: IF( J.NE.PVT ) THEN
294: *
295: * Pivot OK, so can now swap pivot rows and columns
296: *
297: A( PVT, PVT ) = A( J, J )
298: CALL ZSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
299: IF( PVT.LT.N )
300: $ CALL ZSWAP( N-PVT, A( J, PVT+1 ), LDA,
301: $ A( PVT, PVT+1 ), LDA )
302: DO 140 I = J + 1, PVT - 1
303: ZTEMP = DCONJG( A( J, I ) )
304: A( J, I ) = DCONJG( A( I, PVT ) )
305: A( I, PVT ) = ZTEMP
306: 140 CONTINUE
307: A( J, PVT ) = DCONJG( A( J, PVT ) )
308: *
309: * Swap dot products and PIV
310: *
311: DTEMP = WORK( J )
312: WORK( J ) = WORK( PVT )
313: WORK( PVT ) = DTEMP
314: ITEMP = PIV( PVT )
315: PIV( PVT ) = PIV( J )
316: PIV( J ) = ITEMP
317: END IF
318: *
319: AJJ = SQRT( AJJ )
320: A( J, J ) = AJJ
321: *
322: * Compute elements J+1:N of row J.
323: *
324: IF( J.LT.N ) THEN
325: CALL ZLACGV( J-1, A( 1, J ), 1 )
326: CALL ZGEMV( 'Trans', J-K, N-J, -CONE, A( K, J+1 ),
327: $ LDA, A( K, J ), 1, CONE, A( J, J+1 ),
328: $ LDA )
329: CALL ZLACGV( J-1, A( 1, J ), 1 )
330: CALL ZDSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
331: END IF
332: *
333: 150 CONTINUE
334: *
335: * Update trailing matrix, J already incremented
336: *
337: IF( K+JB.LE.N ) THEN
338: CALL ZHERK( 'Upper', 'Conj Trans', N-J+1, JB, -ONE,
339: $ A( K, J ), LDA, ONE, A( J, J ), LDA )
340: END IF
341: *
342: 160 CONTINUE
343: *
344: ELSE
345: *
346: * Compute the Cholesky factorization P**T * A * P = L * L**H
347: *
348: DO 210 K = 1, N, NB
349: *
350: * Account for last block not being NB wide
351: *
352: JB = MIN( NB, N-K+1 )
353: *
354: * Set relevant part of first half of WORK to zero,
355: * holds dot products
356: *
357: DO 170 I = K, N
358: WORK( I ) = 0
359: 170 CONTINUE
360: *
361: DO 200 J = K, K + JB - 1
362: *
363: * Find pivot, test for exit, else swap rows and columns
364: * Update dot products, compute possible pivots which are
365: * stored in the second half of WORK
366: *
367: DO 180 I = J, N
368: *
369: IF( J.GT.K ) THEN
370: WORK( I ) = WORK( I ) +
371: $ DBLE( DCONJG( A( I, J-1 ) )*
372: $ A( I, J-1 ) )
373: END IF
374: WORK( N+I ) = DBLE( A( I, I ) ) - WORK( I )
375: *
376: 180 CONTINUE
377: *
378: IF( J.GT.1 ) THEN
379: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
380: PVT = ITEMP + J - 1
381: AJJ = WORK( N+PVT )
382: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
383: A( J, J ) = AJJ
384: GO TO 220
385: END IF
386: END IF
387: *
388: IF( J.NE.PVT ) THEN
389: *
390: * Pivot OK, so can now swap pivot rows and columns
391: *
392: A( PVT, PVT ) = A( J, J )
393: CALL ZSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
394: IF( PVT.LT.N )
395: $ CALL ZSWAP( N-PVT, A( PVT+1, J ), 1,
396: $ A( PVT+1, PVT ), 1 )
397: DO 190 I = J + 1, PVT - 1
398: ZTEMP = DCONJG( A( I, J ) )
399: A( I, J ) = DCONJG( A( PVT, I ) )
400: A( PVT, I ) = ZTEMP
401: 190 CONTINUE
402: A( PVT, J ) = DCONJG( A( PVT, J ) )
403: *
404: *
405: * Swap dot products and PIV
406: *
407: DTEMP = WORK( J )
408: WORK( J ) = WORK( PVT )
409: WORK( PVT ) = DTEMP
410: ITEMP = PIV( PVT )
411: PIV( PVT ) = PIV( J )
412: PIV( J ) = ITEMP
413: END IF
414: *
415: AJJ = SQRT( AJJ )
416: A( J, J ) = AJJ
417: *
418: * Compute elements J+1:N of column J.
419: *
420: IF( J.LT.N ) THEN
421: CALL ZLACGV( J-1, A( J, 1 ), LDA )
422: CALL ZGEMV( 'No Trans', N-J, J-K, -CONE,
423: $ A( J+1, K ), LDA, A( J, K ), LDA, CONE,
424: $ A( J+1, J ), 1 )
425: CALL ZLACGV( J-1, A( J, 1 ), LDA )
426: CALL ZDSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
427: END IF
428: *
429: 200 CONTINUE
430: *
431: * Update trailing matrix, J already incremented
432: *
433: IF( K+JB.LE.N ) THEN
434: CALL ZHERK( 'Lower', 'No Trans', N-J+1, JB, -ONE,
435: $ A( J, K ), LDA, ONE, A( J, J ), LDA )
436: END IF
437: *
438: 210 CONTINUE
439: *
440: END IF
441: END IF
442: *
443: * Ran to completion, A has full rank
444: *
445: RANK = N
446: *
447: GO TO 230
448: 220 CONTINUE
449: *
450: * Rank is the number of steps completed. Set INFO = 1 to signal
451: * that the factorization cannot be used to solve a system.
452: *
453: RANK = J - 1
454: INFO = 1
455: *
456: 230 CONTINUE
457: RETURN
458: *
459: * End of ZPSTRF
460: *
461: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zpstrf.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack