LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dpstrf.f
Go to the documentation of this file.
00001 *> \brief \b DPSTRF
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DPSTRF + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dpstrf.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dpstrf.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dpstrf.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DPSTRF( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       DOUBLE PRECISION   TOL
00025 *       INTEGER            INFO, LDA, N, RANK
00026 *       CHARACTER          UPLO
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       DOUBLE PRECISION   A( LDA, * ), WORK( 2*N )
00030 *       INTEGER            PIV( N )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> DPSTRF computes the Cholesky factorization with complete
00040 *> pivoting of a real symmetric positive semidefinite matrix A.
00041 *>
00042 *> The factorization has the form
00043 *>    P**T * A * P = U**T * U ,  if UPLO = 'U',
00044 *>    P**T * A * P = L  * L**T,  if UPLO = 'L',
00045 *> where U is an upper triangular matrix and L is lower triangular, and
00046 *> P is stored as vector PIV.
00047 *>
00048 *> This algorithm does not attempt to check that A is positive
00049 *> semidefinite. This version of the algorithm calls level 3 BLAS.
00050 *> \endverbatim
00051 *
00052 *  Arguments:
00053 *  ==========
00054 *
00055 *> \param[in] UPLO
00056 *> \verbatim
00057 *>          UPLO is CHARACTER*1
00058 *>          Specifies whether the upper or lower triangular part of the
00059 *>          symmetric matrix A is stored.
00060 *>          = 'U':  Upper triangular
00061 *>          = 'L':  Lower triangular
00062 *> \endverbatim
00063 *>
00064 *> \param[in] N
00065 *> \verbatim
00066 *>          N is INTEGER
00067 *>          The order of the matrix A.  N >= 0.
00068 *> \endverbatim
00069 *>
00070 *> \param[in,out] A
00071 *> \verbatim
00072 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
00073 *>          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
00074 *>          n by n upper triangular part of A contains the upper
00075 *>          triangular part of the matrix A, and the strictly lower
00076 *>          triangular part of A is not referenced.  If UPLO = 'L', the
00077 *>          leading n by n lower triangular part of A contains the lower
00078 *>          triangular part of the matrix A, and the strictly upper
00079 *>          triangular part of A is not referenced.
00080 *>
00081 *>          On exit, if INFO = 0, the factor U or L from the Cholesky
00082 *>          factorization as above.
00083 *> \endverbatim
00084 *>
00085 *> \param[in] LDA
00086 *> \verbatim
00087 *>          LDA is INTEGER
00088 *>          The leading dimension of the array A.  LDA >= max(1,N).
00089 *> \endverbatim
00090 *>
00091 *> \param[out] PIV
00092 *> \verbatim
00093 *>          PIV is INTEGER array, dimension (N)
00094 *>          PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
00095 *> \endverbatim
00096 *>
00097 *> \param[out] RANK
00098 *> \verbatim
00099 *>          RANK is INTEGER
00100 *>          The rank of A given by the number of steps the algorithm
00101 *>          completed.
00102 *> \endverbatim
00103 *>
00104 *> \param[in] TOL
00105 *> \verbatim
00106 *>          TOL is DOUBLE PRECISION
00107 *>          User defined tolerance. If TOL < 0, then N*U*MAX( A(K,K) )
00108 *>          will be used. The algorithm terminates at the (K-1)st step
00109 *>          if the pivot <= TOL.
00110 *> \endverbatim
00111 *>
00112 *> \param[out] WORK
00113 *> \verbatim
00114 *>          WORK is DOUBLE PRECISION array, dimension (2*N)
00115 *>          Work space.
00116 *> \endverbatim
00117 *>
00118 *> \param[out] INFO
00119 *> \verbatim
00120 *>          INFO is INTEGER
00121 *>          < 0: If INFO = -K, the K-th argument had an illegal value,
00122 *>          = 0: algorithm completed successfully, and
00123 *>          > 0: the matrix A is either rank deficient with computed rank
00124 *>               as returned in RANK, or is indefinite.  See Section 7 of
00125 *>               LAPACK Working Note #161 for further information.
00126 *> \endverbatim
00127 *
00128 *  Authors:
00129 *  ========
00130 *
00131 *> \author Univ. of Tennessee 
00132 *> \author Univ. of California Berkeley 
00133 *> \author Univ. of Colorado Denver 
00134 *> \author NAG Ltd. 
00135 *
00136 *> \date November 2011
00137 *
00138 *> \ingroup doubleOTHERcomputational
00139 *
00140 *  =====================================================================
00141       SUBROUTINE DPSTRF( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
00142 *
00143 *  -- LAPACK computational routine (version 3.4.0) --
00144 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00145 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00146 *     November 2011
00147 *
00148 *     .. Scalar Arguments ..
00149       DOUBLE PRECISION   TOL
00150       INTEGER            INFO, LDA, N, RANK
00151       CHARACTER          UPLO
00152 *     ..
00153 *     .. Array Arguments ..
00154       DOUBLE PRECISION   A( LDA, * ), WORK( 2*N )
00155       INTEGER            PIV( N )
00156 *     ..
00157 *
00158 *  =====================================================================
00159 *
00160 *     .. Parameters ..
00161       DOUBLE PRECISION   ONE, ZERO
00162       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00163 *     ..
00164 *     .. Local Scalars ..
00165       DOUBLE PRECISION   AJJ, DSTOP, DTEMP
00166       INTEGER            I, ITEMP, J, JB, K, NB, PVT
00167       LOGICAL            UPPER
00168 *     ..
00169 *     .. External Functions ..
00170       DOUBLE PRECISION   DLAMCH
00171       INTEGER            ILAENV
00172       LOGICAL            LSAME, DISNAN
00173       EXTERNAL           DLAMCH, ILAENV, LSAME, DISNAN
00174 *     ..
00175 *     .. External Subroutines ..
00176       EXTERNAL           DGEMV, DPSTF2, DSCAL, DSWAP, DSYRK, XERBLA
00177 *     ..
00178 *     .. Intrinsic Functions ..
00179       INTRINSIC          MAX, MIN, SQRT, MAXLOC
00180 *     ..
00181 *     .. Executable Statements ..
00182 *
00183 *     Test the input parameters.
00184 *
00185       INFO = 0
00186       UPPER = LSAME( UPLO, 'U' )
00187       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00188          INFO = -1
00189       ELSE IF( N.LT.0 ) THEN
00190          INFO = -2
00191       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00192          INFO = -4
00193       END IF
00194       IF( INFO.NE.0 ) THEN
00195          CALL XERBLA( 'DPSTRF', -INFO )
00196          RETURN
00197       END IF
00198 *
00199 *     Quick return if possible
00200 *
00201       IF( N.EQ.0 )
00202      $   RETURN
00203 *
00204 *     Get block size
00205 *
00206       NB = ILAENV( 1, 'DPOTRF', UPLO, N, -1, -1, -1 )
00207       IF( NB.LE.1 .OR. NB.GE.N ) THEN
00208 *
00209 *        Use unblocked code
00210 *
00211          CALL DPSTF2( UPLO, N, A( 1, 1 ), LDA, PIV, RANK, TOL, WORK,
00212      $                INFO )
00213          GO TO 200
00214 *
00215       ELSE
00216 *
00217 *     Initialize PIV
00218 *
00219          DO 100 I = 1, N
00220             PIV( I ) = I
00221   100    CONTINUE
00222 *
00223 *     Compute stopping value
00224 *
00225          PVT = 1
00226          AJJ = A( PVT, PVT )
00227          DO I = 2, N
00228             IF( A( I, I ).GT.AJJ ) THEN
00229                PVT = I
00230                AJJ = A( PVT, PVT )
00231             END IF
00232          END DO
00233          IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
00234             RANK = 0
00235             INFO = 1
00236             GO TO 200
00237          END IF
00238 *
00239 *     Compute stopping value if not supplied
00240 *
00241          IF( TOL.LT.ZERO ) THEN
00242             DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ
00243          ELSE
00244             DSTOP = TOL
00245          END IF
00246 *
00247 *
00248          IF( UPPER ) THEN
00249 *
00250 *           Compute the Cholesky factorization P**T * A * P = U**T * U
00251 *
00252             DO 140 K = 1, N, NB
00253 *
00254 *              Account for last block not being NB wide
00255 *
00256                JB = MIN( NB, N-K+1 )
00257 *
00258 *              Set relevant part of first half of WORK to zero,
00259 *              holds dot products
00260 *
00261                DO 110 I = K, N
00262                   WORK( I ) = 0
00263   110          CONTINUE
00264 *
00265                DO 130 J = K, K + JB - 1
00266 *
00267 *              Find pivot, test for exit, else swap rows and columns
00268 *              Update dot products, compute possible pivots which are
00269 *              stored in the second half of WORK
00270 *
00271                   DO 120 I = J, N
00272 *
00273                      IF( J.GT.K ) THEN
00274                         WORK( I ) = WORK( I ) + A( J-1, I )**2
00275                      END IF
00276                      WORK( N+I ) = A( I, I ) - WORK( I )
00277 *
00278   120             CONTINUE
00279 *
00280                   IF( J.GT.1 ) THEN
00281                      ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
00282                      PVT = ITEMP + J - 1
00283                      AJJ = WORK( N+PVT )
00284                      IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
00285                         A( J, J ) = AJJ
00286                         GO TO 190
00287                      END IF
00288                   END IF
00289 *
00290                   IF( J.NE.PVT ) THEN
00291 *
00292 *                    Pivot OK, so can now swap pivot rows and columns
00293 *
00294                      A( PVT, PVT ) = A( J, J )
00295                      CALL DSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
00296                      IF( PVT.LT.N )
00297      $                  CALL DSWAP( N-PVT, A( J, PVT+1 ), LDA,
00298      $                              A( PVT, PVT+1 ), LDA )
00299                      CALL DSWAP( PVT-J-1, A( J, J+1 ), LDA,
00300      $                           A( J+1, PVT ), 1 )
00301 *
00302 *                    Swap dot products and PIV
00303 *
00304                      DTEMP = WORK( J )
00305                      WORK( J ) = WORK( PVT )
00306                      WORK( PVT ) = DTEMP
00307                      ITEMP = PIV( PVT )
00308                      PIV( PVT ) = PIV( J )
00309                      PIV( J ) = ITEMP
00310                   END IF
00311 *
00312                   AJJ = SQRT( AJJ )
00313                   A( J, J ) = AJJ
00314 *
00315 *                 Compute elements J+1:N of row J.
00316 *
00317                   IF( J.LT.N ) THEN
00318                      CALL DGEMV( 'Trans', J-K, N-J, -ONE, A( K, J+1 ),
00319      $                           LDA, A( K, J ), 1, ONE, A( J, J+1 ),
00320      $                           LDA )
00321                      CALL DSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
00322                   END IF
00323 *
00324   130          CONTINUE
00325 *
00326 *              Update trailing matrix, J already incremented
00327 *
00328                IF( K+JB.LE.N ) THEN
00329                   CALL DSYRK( 'Upper', 'Trans', N-J+1, JB, -ONE,
00330      $                        A( K, J ), LDA, ONE, A( J, J ), LDA )
00331                END IF
00332 *
00333   140       CONTINUE
00334 *
00335          ELSE
00336 *
00337 *        Compute the Cholesky factorization P**T * A * P = L * L**T
00338 *
00339             DO 180 K = 1, N, NB
00340 *
00341 *              Account for last block not being NB wide
00342 *
00343                JB = MIN( NB, N-K+1 )
00344 *
00345 *              Set relevant part of first half of WORK to zero,
00346 *              holds dot products
00347 *
00348                DO 150 I = K, N
00349                   WORK( I ) = 0
00350   150          CONTINUE
00351 *
00352                DO 170 J = K, K + JB - 1
00353 *
00354 *              Find pivot, test for exit, else swap rows and columns
00355 *              Update dot products, compute possible pivots which are
00356 *              stored in the second half of WORK
00357 *
00358                   DO 160 I = J, N
00359 *
00360                      IF( J.GT.K ) THEN
00361                         WORK( I ) = WORK( I ) + A( I, J-1 )**2
00362                      END IF
00363                      WORK( N+I ) = A( I, I ) - WORK( I )
00364 *
00365   160             CONTINUE
00366 *
00367                   IF( J.GT.1 ) THEN
00368                      ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
00369                      PVT = ITEMP + J - 1
00370                      AJJ = WORK( N+PVT )
00371                      IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
00372                         A( J, J ) = AJJ
00373                         GO TO 190
00374                      END IF
00375                   END IF
00376 *
00377                   IF( J.NE.PVT ) THEN
00378 *
00379 *                    Pivot OK, so can now swap pivot rows and columns
00380 *
00381                      A( PVT, PVT ) = A( J, J )
00382                      CALL DSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
00383                      IF( PVT.LT.N )
00384      $                  CALL DSWAP( N-PVT, A( PVT+1, J ), 1,
00385      $                              A( PVT+1, PVT ), 1 )
00386                      CALL DSWAP( PVT-J-1, A( J+1, J ), 1, A( PVT, J+1 ),
00387      $                           LDA )
00388 *
00389 *                    Swap dot products and PIV
00390 *
00391                      DTEMP = WORK( J )
00392                      WORK( J ) = WORK( PVT )
00393                      WORK( PVT ) = DTEMP
00394                      ITEMP = PIV( PVT )
00395                      PIV( PVT ) = PIV( J )
00396                      PIV( J ) = ITEMP
00397                   END IF
00398 *
00399                   AJJ = SQRT( AJJ )
00400                   A( J, J ) = AJJ
00401 *
00402 *                 Compute elements J+1:N of column J.
00403 *
00404                   IF( J.LT.N ) THEN
00405                      CALL DGEMV( 'No Trans', N-J, J-K, -ONE,
00406      $                           A( J+1, K ), LDA, A( J, K ), LDA, ONE,
00407      $                           A( J+1, J ), 1 )
00408                      CALL DSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
00409                   END IF
00410 *
00411   170          CONTINUE
00412 *
00413 *              Update trailing matrix, J already incremented
00414 *
00415                IF( K+JB.LE.N ) THEN
00416                   CALL DSYRK( 'Lower', 'No Trans', N-J+1, JB, -ONE,
00417      $                        A( J, K ), LDA, ONE, A( J, J ), LDA )
00418                END IF
00419 *
00420   180       CONTINUE
00421 *
00422          END IF
00423       END IF
00424 *
00425 *     Ran to completion, A has full rank
00426 *
00427       RANK = N
00428 *
00429       GO TO 200
00430   190 CONTINUE
00431 *
00432 *     Rank is the number of steps completed.  Set INFO = 1 to signal
00433 *     that the factorization cannot be used to solve a system.
00434 *
00435       RANK = J - 1
00436       INFO = 1
00437 *
00438   200 CONTINUE
00439       RETURN
00440 *
00441 *     End of DPSTRF
00442 *
00443       END
 All Files Functions