LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dpstf2()

subroutine dpstf2 ( character  UPLO,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
integer, dimension( n )  PIV,
integer  RANK,
double precision  TOL,
double precision, dimension( 2*n )  WORK,
integer  INFO 
)

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

Download DPSTF2 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DPSTF2 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 DOUBLE PRECISION 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 DOUBLE PRECISION
          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 DOUBLE PRECISION 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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 140 of file dpstf2.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 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*
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
logical function disnan(DIN)
DISNAN tests input for NaN.
Definition: disnan.f:59
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:82
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:156
Here is the call graph for this function:
Here is the caller graph for this function: