LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dpstf2.f
Go to the documentation of this file.
1*> \brief \b DPSTF2 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
9*> Download DPSTF2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dpstf2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dpstf2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dpstf2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DPSTF2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
22*
23* .. Scalar Arguments ..
24* DOUBLE PRECISION TOL
25* INTEGER INFO, LDA, N, RANK
26* CHARACTER UPLO
27* ..
28* .. Array Arguments ..
29* DOUBLE PRECISION A( LDA, * ), WORK( 2*N )
30* INTEGER PIV( N )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> DPSTF2 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 2 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 DOUBLE PRECISION 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[out] PIV
86*> \verbatim
87*> PIV is INTEGER array, dimension (N)
88*> PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
89*> \endverbatim
90*>
91*> \param[out] RANK
92*> \verbatim
93*> RANK is INTEGER
94*> The rank of A given by the number of steps the algorithm
95*> completed.
96*> \endverbatim
97*>
98*> \param[in] TOL
99*> \verbatim
100*> TOL is DOUBLE PRECISION
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.
104*> \endverbatim
105*>
106*> \param[in] LDA
107*> \verbatim
108*> LDA is INTEGER
109*> The leading dimension of the array A. LDA >= max(1,N).
110*> \endverbatim
111*>
112*> \param[out] WORK
113*> \verbatim
114*> WORK is DOUBLE PRECISION 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 pstf2
138*
139* =====================================================================
140 SUBROUTINE dpstf2( 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 DOUBLE PRECISION TOL
148 INTEGER INFO, LDA, N, RANK
149 CHARACTER UPLO
150* ..
151* .. Array Arguments ..
152 DOUBLE PRECISION A( LDA, * ), WORK( 2*N )
153 INTEGER PIV( N )
154* ..
155*
156* =====================================================================
157*
158* .. Parameters ..
159 DOUBLE PRECISION ONE, ZERO
160 parameter( one = 1.0d+0, zero = 0.0d+0 )
161* ..
162* .. Local Scalars ..
163 DOUBLE PRECISION AJJ, DSTOP, DTEMP
164 INTEGER I, ITEMP, J, PVT
165 LOGICAL UPPER
166* ..
167* .. External Functions ..
168 DOUBLE PRECISION DLAMCH
169 LOGICAL LSAME, DISNAN
170 EXTERNAL dlamch, lsame, disnan
171* ..
172* .. External Subroutines ..
173 EXTERNAL dgemv, dscal, dswap, 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( 'DPSTF2', -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.disnan( 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 dstop = n * dlamch( 'Epsilon' ) * ajj
227 ELSE
228 dstop = 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.dstop.OR.disnan( 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 dswap( j-1, a( 1, j ), 1, a( 1, pvt ), 1 )
272 IF( pvt.LT.n )
273 $ CALL dswap( n-pvt, a( j, pvt+1 ), lda,
274 $ a( pvt, pvt+1 ), lda )
275 CALL dswap( pvt-j-1, a( j, j+1 ), lda, a( j+1, pvt ), 1 )
276*
277* Swap dot products and PIV
278*
279 dtemp = work( j )
280 work( j ) = work( pvt )
281 work( pvt ) = dtemp
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 dgemv( 'Trans', j-1, n-j, -one, a( 1, j+1 ), lda,
294 $ a( 1, j ), 1, one, a( j, j+1 ), lda )
295 CALL dscal( 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.dstop.OR.disnan( 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 dswap( j-1, a( j, 1 ), lda, a( pvt, 1 ), lda )
335 IF( pvt.LT.n )
336 $ CALL dswap( n-pvt, a( pvt+1, j ), 1, a( pvt+1, pvt ),
337 $ 1 )
338 CALL dswap( pvt-j-1, a( j+1, j ), 1, a( pvt, j+1 ), lda )
339*
340* Swap dot products and PIV
341*
342 dtemp = work( j )
343 work( j ) = work( pvt )
344 work( pvt ) = dtemp
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 dgemv( '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 dscal( 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 DPSTF2
382*
383 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
subroutine dpstf2(uplo, n, a, lda, piv, rank, tol, work, info)
DPSTF2 computes the Cholesky factorization with complete pivoting of a real symmetric positive semide...
Definition dpstf2.f:141
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
subroutine dswap(n, dx, incx, dy, incy)
DSWAP
Definition dswap.f:82