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