LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
spstf2.f
Go to the documentation of this file.
00001 *> \brief \b SPSTF2
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download SPSTF2 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/spstf2.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/spstf2.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/spstf2.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE SPSTF2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       REAL               TOL
00025 *       INTEGER            INFO, LDA, N, RANK
00026 *       CHARACTER          UPLO
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       REAL               A( LDA, * ), WORK( 2*N )
00030 *       INTEGER            PIV( N )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> SPSTF2 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 2 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 REAL 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[out] PIV
00086 *> \verbatim
00087 *>          PIV is INTEGER array, dimension (N)
00088 *>          PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
00089 *> \endverbatim
00090 *>
00091 *> \param[out] RANK
00092 *> \verbatim
00093 *>          RANK is INTEGER
00094 *>          The rank of A given by the number of steps the algorithm
00095 *>          completed.
00096 *> \endverbatim
00097 *>
00098 *> \param[in] TOL
00099 *> \verbatim
00100 *>          TOL is REAL
00101 *>          User defined tolerance. If TOL < 0, then N*U*MAX( A( K,K ) )
00102 *>          will be used. The algorithm terminates at the (K-1)st step
00103 *>          if the pivot <= TOL.
00104 *> \endverbatim
00105 *>
00106 *> \param[in] LDA
00107 *> \verbatim
00108 *>          LDA is INTEGER
00109 *>          The leading dimension of the array A.  LDA >= max(1,N).
00110 *> \endverbatim
00111 *>
00112 *> \param[out] WORK
00113 *> \verbatim
00114 *>          WORK is REAL 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 realOTHERcomputational
00139 *
00140 *  =====================================================================
00141       SUBROUTINE SPSTF2( 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       REAL               TOL
00150       INTEGER            INFO, LDA, N, RANK
00151       CHARACTER          UPLO
00152 *     ..
00153 *     .. Array Arguments ..
00154       REAL               A( LDA, * ), WORK( 2*N )
00155       INTEGER            PIV( N )
00156 *     ..
00157 *
00158 *  =====================================================================
00159 *
00160 *     .. Parameters ..
00161       REAL               ONE, ZERO
00162       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00163 *     ..
00164 *     .. Local Scalars ..
00165       REAL               AJJ, SSTOP, STEMP
00166       INTEGER            I, ITEMP, J, PVT
00167       LOGICAL            UPPER
00168 *     ..
00169 *     .. External Functions ..
00170       REAL               SLAMCH
00171       LOGICAL            LSAME, SISNAN
00172       EXTERNAL           SLAMCH, LSAME, SISNAN
00173 *     ..
00174 *     .. External Subroutines ..
00175       EXTERNAL           SGEMV, SSCAL, SSWAP, XERBLA
00176 *     ..
00177 *     .. Intrinsic Functions ..
00178       INTRINSIC          MAX, SQRT, MAXLOC
00179 *     ..
00180 *     .. Executable Statements ..
00181 *
00182 *     Test the input parameters
00183 *
00184       INFO = 0
00185       UPPER = LSAME( UPLO, 'U' )
00186       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00187          INFO = -1
00188       ELSE IF( N.LT.0 ) THEN
00189          INFO = -2
00190       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00191          INFO = -4
00192       END IF
00193       IF( INFO.NE.0 ) THEN
00194          CALL XERBLA( 'SPSTF2', -INFO )
00195          RETURN
00196       END IF
00197 *
00198 *     Quick return if possible
00199 *
00200       IF( N.EQ.0 )
00201      $   RETURN
00202 *
00203 *     Initialize PIV
00204 *
00205       DO 100 I = 1, N
00206          PIV( I ) = I
00207   100 CONTINUE
00208 *
00209 *     Compute stopping value
00210 *
00211       PVT = 1
00212       AJJ = A( PVT, PVT )
00213       DO I = 2, N
00214          IF( A( I, I ).GT.AJJ ) THEN
00215             PVT = I
00216             AJJ = A( PVT, PVT )
00217          END IF
00218       END DO
00219       IF( AJJ.EQ.ZERO.OR.SISNAN( AJJ ) ) THEN
00220          RANK = 0
00221          INFO = 1
00222          GO TO 170
00223       END IF
00224 *
00225 *     Compute stopping value if not supplied
00226 *
00227       IF( TOL.LT.ZERO ) THEN
00228          SSTOP = N * SLAMCH( 'Epsilon' ) * AJJ
00229       ELSE
00230          SSTOP = TOL
00231       END IF
00232 *
00233 *     Set first half of WORK to zero, holds dot products
00234 *
00235       DO 110 I = 1, N
00236          WORK( I ) = 0
00237   110 CONTINUE
00238 *
00239       IF( UPPER ) THEN
00240 *
00241 *        Compute the Cholesky factorization P**T * A * P = U**T * U
00242 *
00243          DO 130 J = 1, N
00244 *
00245 *        Find pivot, test for exit, else swap rows and columns
00246 *        Update dot products, compute possible pivots which are
00247 *        stored in the second half of WORK
00248 *
00249             DO 120 I = J, N
00250 *
00251                IF( J.GT.1 ) THEN
00252                   WORK( I ) = WORK( I ) + A( J-1, I )**2
00253                END IF
00254                WORK( N+I ) = A( I, I ) - WORK( I )
00255 *
00256   120       CONTINUE
00257 *
00258             IF( J.GT.1 ) THEN
00259                ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
00260                PVT = ITEMP + J - 1
00261                AJJ = WORK( N+PVT )
00262                IF( AJJ.LE.SSTOP.OR.SISNAN( AJJ ) ) THEN
00263                   A( J, J ) = AJJ
00264                   GO TO 160
00265                END IF
00266             END IF
00267 *
00268             IF( J.NE.PVT ) THEN
00269 *
00270 *              Pivot OK, so can now swap pivot rows and columns
00271 *
00272                A( PVT, PVT ) = A( J, J )
00273                CALL SSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
00274                IF( PVT.LT.N )
00275      $            CALL SSWAP( N-PVT, A( J, PVT+1 ), LDA,
00276      $                        A( PVT, PVT+1 ), LDA )
00277                CALL SSWAP( PVT-J-1, A( J, J+1 ), LDA, A( J+1, PVT ), 1 )
00278 *
00279 *              Swap dot products and PIV
00280 *
00281                STEMP = WORK( J )
00282                WORK( J ) = WORK( PVT )
00283                WORK( PVT ) = STEMP
00284                ITEMP = PIV( PVT )
00285                PIV( PVT ) = PIV( J )
00286                PIV( J ) = ITEMP
00287             END IF
00288 *
00289             AJJ = SQRT( AJJ )
00290             A( J, J ) = AJJ
00291 *
00292 *           Compute elements J+1:N of row J
00293 *
00294             IF( J.LT.N ) THEN
00295                CALL SGEMV( 'Trans', J-1, N-J, -ONE, A( 1, J+1 ), LDA,
00296      $                     A( 1, J ), 1, ONE, A( J, J+1 ), LDA )
00297                CALL SSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
00298             END IF
00299 *
00300   130    CONTINUE
00301 *
00302       ELSE
00303 *
00304 *        Compute the Cholesky factorization P**T * A * P = L * L**T
00305 *
00306          DO 150 J = 1, N
00307 *
00308 *        Find pivot, test for exit, else swap rows and columns
00309 *        Update dot products, compute possible pivots which are
00310 *        stored in the second half of WORK
00311 *
00312             DO 140 I = J, N
00313 *
00314                IF( J.GT.1 ) THEN
00315                   WORK( I ) = WORK( I ) + A( I, J-1 )**2
00316                END IF
00317                WORK( N+I ) = A( I, I ) - WORK( I )
00318 *
00319   140       CONTINUE
00320 *
00321             IF( J.GT.1 ) THEN
00322                ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
00323                PVT = ITEMP + J - 1
00324                AJJ = WORK( N+PVT )
00325                IF( AJJ.LE.SSTOP.OR.SISNAN( AJJ ) ) THEN
00326                   A( J, J ) = AJJ
00327                   GO TO 160
00328                END IF
00329             END IF
00330 *
00331             IF( J.NE.PVT ) THEN
00332 *
00333 *              Pivot OK, so can now swap pivot rows and columns
00334 *
00335                A( PVT, PVT ) = A( J, J )
00336                CALL SSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
00337                IF( PVT.LT.N )
00338      $            CALL SSWAP( N-PVT, A( PVT+1, J ), 1, A( PVT+1, PVT ),
00339      $                        1 )
00340                CALL SSWAP( PVT-J-1, A( J+1, J ), 1, A( PVT, J+1 ), LDA )
00341 *
00342 *              Swap dot products and PIV
00343 *
00344                STEMP = WORK( J )
00345                WORK( J ) = WORK( PVT )
00346                WORK( PVT ) = STEMP
00347                ITEMP = PIV( PVT )
00348                PIV( PVT ) = PIV( J )
00349                PIV( J ) = ITEMP
00350             END IF
00351 *
00352             AJJ = SQRT( AJJ )
00353             A( J, J ) = AJJ
00354 *
00355 *           Compute elements J+1:N of column J
00356 *
00357             IF( J.LT.N ) THEN
00358                CALL SGEMV( 'No Trans', N-J, J-1, -ONE, A( J+1, 1 ), LDA,
00359      $                     A( J, 1 ), LDA, ONE, A( J+1, J ), 1 )
00360                CALL SSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
00361             END IF
00362 *
00363   150    CONTINUE
00364 *
00365       END IF
00366 *
00367 *     Ran to completion, A has full rank
00368 *
00369       RANK = N
00370 *
00371       GO TO 170
00372   160 CONTINUE
00373 *
00374 *     Rank is number of steps completed.  Set INFO = 1 to signal
00375 *     that the factorization cannot be used to solve a system.
00376 *
00377       RANK = J - 1
00378       INFO = 1
00379 *
00380   170 CONTINUE
00381       RETURN
00382 *
00383 *     End of SPSTF2
00384 *
00385       END
 All Files Functions