001:       SUBROUTINE ZPSTRF( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
002: *
003: *  -- LAPACK routine (version 3.2.1)                                  --
004: *     
005: *  -- Contributed by Craig Lucas, University of Manchester / NAG Ltd. --
006: *  -- April 2009                                                      --
007: *
008: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
009: *
010: *     .. Scalar Arguments ..
011:       DOUBLE PRECISION   TOL
012:       INTEGER            INFO, LDA, N, RANK
013:       CHARACTER          UPLO
014: *     ..
015: *     .. Array Arguments ..
016:       COMPLEX*16         A( LDA, * )
017:       DOUBLE PRECISION   WORK( 2*N )
018:       INTEGER            PIV( N )
019: *     ..
020: *
021: *  Purpose
022: *  =======
023: *
024: *  ZPSTRF computes the Cholesky factorization with complete
025: *  pivoting of a complex Hermitian positive semidefinite matrix A.
026: *
027: *  The factorization has the form
028: *     P' * A * P = U' * U ,  if UPLO = 'U',
029: *     P' * A * P = L  * L',  if UPLO = 'L',
030: *  where U is an upper triangular matrix and L is lower triangular, and
031: *  P is stored as vector PIV.
032: *
033: *  This algorithm does not attempt to check that A is positive
034: *  semidefinite. This version of the algorithm calls level 3 BLAS.
035: *
036: *  Arguments
037: *  =========
038: *
039: *  UPLO    (input) CHARACTER*1
040: *          Specifies whether the upper or lower triangular part of the
041: *          symmetric matrix A is stored.
042: *          = 'U':  Upper triangular
043: *          = 'L':  Lower triangular
044: *
045: *  N       (input) INTEGER
046: *          The order of the matrix A.  N >= 0.
047: *
048: *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
049: *          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
050: *          n by n upper triangular part of A contains the upper
051: *          triangular part of the matrix A, and the strictly lower
052: *          triangular part of A is not referenced.  If UPLO = 'L', the
053: *          leading n by n lower triangular part of A contains the lower
054: *          triangular part of the matrix A, and the strictly upper
055: *          triangular part of A is not referenced.
056: *
057: *          On exit, if INFO = 0, the factor U or L from the Cholesky
058: *          factorization as above.
059: *
060: *  LDA     (input) INTEGER
061: *          The leading dimension of the array A.  LDA >= max(1,N).
062: *
063: *  PIV     (output) INTEGER array, dimension (N)
064: *          PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
065: *
066: *  RANK    (output) INTEGER
067: *          The rank of A given by the number of steps the algorithm
068: *          completed.
069: *
070: *  TOL     (input) DOUBLE PRECISION
071: *          User defined tolerance. If TOL < 0, then N*U*MAX( A(K,K) )
072: *          will be used. The algorithm terminates at the (K-1)st step
073: *          if the pivot <= TOL.
074: *
075: *  WORK    DOUBLE PRECISION array, dimension (2*N)
076: *          Work space.
077: *
078: *  INFO    (output) INTEGER
079: *          < 0: If INFO = -K, the K-th argument had an illegal value,
080: *          = 0: algorithm completed successfully, and
081: *          > 0: the matrix A is either rank deficient with computed rank
082: *               as returned in RANK, or is indefinite.  See Section 7 of
083: *               LAPACK Working Note #161 for further information.
084: *
085: *  =====================================================================
086: *
087: *     .. Parameters ..
088:       DOUBLE PRECISION   ONE, ZERO
089:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
090:       COMPLEX*16         CONE
091:       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
092: *     ..
093: *     .. Local Scalars ..
094:       COMPLEX*16         ZTEMP
095:       DOUBLE PRECISION   AJJ, DSTOP, DTEMP
096:       INTEGER            I, ITEMP, J, JB, K, NB, PVT
097:       LOGICAL            UPPER
098: *     ..
099: *     .. 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
389: