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