001:       SUBROUTINE CPSTRF( 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:       REAL               TOL
009:       INTEGER            INFO, LDA, N, RANK
010:       CHARACTER          UPLO
011: *     ..
012: *     .. Array Arguments ..
013:       COMPLEX            A( LDA, * )
014:       REAL               WORK( 2*N )
015:       INTEGER            PIV( N )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  CPSTRF 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 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) REAL
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    REAL 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:       REAL               ONE, ZERO
086:       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
087:       COMPLEX            CONE
088:       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
089: *     ..
090: *     .. Local Scalars ..
091:       COMPLEX            CTEMP
092:       REAL               AJJ, SSTOP, STEMP
093:       INTEGER            I, ITEMP, J, JB, K, NB, PVT
094:       LOGICAL            UPPER
095: *     ..
096: *     .. External Functions ..
097:       REAL               SLAMCH
098:       INTEGER            ILAENV
099:       LOGICAL            LSAME, SISNAN
100:       EXTERNAL           SLAMCH, ILAENV, LSAME, SISNAN
101: *     ..
102: *     .. External Subroutines ..
103:       EXTERNAL           CGEMV, CHERK, CLACGV, CPSTF2, CSSCAL, CSWAP,
104:      $                   XERBLA
105: *     ..
106: *     .. Intrinsic Functions ..
107:       INTRINSIC          CONJG, MAX, MIN, REAL, 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( 'CPSTRF', -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, 'CPOTRF', UPLO, N, -1, -1, -1 )
135:       IF( NB.LE.1 .OR. NB.GE.N ) THEN
136: *
137: *        Use unblocked code
138: *
139:          CALL CPSTF2( 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 ) = REAL( A( I, I ) )
155:   110    CONTINUE
156:          PVT = MAXLOC( WORK( 1:N ), 1 )
157:          AJJ = REAL( A( PVT, PVT ) )
158:          IF( AJJ.EQ.ZERO.OR.SISNAN( 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:             SSTOP = N * SLAMCH( 'Epsilon' ) * AJJ
168:          ELSE
169:             SSTOP = 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:      $                              REAL( CONJG( A( J-1, I ) )*
201:      $                                    A( J-1, I ) )
202:                      END IF
203:                      WORK( N+I ) = REAL( 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.SSTOP.OR.SISNAN( 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 CSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
223:                      IF( PVT.LT.N )
224:      $                  CALL CSWAP( N-PVT, A( J, PVT+1 ), LDA,
225:      $                              A( PVT, PVT+1 ), LDA )
226:                      DO 140 I = J + 1, PVT - 1
227:                         CTEMP = CONJG( A( J, I ) )
228:                         A( J, I ) = CONJG( A( I, PVT ) )
229:                         A( I, PVT ) = CTEMP
230:   140                CONTINUE
231:                      A( J, PVT ) = CONJG( A( J, PVT ) )
232: *
233: *                    Swap dot products and PIV
234: *
235:                      STEMP = WORK( J )
236:                      WORK( J ) = WORK( PVT )
237:                      WORK( PVT ) = STEMP
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 CLACGV( J-1, A( 1, J ), 1 )
250:                      CALL CGEMV( '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 CLACGV( J-1, A( 1, J ), 1 )
254:                      CALL CSSCAL( 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 CHERK( '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:      $                              REAL( CONJG( A( I, J-1 ) )*
296:      $                                    A( I, J-1 ) )
297:                      END IF
298:                      WORK( N+I ) = REAL( 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.SSTOP.OR.SISNAN( 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 CSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
318:                      IF( PVT.LT.N )
319:      $                  CALL CSWAP( N-PVT, A( PVT+1, J ), 1,
320:      $                              A( PVT+1, PVT ), 1 )
321:                      DO 190 I = J + 1, PVT - 1
322:                         CTEMP = CONJG( A( I, J ) )
323:                         A( I, J ) = CONJG( A( PVT, I ) )
324:                         A( PVT, I ) = CTEMP
325:   190                CONTINUE
326:                      A( PVT, J ) = CONJG( A( PVT, J ) )
327: *
328: *                    Swap dot products and PIV
329: *
330:                      STEMP = WORK( J )
331:                      WORK( J ) = WORK( PVT )
332:                      WORK( PVT ) = STEMP
333:                      ITEMP = PIV( PVT )
334:                      PIV( PVT ) = PIV( J )
335:                      PIV( J ) = ITEMP
336:                   END IF
337: *
338:                   AJJ = SQRT( AJJ )
339:                   A( J, J ) = AJJ
340: *
341: *                 Compute elements J+1:N of column J.
342: *
343:                   IF( J.LT.N ) THEN
344:                      CALL CLACGV( J-1, A( J, 1 ), LDA )
345:                      CALL CGEMV( 'No Trans', N-J, J-K, -CONE,
346:      $                           A( J+1, K ), LDA, A( J, K ), LDA, CONE,
347:      $                           A( J+1, J ), 1 )
348:                      CALL CLACGV( J-1, A( J, 1 ), LDA )
349:                      CALL CSSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
350:                   END IF
351: *
352:   200          CONTINUE
353: *
354: *              Update trailing matrix, J already incremented
355: *
356:                IF( K+JB.LE.N ) THEN
357:                   CALL CHERK( 'Lower', 'No Trans', N-J+1, JB, -ONE,
358:      $                        A( J, K ), LDA, ONE, A( J, J ), LDA )
359:                END IF
360: *
361:   210       CONTINUE
362: *
363:          END IF
364:       END IF
365: *
366: *     Ran to completion, A has full rank
367: *
368:       RANK = N
369: *
370:       GO TO 230
371:   220 CONTINUE
372: *
373: *     Rank is the number of steps completed.  Set INFO = 1 to signal
374: *     that the factorization cannot be used to solve a system.
375: *
376:       RANK = J - 1
377:       INFO = 1
378: *
379:   230 CONTINUE
380:       RETURN
381: *
382: *     End of CPSTRF
383: *
384:       END
385: