LAPACK: SRC/zpstf2.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZPSTF2( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
00002 *
00003 *  -- LAPACK PROTOTYPE routine (version 3.2.2) --
00004 *     Craig Lucas, University of Manchester / NAG Ltd.
00005 *     October, 2008
00006 *
00007 *     .. Scalar Arguments ..
00008       DOUBLE PRECISION   TOL
00009       INTEGER            INFO, LDA, N, RANK
00010       CHARACTER          UPLO
00011 *     ..
00012 *     .. Array Arguments ..
00013       COMPLEX*16         A( LDA, * )
00014       DOUBLE PRECISION   WORK( 2*N )
00015       INTEGER            PIV( N )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  ZPSTF2 computes the Cholesky factorization with complete
00022 *  pivoting of a complex Hermitian positive semidefinite matrix A.
00023 *
00024 *  The factorization has the form
00025 *     P' * A * P = U' * U ,  if UPLO = 'U',
00026 *     P' * A * P = L  * L',  if UPLO = 'L',
00027 *  where U is an upper triangular matrix and L is lower triangular, and
00028 *  P is stored as vector PIV.
00029 *
00030 *  This algorithm does not attempt to check that A is positive
00031 *  semidefinite. This version of the algorithm calls level 2 BLAS.
00032 *
00033 *  Arguments
00034 *  =========
00035 *
00036 *  UPLO    (input) CHARACTER*1
00037 *          Specifies whether the upper or lower triangular part of the
00038 *          symmetric matrix A is stored.
00039 *          = 'U':  Upper triangular
00040 *          = 'L':  Lower triangular
00041 *
00042 *  N       (input) INTEGER
00043 *          The order of the matrix A.  N >= 0.
00044 *
00045 *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
00046 *          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
00047 *          n by n upper triangular part of A contains the upper
00048 *          triangular part of the matrix A, and the strictly lower
00049 *          triangular part of A is not referenced.  If UPLO = 'L', the
00050 *          leading n by n lower triangular part of A contains the lower
00051 *          triangular part of the matrix A, and the strictly upper
00052 *          triangular part of A is not referenced.
00053 *
00054 *          On exit, if INFO = 0, the factor U or L from the Cholesky
00055 *          factorization as above.
00056 *
00057 *  PIV     (output) INTEGER array, dimension (N)
00058 *          PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
00059 *
00060 *  RANK    (output) INTEGER
00061 *          The rank of A given by the number of steps the algorithm
00062 *          completed.
00063 *
00064 *  TOL     (input) DOUBLE PRECISION
00065 *          User defined tolerance. If TOL < 0, then N*U*MAX( A( K,K ) )
00066 *          will be used. The algorithm terminates at the (K-1)st step
00067 *          if the pivot <= TOL.
00068 *
00069 *  LDA     (input) INTEGER
00070 *          The leading dimension of the array A.  LDA >= max(1,N).
00071 *
00072 *  WORK    (workspace) DOUBLE PRECISION array, dimension (2*N)
00073 *          Work space.
00074 *
00075 *  INFO    (output) INTEGER
00076 *          < 0: If INFO = -K, the K-th argument had an illegal value,
00077 *          = 0: algorithm completed successfully, and
00078 *          > 0: the matrix A is either rank deficient with computed rank
00079 *               as returned in RANK, or is indefinite.  See Section 7 of
00080 *               LAPACK Working Note #161 for further information.
00081 *
00082 *  =====================================================================
00083 *
00084 *     .. Parameters ..
00085       DOUBLE PRECISION   ONE, ZERO
00086       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00087       COMPLEX*16         CONE
00088       PARAMETER          ( CONE = ( 1.0D+0, 0.0D+0 ) )
00089 *     ..
00090 *     .. Local Scalars ..
00091       COMPLEX*16         ZTEMP
00092       DOUBLE PRECISION   AJJ, DSTOP, DTEMP
00093       INTEGER            I, ITEMP, J, PVT
00094       LOGICAL            UPPER
00095 *     ..
00096 *     .. External Functions ..
00097       DOUBLE PRECISION   DLAMCH
00098       LOGICAL            LSAME, DISNAN
00099       EXTERNAL           DLAMCH, LSAME, DISNAN
00100 *     ..
00101 *     .. External Subroutines ..
00102       EXTERNAL           ZDSCAL, ZGEMV, ZLACGV, ZSWAP, XERBLA
00103 *     ..
00104 *     .. Intrinsic Functions ..
00105       INTRINSIC          DBLE, DCONJG, MAX, SQRT
00106 *     ..
00107 *     .. Executable Statements ..
00108 *
00109 *     Test the input parameters
00110 *
00111       INFO = 0
00112       UPPER = LSAME( UPLO, 'U' )
00113       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00114          INFO = -1
00115       ELSE IF( N.LT.0 ) THEN
00116          INFO = -2
00117       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00118          INFO = -4
00119       END IF
00120       IF( INFO.NE.0 ) THEN
00121          CALL XERBLA( 'ZPSTF2', -INFO )
00122          RETURN
00123       END IF
00124 *
00125 *     Quick return if possible
00126 *
00127       IF( N.EQ.0 )
00128      $   RETURN
00129 *
00130 *     Initialize PIV
00131 *
00132       DO 100 I = 1, N
00133          PIV( I ) = I
00134   100 CONTINUE
00135 *
00136 *     Compute stopping value
00137 *
00138       DO 110 I = 1, N
00139          WORK( I ) = DBLE( A( I, I ) )
00140   110 CONTINUE
00141       PVT = MAXLOC( WORK( 1:N ), 1 )
00142       AJJ = DBLE( A( PVT, PVT ) )
00143       IF( AJJ.EQ.ZERO.OR.DISNAN( AJJ ) ) THEN
00144          RANK = 0
00145          INFO = 1
00146          GO TO 200
00147       END IF
00148 *
00149 *     Compute stopping value if not supplied
00150 *
00151       IF( TOL.LT.ZERO ) THEN
00152          DSTOP = N * DLAMCH( 'Epsilon' ) * AJJ
00153       ELSE
00154          DSTOP = TOL
00155       END IF
00156 *
00157 *     Set first half of WORK to zero, holds dot products
00158 *
00159       DO 120 I = 1, N
00160          WORK( I ) = 0
00161   120 CONTINUE
00162 *
00163       IF( UPPER ) THEN
00164 *
00165 *        Compute the Cholesky factorization P' * A * P = U' * U
00166 *
00167          DO 150 J = 1, N
00168 *
00169 *        Find pivot, test for exit, else swap rows and columns
00170 *        Update dot products, compute possible pivots which are
00171 *        stored in the second half of WORK
00172 *
00173             DO 130 I = J, N
00174 *
00175                IF( J.GT.1 ) THEN
00176                   WORK( I ) = WORK( I ) + 
00177      $                        DBLE( DCONJG( A( J-1, I ) )*
00178      $                              A( J-1, I ) )
00179                END IF
00180                WORK( N+I ) = DBLE( A( I, I ) ) - WORK( I )
00181 *
00182   130       CONTINUE
00183 *
00184             IF( J.GT.1 ) THEN
00185                ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
00186                PVT = ITEMP + J - 1
00187                AJJ = WORK( N+PVT )
00188                IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
00189                   A( J, J ) = AJJ
00190                   GO TO 190
00191                END IF
00192             END IF
00193 *
00194             IF( J.NE.PVT ) THEN
00195 *
00196 *              Pivot OK, so can now swap pivot rows and columns
00197 *
00198                A( PVT, PVT ) = A( J, J )
00199                CALL ZSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
00200                IF( PVT.LT.N )
00201      $            CALL ZSWAP( N-PVT, A( J, PVT+1 ), LDA,
00202      $                        A( PVT, PVT+1 ), LDA )
00203                DO 140 I = J + 1, PVT - 1
00204                   ZTEMP = DCONJG( A( J, I ) )
00205                   A( J, I ) = DCONJG( A( I, PVT ) )
00206                   A( I, PVT ) = ZTEMP
00207   140          CONTINUE
00208                A( J, PVT ) = DCONJG( A( J, PVT ) )
00209 *
00210 *              Swap dot products and PIV
00211 *
00212                DTEMP = WORK( J )
00213                WORK( J ) = WORK( PVT )
00214                WORK( PVT ) = DTEMP
00215                ITEMP = PIV( PVT )
00216                PIV( PVT ) = PIV( J )
00217                PIV( J ) = ITEMP
00218             END IF
00219 *
00220             AJJ = SQRT( AJJ )
00221             A( J, J ) = AJJ
00222 *
00223 *           Compute elements J+1:N of row J
00224 *
00225             IF( J.LT.N ) THEN
00226                CALL ZLACGV( J-1, A( 1, J ), 1 )
00227                CALL ZGEMV( 'Trans', J-1, N-J, -CONE, A( 1, J+1 ), LDA,
00228      $                     A( 1, J ), 1, CONE, A( J, J+1 ), LDA )
00229                CALL ZLACGV( J-1, A( 1, J ), 1 )
00230                CALL ZDSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
00231             END IF
00232 *
00233   150    CONTINUE
00234 *
00235       ELSE
00236 *
00237 *        Compute the Cholesky factorization P' * A * P = L * L'
00238 *
00239          DO 180 J = 1, N
00240 *
00241 *        Find pivot, test for exit, else swap rows and columns
00242 *        Update dot products, compute possible pivots which are
00243 *        stored in the second half of WORK
00244 *
00245             DO 160 I = J, N
00246 *
00247                IF( J.GT.1 ) THEN
00248                   WORK( I ) = WORK( I ) + 
00249      $                        DBLE( DCONJG( A( I, J-1 ) )*
00250      $                              A( I, J-1 ) )
00251                END IF
00252                WORK( N+I ) = DBLE( A( I, I ) ) - WORK( I )
00253 *
00254   160       CONTINUE
00255 *
00256             IF( J.GT.1 ) THEN
00257                ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
00258                PVT = ITEMP + J - 1
00259                AJJ = WORK( N+PVT )
00260                IF( AJJ.LE.DSTOP.OR.DISNAN( AJJ ) ) THEN
00261                   A( J, J ) = AJJ
00262                   GO TO 190
00263                END IF
00264             END IF
00265 *
00266             IF( J.NE.PVT ) THEN
00267 *
00268 *              Pivot OK, so can now swap pivot rows and columns
00269 *
00270                A( PVT, PVT ) = A( J, J )
00271                CALL ZSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
00272                IF( PVT.LT.N )
00273      $            CALL ZSWAP( N-PVT, A( PVT+1, J ), 1, A( PVT+1, PVT ),
00274      $                        1 )
00275                DO 170 I = J + 1, PVT - 1
00276                   ZTEMP = DCONJG( A( I, J ) )
00277                   A( I, J ) = DCONJG( A( PVT, I ) )
00278                   A( PVT, I ) = ZTEMP
00279   170          CONTINUE
00280                A( PVT, J ) = DCONJG( A( PVT, J ) )
00281 *
00282 *              Swap dot products and PIV
00283 *
00284                DTEMP = WORK( J )
00285                WORK( J ) = WORK( PVT )
00286                WORK( PVT ) = DTEMP
00287                ITEMP = PIV( PVT )
00288                PIV( PVT ) = PIV( J )
00289                PIV( J ) = ITEMP
00290             END IF
00291 *
00292             AJJ = SQRT( AJJ )
00293             A( J, J ) = AJJ
00294 *
00295 *           Compute elements J+1:N of column J
00296 *
00297             IF( J.LT.N ) THEN
00298                CALL ZLACGV( J-1, A( J, 1 ), LDA )
00299                CALL ZGEMV( 'No Trans', N-J, J-1, -CONE, A( J+1, 1 ),
00300      $                     LDA, A( J, 1 ), LDA, CONE, A( J+1, J ), 1 )
00301                CALL ZLACGV( J-1, A( J, 1 ), LDA )
00302                CALL ZDSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
00303             END IF
00304 *
00305   180    CONTINUE
00306 *
00307       END IF
00308 *
00309 *     Ran to completion, A has full rank
00310 *
00311       RANK = N
00312 *
00313       GO TO 200
00314   190 CONTINUE
00315 *
00316 *     Rank is number of steps completed.  Set INFO = 1 to signal
00317 *     that the factorization cannot be used to solve a system.
00318 *
00319       RANK = J - 1
00320       INFO = 1
00321 *
00322   200 CONTINUE
00323       RETURN
00324 *
00325 *     End of ZPSTF2
00326 *
00327       END