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