LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches
spstrf.f
Go to the documentation of this file.
1*> \brief \b SPSTRF computes the Cholesky factorization with complete pivoting of a real symmetric positive semidefinite matrix.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/spstrf.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/spstrf.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/spstrf.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SPSTRF( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
22*
23* .. Scalar Arguments ..
24* REAL TOL
25* INTEGER INFO, LDA, N, RANK
26* CHARACTER UPLO
27* ..
28* .. Array Arguments ..
29* REAL A( LDA, * ), WORK( 2*N )
30* INTEGER PIV( N )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> SPSTRF computes the Cholesky factorization with complete
40*> pivoting of a real symmetric positive semidefinite matrix A.
41*>
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.
47*>
48*> This algorithm does not attempt to check that A is positive
49*> semidefinite. This version of the algorithm calls level 3 BLAS.
50*> \endverbatim
51*
52* Arguments:
53* ==========
54*
55*> \param[in] UPLO
56*> \verbatim
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
62*> \endverbatim
63*>
64*> \param[in] N
65*> \verbatim
66*> N is INTEGER
67*> The order of the matrix A. N >= 0.
68*> \endverbatim
69*>
70*> \param[in,out] A
71*> \verbatim
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.
80*>
81*> On exit, if INFO = 0, the factor U or L from the Cholesky
82*> factorization as above.
83*> \endverbatim
84*>
85*> \param[in] LDA
86*> \verbatim
87*> LDA is INTEGER
88*> The leading dimension of the array A. LDA >= max(1,N).
89*> \endverbatim
90*>
91*> \param[out] PIV
92*> \verbatim
93*> PIV is INTEGER array, dimension (N)
94*> PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
95*> \endverbatim
96*>
97*> \param[out] RANK
98*> \verbatim
99*> RANK is INTEGER
100*> The rank of A given by the number of steps the algorithm
101*> completed.
102*> \endverbatim
103*>
104*> \param[in] TOL
105*> \verbatim
106*> TOL is REAL
107*> User defined tolerance. If TOL < 0, then N*U*MAX( A(K,K) )
108*> will be used. The algorithm terminates at the (K-1)st step
109*> if the pivot <= TOL.
110*> \endverbatim
111*>
112*> \param[out] WORK
113*> \verbatim
114*> WORK is REAL array, dimension (2*N)
115*> Work space.
116*> \endverbatim
117*>
118*> \param[out] INFO
119*> \verbatim
120*> INFO is INTEGER
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
126*> information.
127*> \endverbatim
128*
129* Authors:
130* ========
131*
132*> \author Univ. of Tennessee
133*> \author Univ. of California Berkeley
134*> \author Univ. of Colorado Denver
135*> \author NAG Ltd.
136*
137*> \ingroup pstrf
138*
139* =====================================================================
140 SUBROUTINE spstrf( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
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, JB, K, NB, PVT
165 LOGICAL UPPER
166* ..
167* .. External Functions ..
168 REAL SLAMCH
169 INTEGER ILAENV
170 LOGICAL LSAME, SISNAN
171 EXTERNAL slamch, ilaenv, lsame, sisnan
172* ..
173* .. External Subroutines ..
174 EXTERNAL sgemv, spstf2, sscal, sswap, ssyrk, xerbla
175* ..
176* .. Intrinsic Functions ..
177 INTRINSIC max, min, sqrt, maxloc
178* ..
179* .. Executable Statements ..
180*
181* Test the input parameters.
182*
183 info = 0
184 upper = lsame( uplo, 'U' )
185 IF( .NOT.upper .AND. .NOT.lsame( uplo, 'L' ) ) THEN
186 info = -1
187 ELSE IF( n.LT.0 ) THEN
188 info = -2
189 ELSE IF( lda.LT.max( 1, n ) ) THEN
190 info = -4
191 END IF
192 IF( info.NE.0 ) THEN
193 CALL xerbla( 'SPSTRF', -info )
194 RETURN
195 END IF
196*
197* Quick return if possible
198*
199 IF( n.EQ.0 )
200 \$ RETURN
201*
202* Get block size
203*
204 nb = ilaenv( 1, 'SPOTRF', uplo, n, -1, -1, -1 )
205 IF( nb.LE.1 .OR. nb.GE.n ) THEN
206*
207* Use unblocked code
208*
209 CALL spstf2( uplo, n, a( 1, 1 ), lda, piv, rank, tol, work,
210 \$ info )
211 GO TO 200
212*
213 ELSE
214*
215* Initialize PIV
216*
217 DO 100 i = 1, n
218 piv( i ) = i
219 100 CONTINUE
220*
221* Compute stopping value
222*
223 pvt = 1
224 ajj = a( pvt, pvt )
225 DO i = 2, n
226 IF( a( i, i ).GT.ajj ) THEN
227 pvt = i
228 ajj = a( pvt, pvt )
229 END IF
230 END DO
231 IF( ajj.LE.zero.OR.sisnan( ajj ) ) THEN
232 rank = 0
233 info = 1
234 GO TO 200
235 END IF
236*
237* Compute stopping value if not supplied
238*
239 IF( tol.LT.zero ) THEN
240 sstop = n * slamch( 'Epsilon' ) * ajj
241 ELSE
242 sstop = tol
243 END IF
244*
245*
246 IF( upper ) THEN
247*
248* Compute the Cholesky factorization P**T * A * P = U**T * U
249*
250 DO 140 k = 1, n, nb
251*
252* Account for last block not being NB wide
253*
254 jb = min( nb, n-k+1 )
255*
256* Set relevant part of first half of WORK to zero,
257* holds dot products
258*
259 DO 110 i = k, n
260 work( i ) = 0
261 110 CONTINUE
262*
263 DO 130 j = k, k + jb - 1
264*
265* Find pivot, test for exit, else swap rows and columns
266* Update dot products, compute possible pivots which are
267* stored in the second half of WORK
268*
269 DO 120 i = j, n
270*
271 IF( j.GT.k ) THEN
272 work( i ) = work( i ) + a( j-1, i )**2
273 END IF
274 work( n+i ) = a( i, i ) - work( i )
275*
276 120 CONTINUE
277*
278 IF( j.GT.1 ) THEN
279 itemp = maxloc( work( (n+j):(2*n) ), 1 )
280 pvt = itemp + j - 1
281 ajj = work( n+pvt )
282 IF( ajj.LE.sstop.OR.sisnan( ajj ) ) THEN
283 a( j, j ) = ajj
284 GO TO 190
285 END IF
286 END IF
287*
288 IF( j.NE.pvt ) THEN
289*
290* Pivot OK, so can now swap pivot rows and columns
291*
292 a( pvt, pvt ) = a( j, j )
293 CALL sswap( j-1, a( 1, j ), 1, a( 1, pvt ), 1 )
294 IF( pvt.LT.n )
295 \$ CALL sswap( n-pvt, a( j, pvt+1 ), lda,
296 \$ a( pvt, pvt+1 ), lda )
297 CALL sswap( pvt-j-1, a( j, j+1 ), lda,
298 \$ a( j+1, pvt ), 1 )
299*
300* Swap dot products and PIV
301*
302 stemp = work( j )
303 work( j ) = work( pvt )
304 work( pvt ) = stemp
305 itemp = piv( pvt )
306 piv( pvt ) = piv( j )
307 piv( j ) = itemp
308 END IF
309*
310 ajj = sqrt( ajj )
311 a( j, j ) = ajj
312*
313* Compute elements J+1:N of row J.
314*
315 IF( j.LT.n ) THEN
316 CALL sgemv( 'Trans', j-k, n-j, -one, a( k, j+1 ),
317 \$ lda, a( k, j ), 1, one, a( j, j+1 ),
318 \$ lda )
319 CALL sscal( n-j, one / ajj, a( j, j+1 ), lda )
320 END IF
321*
322 130 CONTINUE
323*
324* Update trailing matrix, J already incremented
325*
326 IF( k+jb.LE.n ) THEN
327 CALL ssyrk( 'Upper', 'Trans', n-j+1, jb, -one,
328 \$ a( k, j ), lda, one, a( j, j ), lda )
329 END IF
330*
331 140 CONTINUE
332*
333 ELSE
334*
335* Compute the Cholesky factorization P**T * A * P = L * L**T
336*
337 DO 180 k = 1, n, nb
338*
339* Account for last block not being NB wide
340*
341 jb = min( nb, n-k+1 )
342*
343* Set relevant part of first half of WORK to zero,
344* holds dot products
345*
346 DO 150 i = k, n
347 work( i ) = 0
348 150 CONTINUE
349*
350 DO 170 j = k, k + jb - 1
351*
352* Find pivot, test for exit, else swap rows and columns
353* Update dot products, compute possible pivots which are
354* stored in the second half of WORK
355*
356 DO 160 i = j, n
357*
358 IF( j.GT.k ) THEN
359 work( i ) = work( i ) + a( i, j-1 )**2
360 END IF
361 work( n+i ) = a( i, i ) - work( i )
362*
363 160 CONTINUE
364*
365 IF( j.GT.1 ) THEN
366 itemp = maxloc( work( (n+j):(2*n) ), 1 )
367 pvt = itemp + j - 1
368 ajj = work( n+pvt )
369 IF( ajj.LE.sstop.OR.sisnan( ajj ) ) THEN
370 a( j, j ) = ajj
371 GO TO 190
372 END IF
373 END IF
374*
375 IF( j.NE.pvt ) THEN
376*
377* Pivot OK, so can now swap pivot rows and columns
378*
379 a( pvt, pvt ) = a( j, j )
380 CALL sswap( j-1, a( j, 1 ), lda, a( pvt, 1 ), lda )
381 IF( pvt.LT.n )
382 \$ CALL sswap( n-pvt, a( pvt+1, j ), 1,
383 \$ a( pvt+1, pvt ), 1 )
384 CALL sswap( pvt-j-1, a( j+1, j ), 1, a( pvt, j+1 ),
385 \$ lda )
386*
387* Swap dot products and PIV
388*
389 stemp = work( j )
390 work( j ) = work( pvt )
391 work( pvt ) = stemp
392 itemp = piv( pvt )
393 piv( pvt ) = piv( j )
394 piv( j ) = itemp
395 END IF
396*
397 ajj = sqrt( ajj )
398 a( j, j ) = ajj
399*
400* Compute elements J+1:N of column J.
401*
402 IF( j.LT.n ) THEN
403 CALL sgemv( 'No Trans', n-j, j-k, -one,
404 \$ a( j+1, k ), lda, a( j, k ), lda, one,
405 \$ a( j+1, j ), 1 )
406 CALL sscal( n-j, one / ajj, a( j+1, j ), 1 )
407 END IF
408*
409 170 CONTINUE
410*
411* Update trailing matrix, J already incremented
412*
413 IF( k+jb.LE.n ) THEN
414 CALL ssyrk( 'Lower', 'No Trans', n-j+1, jb, -one,
415 \$ a( j, k ), lda, one, a( j, j ), lda )
416 END IF
417*
418 180 CONTINUE
419*
420 END IF
421 END IF
422*
423* Ran to completion, A has full rank
424*
425 rank = n
426*
427 GO TO 200
428 190 CONTINUE
429*
430* Rank is the number of steps completed. Set INFO = 1 to signal
431* that the factorization cannot be used to solve a system.
432*
433 rank = j - 1
434 info = 1
435*
436 200 CONTINUE
437 RETURN
438*
439* End of SPSTRF
440*
441 END
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
subroutine ssyrk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
SSYRK
Definition ssyrk.f:169
subroutine spstf2(uplo, n, a, lda, piv, rank, tol, work, info)
SPSTF2 computes the Cholesky factorization with complete pivoting of a real symmetric positive semide...
Definition spstf2.f:141
subroutine spstrf(uplo, n, a, lda, piv, rank, tol, work, info)
SPSTRF computes the Cholesky factorization with complete pivoting of a real symmetric positive semide...
Definition spstrf.f:141
subroutine sscal(n, sa, sx, incx)
SSCAL
Definition sscal.f:79
subroutine sswap(n, sx, incx, sy, incy)
SSWAP
Definition sswap.f:82