1 *> \brief \b SPSTF2 computes the Cholesky factorization with complete pivoting of a real symmetric positive semidefinite matrix.
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download SPSTF2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/spstf2.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/spstf2.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/spstf2.f">
21 * SUBROUTINE SPSTF2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
23 * .. Scalar Arguments ..
25 * INTEGER INFO, LDA, N, RANK
28 * .. Array Arguments ..
29 * REAL A( LDA, * ), WORK( 2*N )
39 *> SPSTF2 computes the Cholesky factorization with complete
40 *> pivoting of a real symmetric positive semidefinite matrix A.
42 *> The factorization has the form
43 *> P**T * A * P = U**T * U , if UPLO = 'U',
44 *> P**T * A * P = L * L**T, if UPLO = 'L',
45 *> where U is an upper triangular matrix and L is lower triangular, and
46 *> P is stored as vector PIV.
48 *> This algorithm does not attempt to check that A is positive
49 *> semidefinite. This version of the algorithm calls level 2 BLAS.
57 *> UPLO is CHARACTER*1
58 *> Specifies whether the upper or lower triangular part of the
59 *> symmetric matrix A is stored.
60 *> = 'U': Upper triangular
61 *> = 'L': Lower triangular
67 *> The order of the matrix A. N >= 0.
72 *> A is REAL array, dimension (LDA,N)
73 *> On entry, the symmetric matrix A. If UPLO = 'U', the leading
74 *> n by n upper triangular part of A contains the upper
75 *> triangular part of the matrix A, and the strictly lower
76 *> triangular part of A is not referenced. If UPLO = 'L', the
77 *> leading n by n lower triangular part of A contains the lower
78 *> triangular part of the matrix A, and the strictly upper
79 *> triangular part of A is not referenced.
81 *> On exit, if INFO = 0, the factor U or L from the Cholesky
82 *> factorization as above.
87 *> PIV is INTEGER array, dimension (N)
88 *> PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
94 *> The rank of A given by the number of steps the algorithm
101 *> User defined tolerance. If TOL < 0, then N*U*MAX( A( K,K ) )
102 *> will be used. The algorithm terminates at the (K-1)st step
103 *> if the pivot <= TOL.
109 *> The leading dimension of the array A. LDA >= max(1,N).
114 *> WORK is REAL array, dimension (2*N)
121 *> < 0: If INFO = -K, the K-th argument had an illegal value,
122 *> = 0: algorithm completed successfully, and
123 *> > 0: the matrix A is either rank deficient with computed rank
124 *> as returned in RANK, or is not positive semidefinite. See
125 *> Section 7 of LAPACK Working Note #161 for further
132 *> \author Univ. of Tennessee
133 *> \author Univ. of California Berkeley
134 *> \author Univ. of Colorado Denver
137 *> \date November 2015
139 *> \ingroup realOTHERcomputational
141 * =====================================================================
142 SUBROUTINE SPSTF2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
144 * -- LAPACK computational routine (version 3.6.0) --
145 * -- LAPACK is a software package provided by Univ. of Tennessee, --
146 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
149 * .. Scalar Arguments ..
151 INTEGER INFO, LDA, N, RANK
154 * .. Array Arguments ..
155 REAL A( LDA, * ), WORK( 2*N )
159 * =====================================================================
163 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
165 * .. Local Scalars ..
166 REAL AJJ, SSTOP, STEMP
167 INTEGER I, ITEMP, J, PVT
170 * .. External Functions ..
172 LOGICAL LSAME, SISNAN
173 EXTERNAL SLAMCH, LSAME, SISNAN
175 * .. External Subroutines ..
176 EXTERNAL SGEMV, SSCAL, SSWAP, XERBLA
178 * .. Intrinsic Functions ..
179 INTRINSIC MAX, SQRT, MAXLOC
181 * .. Executable Statements ..
183 * Test the input parameters
186 UPPER = LSAME( UPLO, 'U' )
187 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
189 ELSE IF( N.LT.0 ) THEN
191 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
195 CALL XERBLA( 'SPSTF2', -INFO )
199 * Quick return if possible
210 * Compute stopping value
215 IF( A( I, I ).GT.AJJ ) THEN
220 IF( AJJ.LE.ZERO.OR.SISNAN( AJJ ) ) THEN
226 * Compute stopping value if not supplied
228 IF( TOL.LT.ZERO ) THEN
229 SSTOP = N * SLAMCH( 'Epsilon' ) * AJJ
234 * Set first half of WORK to zero, holds dot products
242 * Compute the Cholesky factorization P**T * A * P = U**T * U
246 * Find pivot, test for exit, else swap rows and columns
247 * Update dot products, compute possible pivots which are
248 * stored in the second half of WORK
253 WORK( I ) = WORK( I ) + A( J-1, I )**2
255 WORK( N+I ) = A( I, I ) - WORK( I )
260 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
263 IF( AJJ.LE.SSTOP.OR.SISNAN( AJJ ) ) THEN
271 * Pivot OK, so can now swap pivot rows and columns
273 A( PVT, PVT ) = A( J, J )
274 CALL SSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
276 $ CALL SSWAP( N-PVT, A( J, PVT+1 ), LDA,
277 $ A( PVT, PVT+1 ), LDA )
278 CALL SSWAP( PVT-J-1, A( J, J+1 ), LDA, A( J+1, PVT ), 1 )
280 * Swap dot products and PIV
283 WORK( J ) = WORK( PVT )
286 PIV( PVT ) = PIV( J )
293 * Compute elements J+1:N of row J
296 CALL SGEMV( 'Trans', J-1, N-J, -ONE, A( 1, J+1 ), LDA,
297 $ A( 1, J ), 1, ONE, A( J, J+1 ), LDA )
298 CALL SSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
305 * Compute the Cholesky factorization P**T * A * P = L * L**T
309 * Find pivot, test for exit, else swap rows and columns
310 * Update dot products, compute possible pivots which are
311 * stored in the second half of WORK
316 WORK( I ) = WORK( I ) + A( I, J-1 )**2
318 WORK( N+I ) = A( I, I ) - WORK( I )
323 ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
326 IF( AJJ.LE.SSTOP.OR.SISNAN( AJJ ) ) THEN
334 * Pivot OK, so can now swap pivot rows and columns
336 A( PVT, PVT ) = A( J, J )
337 CALL SSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
339 $ CALL SSWAP( N-PVT, A( PVT+1, J ), 1, A( PVT+1, PVT ),
341 CALL SSWAP( PVT-J-1, A( J+1, J ), 1, A( PVT, J+1 ), LDA )
343 * Swap dot products and PIV
346 WORK( J ) = WORK( PVT )
349 PIV( PVT ) = PIV( J )
356 * Compute elements J+1:N of column J
359 CALL SGEMV( 'No Trans', N-J, J-1, -ONE, A( J+1, 1 ), LDA,
360 $ A( J, 1 ), LDA, ONE, A( J+1, J ), 1 )
361 CALL SSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
368 * Ran to completion, A has full rank
375 * Rank is number of steps completed. Set INFO = 1 to signal
376 * that the factorization cannot be used to solve a system.