LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ spstf2()

 subroutine spstf2 ( character uplo, integer n, real, dimension( lda, * ) a, integer lda, integer, dimension( n ) piv, integer rank, real tol, real, dimension( 2*n ) work, integer info )

SPSTF2 computes the Cholesky factorization with complete pivoting of a real symmetric positive semidefinite matrix.

Purpose:
``` SPSTF2 computes the Cholesky factorization with complete
pivoting of a real symmetric positive semidefinite matrix A.

The factorization has the form
P**T * A * P = U**T * U ,  if UPLO = 'U',
P**T * A * P = L  * L**T,  if UPLO = 'L',
where U is an upper triangular matrix and L is lower triangular, and
P is stored as vector PIV.

This algorithm does not attempt to check that A is positive
semidefinite. This version of the algorithm calls level 2 BLAS.```
Parameters
 [in] UPLO ``` UPLO is CHARACTER*1 Specifies whether the upper or lower triangular part of the symmetric matrix A is stored. = 'U': Upper triangular = 'L': Lower triangular``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in,out] A ``` A is REAL array, dimension (LDA,N) On entry, the symmetric matrix A. If UPLO = 'U', the leading n by n upper triangular part of A contains the upper triangular part of the matrix A, and the strictly lower triangular part of A is not referenced. If UPLO = 'L', the leading n by n lower triangular part of A contains the lower triangular part of the matrix A, and the strictly upper triangular part of A is not referenced. On exit, if INFO = 0, the factor U or L from the Cholesky factorization as above.``` [out] PIV ``` PIV is INTEGER array, dimension (N) PIV is such that the nonzero entries are P( PIV(K), K ) = 1.``` [out] RANK ``` RANK is INTEGER The rank of A given by the number of steps the algorithm completed.``` [in] TOL ``` TOL is REAL User defined tolerance. If TOL < 0, then N*U*MAX( A( K,K ) ) will be used. The algorithm terminates at the (K-1)st step if the pivot <= TOL.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).``` [out] WORK ``` WORK is REAL array, dimension (2*N) Work space.``` [out] INFO ``` INFO is INTEGER < 0: If INFO = -K, the K-th argument had an illegal value, = 0: algorithm completed successfully, and > 0: the matrix A is either rank deficient with computed rank as returned in RANK, or is not positive semidefinite. See Section 7 of LAPACK Working Note #161 for further information.```

Definition at line 140 of file spstf2.f.

