001:       SUBROUTINE ZPSTRF( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
002: *
003: *  -- LAPACK routine (version 3.2) --
004: *     Craig Lucas, University of Manchester / NAG Ltd.
005: *     October, 2008
006: *
007: *     .. Scalar Arguments ..
008:       DOUBLE PRECISION   TOL
009:       INTEGER            INFO, LDA, N, RANK
010:       CHARACTER          UPLO
011: *     ..
012: *     .. Array Arguments ..
013:       COMPLEX*16         A( LDA, * )
014:       DOUBLE PRECISION   WORK( 2*N )
015:       INTEGER            PIV( N )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  ZPSTRF computes the Cholesky factorization with complete
022: *  pivoting of a complex Hermitian positive semidefinite matrix A.
023: *
024: *  The factorization has the form
025: *     P' * A * P = U' * U ,  if UPLO = 'U',
026: *     P' * A * P = L  * L',  if UPLO = 'L',
027: *  where U is an upper triangular matrix and L is lower triangular, and
028: *  P is stored as vector PIV.
029: *
030: *  This algorithm does not attempt to check that A is positive
031: *  semidefinite. This version of the algorithm calls level 3 BLAS.
032: *
033: *  Arguments
034: *  =========
035: *
036: *  UPLO    (input) CHARACTER*1
037: *          Specifies whether the upper or lower triangular part of the
038: *          symmetric matrix A is stored.
039: *          = 'U':  Upper triangular
040: *          = 'L':  Lower triangular
041: *
042: *  N       (input) INTEGER
043: *          The order of the matrix A.  N >= 0.
044: *
045: *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
046: *          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
047: *          n by n upper triangular part of A contains the upper
048: *          triangular part of the matrix A, and the strictly lower
049: *          triangular part of A is not referenced.  If UPLO = 'L', the
050: *          leading n by n lower triangular part of A contains the lower
051: *          triangular part of the matrix A, and the strictly upper
052: *          triangular part of A is not referenced.
053: *
054: *          On exit, if INFO = 0, the factor U or L from the Cholesky
055: *          factorization as above.
056: *
057: *  PIV     (output) INTEGER array, dimension (N)
058: *          PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
059: *
060: *  RANK    (output) INTEGER
061: *          The rank of A given by the number of steps the algorithm
062: *          completed.
063: *
064: *  TOL     (input) DOUBLE PRECISION
065: *          User defined tolerance. If TOL < 0, then N*U*MAX( A(K,K) )
066: *          will be used. The algorithm terminates at the (K-1)st step
067: *          if the pivot <= TOL.
068: *
069: *  LDA     (input) INTEGER
070: *          The leading dimension of the array A.  LDA >= max(1,N).
071: *
072: *  WORK    DOUBLE PRECISION array, dimension (2*N)
073: *          Work space.
074: *
075: *  INFO    (output) INTEGER
076: *          < 0: If INFO = -K, the K-th argument had an illegal value,
077: *          = 0: algorithm completed successfully, and
078: *          > 0: the matrix A is either rank deficient with computed rank
079: *               as returned in RANK, or is indefinite.  See Section 7 of
080: *               LAPACK Working Note #161 for further information.
081: *
082: *  =====================================================================
083: *
084: *     .. Parameters ..
085:       DOUBLE PRECISION   ONE, ZERO
086:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
087:       COMPLEX*16         CONE
088:       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
089: *     ..
090: *     .. Local Scalars ..
091:       COMPLEX*16         ZTEMP
092:       DOUBLE PRECISION   AJJ, DSTOP, DTEMP
093:       INTEGER            I, ITEMP, J, JB, K, NB, PVT
094:       LOGICAL            UPPER
095: *     ..
096: *     .. External Functions ..
097:       DOUBLE PRECISION   DLAMCH
098:       INTEGER            ILAENV
099:       LOGICAL            LSAME, DISNAN
100:       EXTERNAL           DLAMCH, ILAENV, LSAME, DISNAN
101: *     ..
102: *     .. External Subroutines ..
103:       EXTERNAL           ZDSCAL, ZGEMV, ZHERK, ZLACGV, ZPSTF2, ZSWAP
104:      $                   XERBLA
105: *     ..
106: *     .. Intrinsic Functions ..
107:       INTRINSIC          DBLE, DCONJG, MAX, MIN, SQRT, MAXLOC
108: *     ..
109: *     .. Executable Statements ..
110: *
111: *     Test the input parameters.
112: *
113:       INFO = 0
114:       UPPER = LSAME( UPLO, 'U' )
115:       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
116:          INFO = -1
117:       ELSE IF( N.LT.0 ) THEN
118:          INFO = -2
119:       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
120:          INFO = -4
121:       END IF
122:       IF( INFO.NE.0 ) THEN
123:          CALL XERBLA( 'ZPSTRF', -INFO )
124:          RETURN
125:       END IF
126: *
127: *     Quick return if possible
128: *
129:       IF( N.EQ.0 )
130:      $   RETURN
131: *
132: *     Get block size
133: *
134:       NB = ILAENV( 1, 'ZPOTRF', UPLO, N, -1, -1, -1 )
135:       IF( NB.LE.1 .OR. NB.GE.N ) THEN
136: *
137: *        Use unblocked code
138: *
139:          CALL ZPSTF2( UPLO, N, A( 1, 1 ), LDA, PIV, RANK, TOL, WORK,
140:      $                INFO )
141:          GO TO 230
142: *
143:       ELSE
144: *
145: *     Initialize PIV
146: *
147:          DO 100 I = 1, N
148:             PIV( I ) = I
149:   100    CONTINUE
150: *
151: *     Compute stopping value
152: *
153:          DO 110 I = 1, N
154:             WORK( I ) = DBLE( A( I, I ) )
155:   110    CONTINUE
156:          PVT = MAXLOC( WORK( 1:N ), 1 )
157:          AJJ = DBLE( A( PVT, PVT ) )
158:          IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
159:             RANK = 0
160:             INFO = 1
161:             GO TO 230
162:          END IF
163: *
164: *     Compute stopping value if not supplied
165: *
166:          IF( TOL.LT.ZERO ) THEN
167:             DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ
168:          ELSE
169:             DSTOP = TOL
170:          END IF
171: *
172: *
173:          IF( UPPER ) THEN
174: *
175: *           Compute the Cholesky factorization P' * A * P = U' * U
176: *
177:             DO 160 K = 1, N, NB
178: *
179: *              Account for last block not being NB wide
180: *
181:                JB = MIN( NB, N-K+1 )
182: *
183: *              Set relevant part of first half of WORK to zero,
184: *              holds dot products
185: *
186:                DO 120 I = K, N
187:                   WORK( I ) = 0
188:   120          CONTINUE
189: *
190:                DO 150 J = K, K + JB - 1
191: *
192: *              Find pivot, test for exit, else swap rows and columns
193: *              Update dot products, compute possible pivots which are
194: *              stored in the second half of WORK
195: *
196:                   DO 130 I = J, N
197: *
198:                      IF( J.GT.K ) THEN
199:                         WORK( I ) = WORK( I ) +
200:      $                              DBLE( DCONJG( A( J-1, I ) )*
201:      $                                    A( J-1, I ) )
202:                      END IF
203:                      WORK( N+I ) = DBLE( A( I, I ) ) - WORK( I )
204: *
205:   130             CONTINUE
206: *
207:                   IF( J.GT.1 ) THEN
208:                      ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
209:                      PVT = ITEMP + J - 1
210:                      AJJ = WORK( N+PVT )
211:                      IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
212:                         A( J, J ) = AJJ
213:                         GO TO 220
214:                      END IF
215:                   END IF
216: *
217:                   IF( J.NE.PVT ) THEN
218: *
219: *                    Pivot OK, so can now swap pivot rows and columns
220: *
221:                      A( PVT, PVT ) = A( J, J )
222:                      CALL ZSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
223:                      IF( PVT.LT.N )
224:      $                  CALL ZSWAP( N-PVT, A( J, PVT+1 ), LDA,
225:      $                              A( PVT, PVT+1 ), LDA )
226:                      DO 140 I = J + 1, PVT - 1
227:                         ZTEMP = DCONJG( A( J, I ) )
228:                         A( J, I ) = DCONJG( A( I, PVT ) )
229:                         A( I, PVT ) = ZTEMP
230:   140                CONTINUE
231:                      A( J, PVT ) = DCONJG( A( J, PVT ) )
232: *
233: *                    Swap dot products and PIV
234: *
235:                      DTEMP = WORK( J )
236:                      WORK( J ) = WORK( PVT )
237:                      WORK( PVT ) = DTEMP
238:                      ITEMP = PIV( PVT )
239:                      PIV( PVT ) = PIV( J )
240:                      PIV( J ) = ITEMP
241:                   END IF
242: *
243:                   AJJ = SQRT( AJJ )
244:                   A( J, J ) = AJJ
245: *
246: *                 Compute elements J+1:N of row J.
247: *
248:                   IF( J.LT.N ) THEN
249:                      CALL ZLACGV( J-1, A( 1, J ), 1 )
250:                      CALL ZGEMV( 'Trans', J-K, N-J, -CONE, A( K, J+1 ),
251:      $                           LDA, A( K, J ), 1, CONE, A( J, J+1 ),
252:      $                           LDA )
253:                      CALL ZLACGV( J-1, A( 1, J ), 1 )
254:                      CALL ZDSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
255:                   END IF
256: *
257:   150          CONTINUE
258: *
259: *              Update trailing matrix, J already incremented
260: *
261:                IF( K+JB.LE.N ) THEN
262:                   CALL ZHERK( 'Upper', 'Conj Trans', N-J+1, JB, -ONE,
263:      $                        A( K, J ), LDA, ONE, A( J, J ), LDA )
264:                END IF
265: *
266:   160       CONTINUE
267: *
268:          ELSE
269: *
270: *        Compute the Cholesky factorization P' * A * P = L * L'
271: *
272:             DO 210 K = 1, N, NB
273: *
274: *              Account for last block not being NB wide
275: *
276:                JB = MIN( NB, N-K+1 )
277: *
278: *              Set relevant part of first half of WORK to zero,
279: *              holds dot products
280: *
281:                DO 170 I = K, N
282:                   WORK( I ) = 0
283:   170          CONTINUE
284: *
285:                DO 200 J = K, K + JB - 1
286: *
287: *              Find pivot, test for exit, else swap rows and columns
288: *              Update dot products, compute possible pivots which are
289: *              stored in the second half of WORK
290: *
291:                   DO 180 I = J, N
292: *
293:                      IF( J.GT.K ) THEN
294:                         WORK( I ) = WORK( I ) +
295:      $                              DBLE( DCONJG( A( I, J-1 ) )*
296:      $                                    A( I, J-1 ) )
297:                      END IF
298:                      WORK( N+I ) = DBLE( A( I, I ) ) - WORK( I )
299: *
300:   180             CONTINUE
301: *
302:                   IF( J.GT.1 ) THEN
303:                      ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
304:                      PVT = ITEMP + J - 1
305:                      AJJ = WORK( N+PVT )
306:                      IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
307:                         A( J, J ) = AJJ
308:                         GO TO 220
309:                      END IF
310:                   END IF
311: *
312:                   IF( J.NE.PVT ) THEN
313: *
314: *                    Pivot OK, so can now swap pivot rows and columns
315: *
316:                      A( PVT, PVT ) = A( J, J )
317:                      CALL ZSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
318:                      IF( PVT.LT.N )
319:      $                  CALL ZSWAP( N-PVT, A( PVT+1, J ), 1,
320:      $                              A( PVT+1, PVT ), 1 )
321:                      DO 190 I = J + 1, PVT - 1
322:                         ZTEMP = DCONJG( A( I, J ) )
323:                         A( I, J ) = DCONJG( A( PVT, I ) )
324:                         A( PVT, I ) = ZTEMP
325:   190                CONTINUE
326:                      A( PVT, J ) = DCONJG( A( PVT, J ) )
327: *
328: *
329: *                    Swap dot products and PIV
330: *
331:                      DTEMP = WORK( J )
332:                      WORK( J ) = WORK( PVT )
333:                      WORK( PVT ) = DTEMP
334:                      ITEMP = PIV( PVT )
335:                      PIV( PVT ) = PIV( J )
336:                      PIV( J ) = ITEMP
337:                   END IF
338: *
339:                   AJJ = SQRT( AJJ )
340:                   A( J, J ) = AJJ
341: *
342: *                 Compute elements J+1:N of column J.
343: *
344:                   IF( J.LT.N ) THEN
345:                      CALL ZLACGV( J-1, A( J, 1 ), LDA )
346:                      CALL ZGEMV( 'No Trans', N-J, J-K, -CONE,
347:      $                           A( J+1, K ), LDA, A( J, K ), LDA, CONE,
348:      $                           A( J+1, J ), 1 )
349:                      CALL ZLACGV( J-1, A( J, 1 ), LDA )
350:                      CALL ZDSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
351:                   END IF
352: *
353:   200          CONTINUE
354: *
355: *              Update trailing matrix, J already incremented
356: *
357:                IF( K+JB.LE.N ) THEN
358:                   CALL ZHERK( 'Lower', 'No Trans', N-J+1, JB, -ONE,
359:      $                        A( J, K ), LDA, ONE, A( J, J ), LDA )
360:                END IF
361: *
362:   210       CONTINUE
363: *
364:          END IF
365:       END IF
366: *
367: *     Ran to completion, A has full rank
368: *
369:       RANK = N
370: *
371:       GO TO 230
372:   220 CONTINUE
373: *
374: *     Rank is the number of steps completed.  Set INFO = 1 to signal
375: *     that the factorization cannot be used to solve a system.
376: *
377:       RANK = J - 1
378:       INFO = 1
379: *
380:   230 CONTINUE
381:       RETURN
382: *
383: *     End of ZPSTRF
384: *
385:       END
386: