LAPACK 3.3.1
Linear Algebra PACKage

cpstrf.f

Go to the documentation of this file.
00001       SUBROUTINE CPSTRF( UPLO, N, A, LDA, PIV, RANK, TOL, WORK, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2.2) --
00004 *     Craig Lucas, University of Manchester / NAG Ltd.
00005 *     October, 2008
00006 *
00007 *     .. Scalar Arguments ..
00008       REAL               TOL
00009       INTEGER            INFO, LDA, N, RANK
00010       CHARACTER          UPLO
00011 *     ..
00012 *     .. Array Arguments ..
00013       COMPLEX            A( LDA, * )
00014       REAL               WORK( 2*N )
00015       INTEGER            PIV( N )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  CPSTRF 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**T * A * P = U**H * U ,  if UPLO = 'U',
00026 *     P**T * A * P = L  * L**H,  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 3 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 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 *  LDA     (input) INTEGER
00058 *          The leading dimension of the array A.  LDA >= max(1,N).
00059 *
00060 *  PIV     (output) INTEGER array, dimension (N)
00061 *          PIV is such that the nonzero entries are P( PIV(K), K ) = 1.
00062 *
00063 *  RANK    (output) INTEGER
00064 *          The rank of A given by the number of steps the algorithm
00065 *          completed.
00066 *
00067 *  TOL     (input) REAL
00068 *          User defined tolerance. If TOL < 0, then N*U*MAX( A(K,K) )
00069 *          will be used. The algorithm terminates at the (K-1)st step
00070 *          if the pivot <= TOL.
00071 *
00072 *  WORK    (workspace) REAL 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       REAL               ONE, ZERO
00086       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00087       COMPLEX            CONE
00088       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
00089 *     ..
00090 *     .. Local Scalars ..
00091       COMPLEX            CTEMP
00092       REAL               AJJ, SSTOP, STEMP
00093       INTEGER            I, ITEMP, J, JB, K, NB, PVT
00094       LOGICAL            UPPER
00095 *     ..
00096 *     .. External Functions ..
00097       REAL               SLAMCH
00098       INTEGER            ILAENV
00099       LOGICAL            LSAME, SISNAN
00100       EXTERNAL           SLAMCH, ILAENV, LSAME, SISNAN
00101 *     ..
00102 *     .. External Subroutines ..
00103       EXTERNAL           CGEMV, CHERK, CLACGV, CPSTF2, CSSCAL, CSWAP,
00104      $                   XERBLA
00105 *     ..
00106 *     .. Intrinsic Functions ..
00107       INTRINSIC          CONJG, MAX, MIN, REAL, SQRT, MAXLOC
00108 *     ..
00109 *     .. Executable Statements ..
00110 *
00111 *     Test the input parameters.
00112 *
00113       INFO = 0
00114       UPPER = LSAME( UPLO, 'U' )
00115       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00116          INFO = -1
00117       ELSE IF( N.LT.0 ) THEN
00118          INFO = -2
00119       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00120          INFO = -4
00121       END IF
00122       IF( INFO.NE.0 ) THEN
00123          CALL XERBLA( 'CPSTRF', -INFO )
00124          RETURN
00125       END IF
00126 *
00127 *     Quick return if possible
00128 *
00129       IF( N.EQ.0 )
00130      $   RETURN
00131 *
00132 *     Get block size
00133 *
00134       NB = ILAENV( 1, 'CPOTRF', UPLO, N, -1, -1, -1 )
00135       IF( NB.LE.1 .OR. NB.GE.N ) THEN
00136 *
00137 *        Use unblocked code
00138 *
00139          CALL CPSTF2( UPLO, N, A( 1, 1 ), LDA, PIV, RANK, TOL, WORK,
00140      $                INFO )
00141          GO TO 230
00142 *
00143       ELSE
00144 *
00145 *     Initialize PIV
00146 *
00147          DO 100 I = 1, N
00148             PIV( I ) = I
00149   100    CONTINUE
00150 *
00151 *     Compute stopping value
00152 *
00153          DO 110 I = 1, N
00154             WORK( I ) = REAL( A( I, I ) )
00155   110    CONTINUE
00156          PVT = MAXLOC( WORK( 1:N ), 1 )
00157          AJJ = REAL( A( PVT, PVT ) )
00158          IF( AJJ.EQ.ZERO.OR.SISNAN( AJJ ) ) THEN
00159             RANK = 0
00160             INFO = 1
00161             GO TO 230
00162          END IF
00163 *
00164 *     Compute stopping value if not supplied
00165 *
00166          IF( TOL.LT.ZERO ) THEN
00167             SSTOP = N * SLAMCH( 'Epsilon' ) * AJJ
00168          ELSE
00169             SSTOP = TOL
00170          END IF
00171 *
00172 *
00173          IF( UPPER ) THEN
00174 *
00175 *           Compute the Cholesky factorization P**T * A * P = U**H * U
00176 *
00177             DO 160 K = 1, N, NB
00178 *
00179 *              Account for last block not being NB wide
00180 *
00181                JB = MIN( NB, N-K+1 )
00182 *
00183 *              Set relevant part of first half of WORK to zero,
00184 *              holds dot products
00185 *
00186                DO 120 I = K, N
00187                   WORK( I ) = 0
00188   120          CONTINUE
00189 *
00190                DO 150 J = K, K + JB - 1
00191 *
00192 *              Find pivot, test for exit, else swap rows and columns
00193 *              Update dot products, compute possible pivots which are
00194 *              stored in the second half of WORK
00195 *
00196                   DO 130 I = J, N
00197 *
00198                      IF( J.GT.K ) THEN
00199                         WORK( I ) = WORK( I ) +
00200      $                              REAL( CONJG( A( J-1, I ) )*
00201      $                                    A( J-1, I ) )
00202                      END IF
00203                      WORK( N+I ) = REAL( A( I, I ) ) - WORK( I )
00204 *
00205   130             CONTINUE
00206 *
00207                   IF( J.GT.1 ) THEN
00208                      ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
00209                      PVT = ITEMP + J - 1
00210                      AJJ = WORK( N+PVT )
00211                      IF( AJJ.LE.SSTOP.OR.SISNAN( AJJ ) ) THEN
00212                         A( J, J ) = AJJ
00213                         GO TO 220
00214                      END IF
00215                   END IF
00216 *
00217                   IF( J.NE.PVT ) THEN
00218 *
00219 *                    Pivot OK, so can now swap pivot rows and columns
00220 *
00221                      A( PVT, PVT ) = A( J, J )
00222                      CALL CSWAP( J-1, A( 1, J ), 1, A( 1, PVT ), 1 )
00223                      IF( PVT.LT.N )
00224      $                  CALL CSWAP( N-PVT, A( J, PVT+1 ), LDA,
00225      $                              A( PVT, PVT+1 ), LDA )
00226                      DO 140 I = J + 1, PVT - 1
00227                         CTEMP = CONJG( A( J, I ) )
00228                         A( J, I ) = CONJG( A( I, PVT ) )
00229                         A( I, PVT ) = CTEMP
00230   140                CONTINUE
00231                      A( J, PVT ) = CONJG( A( J, PVT ) )
00232 *
00233 *                    Swap dot products and PIV
00234 *
00235                      STEMP = WORK( J )
00236                      WORK( J ) = WORK( PVT )
00237                      WORK( PVT ) = STEMP
00238                      ITEMP = PIV( PVT )
00239                      PIV( PVT ) = PIV( J )
00240                      PIV( J ) = ITEMP
00241                   END IF
00242 *
00243                   AJJ = SQRT( AJJ )
00244                   A( J, J ) = AJJ
00245 *
00246 *                 Compute elements J+1:N of row J.
00247 *
00248                   IF( J.LT.N ) THEN
00249                      CALL CLACGV( J-1, A( 1, J ), 1 )
00250                      CALL CGEMV( 'Trans', J-K, N-J, -CONE, A( K, J+1 ),
00251      $                           LDA, A( K, J ), 1, CONE, A( J, J+1 ),
00252      $                           LDA )
00253                      CALL CLACGV( J-1, A( 1, J ), 1 )
00254                      CALL CSSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
00255                   END IF
00256 *
00257   150          CONTINUE
00258 *
00259 *              Update trailing matrix, J already incremented
00260 *
00261                IF( K+JB.LE.N ) THEN
00262                   CALL CHERK( 'Upper', 'Conj Trans', N-J+1, JB, -ONE,
00263      $                        A( K, J ), LDA, ONE, A( J, J ), LDA )
00264                END IF
00265 *
00266   160       CONTINUE
00267 *
00268          ELSE
00269 *
00270 *        Compute the Cholesky factorization P**T * A * P = L * L**H
00271 *
00272             DO 210 K = 1, N, NB
00273 *
00274 *              Account for last block not being NB wide
00275 *
00276                JB = MIN( NB, N-K+1 )
00277 *
00278 *              Set relevant part of first half of WORK to zero,
00279 *              holds dot products
00280 *
00281                DO 170 I = K, N
00282                   WORK( I ) = 0
00283   170          CONTINUE
00284 *
00285                DO 200 J = K, K + JB - 1
00286 *
00287 *              Find pivot, test for exit, else swap rows and columns
00288 *              Update dot products, compute possible pivots which are
00289 *              stored in the second half of WORK
00290 *
00291                   DO 180 I = J, N
00292 *
00293                      IF( J.GT.K ) THEN
00294                         WORK( I ) = WORK( I ) +
00295      $                              REAL( CONJG( A( I, J-1 ) )*
00296      $                                    A( I, J-1 ) )
00297                      END IF
00298                      WORK( N+I ) = REAL( A( I, I ) ) - WORK( I )
00299 *
00300   180             CONTINUE
00301 *
00302                   IF( J.GT.1 ) THEN
00303                      ITEMP = MAXLOC( WORK( (N+J):(2*N) ), 1 )
00304                      PVT = ITEMP + J - 1
00305                      AJJ = WORK( N+PVT )
00306                      IF( AJJ.LE.SSTOP.OR.SISNAN( AJJ ) ) THEN
00307                         A( J, J ) = AJJ
00308                         GO TO 220
00309                      END IF
00310                   END IF
00311 *
00312                   IF( J.NE.PVT ) THEN
00313 *
00314 *                    Pivot OK, so can now swap pivot rows and columns
00315 *
00316                      A( PVT, PVT ) = A( J, J )
00317                      CALL CSWAP( J-1, A( J, 1 ), LDA, A( PVT, 1 ), LDA )
00318                      IF( PVT.LT.N )
00319      $                  CALL CSWAP( N-PVT, A( PVT+1, J ), 1,
00320      $                              A( PVT+1, PVT ), 1 )
00321                      DO 190 I = J + 1, PVT - 1
00322                         CTEMP = CONJG( A( I, J ) )
00323                         A( I, J ) = CONJG( A( PVT, I ) )
00324                         A( PVT, I ) = CTEMP
00325   190                CONTINUE
00326                      A( PVT, J ) = CONJG( A( PVT, J ) )
00327 *
00328 *                    Swap dot products and PIV
00329 *
00330                      STEMP = WORK( J )
00331                      WORK( J ) = WORK( PVT )
00332                      WORK( PVT ) = STEMP
00333                      ITEMP = PIV( PVT )
00334                      PIV( PVT ) = PIV( J )
00335                      PIV( J ) = ITEMP
00336                   END IF
00337 *
00338                   AJJ = SQRT( AJJ )
00339                   A( J, J ) = AJJ
00340 *
00341 *                 Compute elements J+1:N of column J.
00342 *
00343                   IF( J.LT.N ) THEN
00344                      CALL CLACGV( J-1, A( J, 1 ), LDA )
00345                      CALL CGEMV( 'No Trans', N-J, J-K, -CONE,
00346      $                           A( J+1, K ), LDA, A( J, K ), LDA, CONE,
00347      $                           A( J+1, J ), 1 )
00348                      CALL CLACGV( J-1, A( J, 1 ), LDA )
00349                      CALL CSSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
00350                   END IF
00351 *
00352   200          CONTINUE
00353 *
00354 *              Update trailing matrix, J already incremented
00355 *
00356                IF( K+JB.LE.N ) THEN
00357                   CALL CHERK( 'Lower', 'No Trans', N-J+1, JB, -ONE,
00358      $                        A( J, K ), LDA, ONE, A( J, J ), LDA )
00359                END IF
00360 *
00361   210       CONTINUE
00362 *
00363          END IF
00364       END IF
00365 *
00366 *     Ran to completion, A has full rank
00367 *
00368       RANK = N
00369 *
00370       GO TO 230
00371   220 CONTINUE
00372 *
00373 *     Rank is the number of steps completed.  Set INFO = 1 to signal
00374 *     that the factorization cannot be used to solve a system.
00375 *
00376       RANK = J - 1
00377       INFO = 1
00378 *
00379   230 CONTINUE
00380       RETURN
00381 *
00382 *     End of CPSTRF
00383 *
00384       END
 All Files Functions