![[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|
Return to zpstrf.f CVS log | Up to [local] / rpl / lapack / lapack |
1.11 bertrand 1: *> \brief \b ZPSTRF computes the Cholesky factorization with complete pivoting of a complex Hermitian positive semidefinite matrix.
1.6 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.13 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.6 bertrand 7: *
8: *> \htmlonly
1.13 bertrand 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">
1.6 bertrand 15: *> [TXT]</a>
1.13 bertrand 16: *> \endhtmlonly
1.6 bertrand 17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZPSTRF( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
1.13 bertrand 22: *
1.6 bertrand 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: * ..
1.13 bertrand 33: *
1.6 bertrand 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
1.11 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: *
1.13 bertrand 133: *> \author Univ. of Tennessee
134: *> \author Univ. of California Berkeley
135: *> \author Univ. of Colorado Denver
136: *> \author NAG Ltd.
1.6 bertrand 137: *
1.13 bertrand 138: *> \date December 2016
1.6 bertrand 139: *
140: *> \ingroup complex16OTHERcomputational
141: *
142: * =====================================================================
1.1 bertrand 143: SUBROUTINE ZPSTRF( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
144: *
1.13 bertrand 145: * -- LAPACK computational routine (version 3.7.0) --
1.1 bertrand 146: * -- LAPACK is a software package provided by Univ. of Tennessee, --
1.6 bertrand 147: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.13 bertrand 148: * December 2016
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, JB, K, NB, PVT
173: LOGICAL UPPER
174: * ..
175: * .. External Functions ..
176: DOUBLE PRECISION DLAMCH
177: INTEGER ILAENV
178: LOGICAL LSAME, DISNAN
179: EXTERNAL DLAMCH, ILAENV, LSAME, DISNAN
180: * ..
181: * .. External Subroutines ..
182: EXTERNAL ZDSCAL, ZGEMV, ZHERK, ZLACGV, ZPSTF2, ZSWAP,
183: $ XERBLA
184: * ..
185: * .. Intrinsic Functions ..
186: INTRINSIC DBLE, DCONJG, MAX, MIN, SQRT, MAXLOC
187: * ..
188: * .. Executable Statements ..
189: *
190: * Test the input parameters.
191: *
192: INFO = 0
193: UPPER = LSAME( UPLO, 'U' )
194: IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
195: INFO = -1
196: ELSE IF( N.LT.0 ) THEN
197: INFO = -2
198: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
199: INFO = -4
200: END IF
201: IF( INFO.NE.0 ) THEN
202: CALL XERBLA( 'ZPSTRF', -INFO )
203: RETURN
204: END IF
205: *
206: * Quick return if possible
207: *
208: IF( N.EQ.0 )
209: $ RETURN
210: *
211: * Get block size
212: *
213: NB = ILAENV( 1, 'ZPOTRF', UPLO, N, -1, -1, -1 )
214: IF( NB.LE.1 .OR. NB.GE.N ) THEN
215: *
216: * Use unblocked code
217: *
218: CALL ZPSTF2( UPLO, N, A( 1, 1 ), LDA, PIV, RANK, TOL, WORK,
219: $ INFO )
220: GO TO 230
221: *
222: ELSE
223: *
224: * Initialize PIV
225: *
226: DO 100 I = 1, N
227: PIV( I ) = I
228: 100 CONTINUE
229: *
230: * Compute stopping value
231: *
232: DO 110 I = 1, N
233: WORK( I ) = DBLE( A( I, I ) )
234: 110 CONTINUE
235: PVT = MAXLOC( WORK( 1:N ), 1 )
236: AJJ = DBLE( A( PVT, PVT ) )
1.11 bertrand 237: IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN
1.1 bertrand 238: RANK = 0
239: INFO = 1
240: GO TO 230
241: END IF
242: *
243: * Compute stopping value if not supplied
244: *
245: IF( TOL.LT.ZERO ) THEN
246: DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ
247: ELSE
248: DSTOP = TOL
249: END IF
250: *
251: *
252: IF( UPPER ) THEN
253: *
1.5 bertrand 254: * Compute the Cholesky factorization P**T * A * P = U**H * U
1.1 bertrand 255: *
256: DO 160 K = 1, N, NB
257: *
258: * Account for last block not being NB wide
259: *
260: JB = MIN( NB, N-K+1 )
261: *
262: * Set relevant part of first half of WORK to zero,
263: * holds dot products
264: *
265: DO 120 I = K, N
266: WORK( I ) = 0
267: 120 CONTINUE
268: *
269: DO 150 J = K, K + JB - 1
270: *
271: * Find pivot, test for exit, else swap rows and columns
272: * Update dot products, compute possible pivots which are
273: * stored in the second half of WORK
274: *
275: DO 130 I = J, N
276: *
277: IF( J.GT.K ) THEN
278: WORK( I ) = WORK( I ) +
279: $ DBLE( DCONJG( A( J-1, I ) )*
280: $ A( J-1, I ) )
281: END IF
282: WORK( N+I ) = DBLE( A( I, I ) ) - WORK( I )
283: *
284: 130 CONTINUE
285: *
286: IF( J.GT.1 ) THEN
287: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
288: PVT = ITEMP + J - 1
289: AJJ = WORK( N+PVT )
290: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
291: A( J, J ) = AJJ
292: GO TO 220
293: END IF
294: END IF
295: *
296: IF( J.NE.PVT ) THEN
297: *
298: * Pivot OK, so can now swap pivot rows and columns
299: *
300: A( PVT, PVT ) = A( J, J )
301: CALL ZSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
302: IF( PVT.LT.N )
303: $ CALL ZSWAP( N-PVT, A( J, PVT+1 ), LDA,
304: $ A( PVT, PVT+1 ), LDA )
305: DO 140 I = J + 1, PVT - 1
306: ZTEMP = DCONJG( A( J, I ) )
307: A( J, I ) = DCONJG( A( I, PVT ) )
308: A( I, PVT ) = ZTEMP
309: 140 CONTINUE
310: A( J, PVT ) = DCONJG( A( J, PVT ) )
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 row J.
326: *
327: IF( J.LT.N ) THEN
328: CALL ZLACGV( J-1, A( 1, J ), 1 )
329: CALL ZGEMV( 'Trans', J-K, N-J, -CONE, A( K, J+1 ),
330: $ LDA, A( K, J ), 1, CONE, A( J, J+1 ),
331: $ LDA )
332: CALL ZLACGV( J-1, A( 1, J ), 1 )
333: CALL ZDSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
334: END IF
335: *
336: 150 CONTINUE
337: *
338: * Update trailing matrix, J already incremented
339: *
340: IF( K+JB.LE.N ) THEN
341: CALL ZHERK( 'Upper', 'Conj Trans', N-J+1, JB, -ONE,
342: $ A( K, J ), LDA, ONE, A( J, J ), LDA )
343: END IF
344: *
345: 160 CONTINUE
346: *
347: ELSE
348: *
1.5 bertrand 349: * Compute the Cholesky factorization P**T * A * P = L * L**H
1.1 bertrand 350: *
351: DO 210 K = 1, N, NB
352: *
353: * Account for last block not being NB wide
354: *
355: JB = MIN( NB, N-K+1 )
356: *
357: * Set relevant part of first half of WORK to zero,
358: * holds dot products
359: *
360: DO 170 I = K, N
361: WORK( I ) = 0
362: 170 CONTINUE
363: *
364: DO 200 J = K, K + JB - 1
365: *
366: * Find pivot, test for exit, else swap rows and columns
367: * Update dot products, compute possible pivots which are
368: * stored in the second half of WORK
369: *
370: DO 180 I = J, N
371: *
372: IF( J.GT.K ) THEN
373: WORK( I ) = WORK( I ) +
374: $ DBLE( DCONJG( A( I, J-1 ) )*
375: $ A( I, J-1 ) )
376: END IF
377: WORK( N+I ) = DBLE( A( I, I ) ) - WORK( I )
378: *
379: 180 CONTINUE
380: *
381: IF( J.GT.1 ) THEN
382: ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
383: PVT = ITEMP + J - 1
384: AJJ = WORK( N+PVT )
385: IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
386: A( J, J ) = AJJ
387: GO TO 220
388: END IF
389: END IF
390: *
391: IF( J.NE.PVT ) THEN
392: *
393: * Pivot OK, so can now swap pivot rows and columns
394: *
395: A( PVT, PVT ) = A( J, J )
396: CALL ZSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
397: IF( PVT.LT.N )
398: $ CALL ZSWAP( N-PVT, A( PVT+1, J ), 1,
399: $ A( PVT+1, PVT ), 1 )
400: DO 190 I = J + 1, PVT - 1
401: ZTEMP = DCONJG( A( I, J ) )
402: A( I, J ) = DCONJG( A( PVT, I ) )
403: A( PVT, I ) = ZTEMP
404: 190 CONTINUE
405: A( PVT, J ) = DCONJG( A( PVT, J ) )
406: *
407: *
408: * Swap dot products and PIV
409: *
410: DTEMP = WORK( J )
411: WORK( J ) = WORK( PVT )
412: WORK( PVT ) = DTEMP
413: ITEMP = PIV( PVT )
414: PIV( PVT ) = PIV( J )
415: PIV( J ) = ITEMP
416: END IF
417: *
418: AJJ = SQRT( AJJ )
419: A( J, J ) = AJJ
420: *
421: * Compute elements J+1:N of column J.
422: *
423: IF( J.LT.N ) THEN
424: CALL ZLACGV( J-1, A( J, 1 ), LDA )
425: CALL ZGEMV( 'No Trans', N-J, J-K, -CONE,
426: $ A( J+1, K ), LDA, A( J, K ), LDA, CONE,
427: $ A( J+1, J ), 1 )
428: CALL ZLACGV( J-1, A( J, 1 ), LDA )
429: CALL ZDSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
430: END IF
431: *
432: 200 CONTINUE
433: *
434: * Update trailing matrix, J already incremented
435: *
436: IF( K+JB.LE.N ) THEN
437: CALL ZHERK( 'Lower', 'No Trans', N-J+1, JB, -ONE,
438: $ A( J, K ), LDA, ONE, A( J, J ), LDA )
439: END IF
440: *
441: 210 CONTINUE
442: *
443: END IF
444: END IF
445: *
446: * Ran to completion, A has full rank
447: *
448: RANK = N
449: *
450: GO TO 230
451: 220 CONTINUE
452: *
453: * Rank is the number of steps completed. Set INFO = 1 to signal
454: * that the factorization cannot be used to solve a system.
455: *
456: RANK = J - 1
457: INFO = 1
458: *
459: 230 CONTINUE
460: RETURN
461: *
462: * End of ZPSTRF
463: *
464: END