141*
142* -- LAPACK computational routine --
143* -- LAPACK is a software package provided by Univ. of Tennessee, --
144* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
145*
146* .. Scalar Arguments ..
147 REAL TOL
148 INTEGER INFO, LDA, N, RANK
149 CHARACTER UPLO
150* ..
151* .. Array Arguments ..
152 REAL A( LDA, * ), WORK( 2*N )
153 INTEGER PIV( N )
154* ..
155*
156* =====================================================================
157*
158* .. Parameters ..
159 REAL ONE, ZERO
160 parameter( one = 1.0e+0, zero = 0.0e+0 )
161* ..
162* .. Local Scalars ..
163 REAL AJJ, SSTOP, STEMP
164 INTEGER I, ITEMP, J, PVT
165 LOGICAL UPPER
166* ..
167* .. External Functions ..
168 REAL SLAMCH
169 LOGICAL LSAME, SISNAN
170 EXTERNAL slamch, lsame, sisnan
171* ..
172* .. External Subroutines ..
173 EXTERNAL sgemv, sscal, sswap, xerbla
174* ..
175* .. Intrinsic Functions ..
176 INTRINSIC max, sqrt, maxloc
177* ..
178* .. Executable Statements ..
179*
180* Test the input parameters
181*
182 info = 0
183 upper = lsame( uplo, 'U' )
184 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
185 info = -1
186 ELSE IF( n.LT.0 ) THEN
187 info = -2
188 ELSE IF( lda.LT.max( 1, n ) ) THEN
189 info = -4
190 END IF
191 IF( info.NE.0 ) THEN
192 CALL xerbla( 'SPSTF2', -info )
193 RETURN
194 END IF
195*
196* Quick return if possible
197*
198 IF( n.EQ.0 )
199 \$ RETURN
200*
201* Initialize PIV
202*
203 DO 100 i = 1, n
204 piv( i ) = i
205 100 CONTINUE
206*
207* Compute stopping value
208*
209 pvt = 1
210 ajj = a( pvt, pvt )
211 DO i = 2, n
212 IF( a( i, i ).GT.ajj ) THEN
213 pvt = i
214 ajj = a( pvt, pvt )
215 END IF
216 END DO
217 IF( ajj.LE.zero.OR.sisnan( ajj ) ) THEN
218 rank = 0
219 info = 1
220 GO TO 170
221 END IF
222*
223* Compute stopping value if not supplied
224*
225 IF( tol.LT.zero ) THEN
226 sstop = n * slamch( 'Epsilon' ) * ajj
227 ELSE
228 sstop = tol
229 END IF
230*
231* Set first half of WORK to zero, holds dot products
232*
233 DO 110 i = 1, n
234 work( i ) = 0
235 110 CONTINUE
236*
237 IF( upper ) THEN
238*
239* Compute the Cholesky factorization P**T * A * P = U**T * U
240*
241 DO 130 j = 1, n
242*
243* Find pivot, test for exit, else swap rows and columns
244* Update dot products, compute possible pivots which are
245* stored in the second half of WORK
246*
247 DO 120 i = j, n
248*
249 IF( j.GT.1 ) THEN
250 work( i ) = work( i ) + a( j-1, i )**2
251 END IF
252 work( n+i ) = a( i, i ) - work( i )
253*
254 120 CONTINUE
255*
256 IF( j.GT.1 ) THEN
257 itemp = maxloc( work( (n+j):(2*n) ), 1 )
258 pvt = itemp + j - 1
259 ajj = work( n+pvt )
260 IF( ajj.LE.sstop.OR.sisnan( ajj ) ) THEN
261 a( j, j ) = ajj
262 GO TO 160
263 END IF
264 END IF
265*
266 IF( j.NE.pvt ) THEN
267*
268* Pivot OK, so can now swap pivot rows and columns
269*
270 a( pvt, pvt ) = a( j, j )
271 CALL sswap( j-1, a( 1, j ), 1, a( 1, pvt ), 1 )
272 IF( pvt.LT.n )
273 \$ CALL sswap( n-pvt, a( j, pvt+1 ), lda,
274 \$ a( pvt, pvt+1 ), lda )
275 CALL sswap( pvt-j-1, a( j, j+1 ), lda, a( j+1, pvt ), 1 )
276*
277* Swap dot products and PIV
278*
279 stemp = work( j )
280 work( j ) = work( pvt )
281 work( pvt ) = stemp
282 itemp = piv( pvt )
283 piv( pvt ) = piv( j )
284 piv( j ) = itemp
285 END IF
286*
287 ajj = sqrt( ajj )
288 a( j, j ) = ajj
289*
290* Compute elements J+1:N of row J
291*
292 IF( j.LT.n ) THEN
293 CALL sgemv( 'Trans', j-1, n-j, -one, a( 1, j+1 ), lda,
294 \$ a( 1, j ), 1, one, a( j, j+1 ), lda )
295 CALL sscal( n-j, one / ajj, a( j, j+1 ), lda )
296 END IF
297*
298 130 CONTINUE
299*
300 ELSE
301*
302* Compute the Cholesky factorization P**T * A * P = L * L**T
303*
304 DO 150 j = 1, n
305*
306* Find pivot, test for exit, else swap rows and columns
307* Update dot products, compute possible pivots which are
308* stored in the second half of WORK
309*
310 DO 140 i = j, n
311*
312 IF( j.GT.1 ) THEN
313 work( i ) = work( i ) + a( i, j-1 )**2
314 END IF
315 work( n+i ) = a( i, i ) - work( i )
316*
317 140 CONTINUE
318*
319 IF( j.GT.1 ) THEN
320 itemp = maxloc( work( (n+j):(2*n) ), 1 )
321 pvt = itemp + j - 1
322 ajj = work( n+pvt )
323 IF( ajj.LE.sstop.OR.sisnan( ajj ) ) THEN
324 a( j, j ) = ajj
325 GO TO 160
326 END IF
327 END IF
328*
329 IF( j.NE.pvt ) THEN
330*
331* Pivot OK, so can now swap pivot rows and columns
332*
333 a( pvt, pvt ) = a( j, j )
334 CALL sswap( j-1, a( j, 1 ), lda, a( pvt, 1 ), lda )
335 IF( pvt.LT.n )
336 \$ CALL sswap( n-pvt, a( pvt+1, j ), 1, a( pvt+1, pvt ),
337 \$ 1 )
338 CALL sswap( pvt-j-1, a( j+1, j ), 1, a( pvt, j+1 ), lda )
339*
340* Swap dot products and PIV
341*
342 stemp = work( j )
343 work( j ) = work( pvt )
344 work( pvt ) = stemp
345 itemp = piv( pvt )
346 piv( pvt ) = piv( j )
347 piv( j ) = itemp
348 END IF
349*
350 ajj = sqrt( ajj )
351 a( j, j ) = ajj
352*
353* Compute elements J+1:N of column J
354*
355 IF( j.LT.n ) THEN
356 CALL sgemv( 'No Trans', n-j, j-1, -one, a( j+1, 1 ), lda,
357 \$ a( j, 1 ), lda, one, a( j+1, j ), 1 )
358 CALL sscal( n-j, one / ajj, a( j+1, j ), 1 )
359 END IF
360*
361 150 CONTINUE
362*
363 END IF
364*
365* Ran to completion, A has full rank
366*
367 rank = n
368*
369 GO TO 170
370 160 CONTINUE
371*
372* Rank is number of steps completed. Set INFO = 1 to signal
373* that the factorization cannot be used to solve a system.
374*
375 rank = j - 1
376 info = 1
377*
378 170 CONTINUE
379 RETURN
380*
381* End of SPSTF2
382*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:158
logical function sisnan(sin)
SISNAN tests input for NaN.
Definition sisnan.f:59
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82
Here is the call graph for this function:
Here is the caller graph for this function: