SRC/cpotf2.f

   1 *> \brief \b CPOTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite matrix (unblocked algorithm).
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CPOTF2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cpotf2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cpotf2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cpotf2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CPOTF2( UPLO, N, A, LDA, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          UPLO
  25 *       INTEGER            INFO, LDA, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       COMPLEX            A( LDA, * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> CPOTF2 computes the Cholesky factorization of a complex Hermitian
  38 *> positive definite matrix A.
  39 *>
  40 *> The factorization has the form
  41 *>    A = U**H * U ,  if UPLO = 'U', or
  42 *>    A = L  * L**H,  if UPLO = 'L',
  43 *> where U is an upper triangular matrix and L is lower triangular.
  44 *>
  45 *> This is the unblocked version of the algorithm, calling Level 2 BLAS.
  46 *> \endverbatim
  47 *
  48 *  Arguments:
  49 *  ==========
  50 *
  51 *> \param[in] UPLO
  52 *> \verbatim
  53 *>          UPLO is CHARACTER*1
  54 *>          Specifies whether the upper or lower triangular part of the
  55 *>          Hermitian matrix A is stored.
  56 *>          = 'U':  Upper triangular
  57 *>          = 'L':  Lower triangular
  58 *> \endverbatim
  59 *>
  60 *> \param[in] N
  61 *> \verbatim
  62 *>          N is INTEGER
  63 *>          The order of the matrix A.  N >= 0.
  64 *> \endverbatim
  65 *>
  66 *> \param[in,out] A
  67 *> \verbatim
  68 *>          A is COMPLEX array, dimension (LDA,N)
  69 *>          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
  70 *>          n by n upper triangular part of A contains the upper
  71 *>          triangular part of the matrix A, and the strictly lower
  72 *>          triangular part of A is not referenced.  If UPLO = 'L', the
  73 *>          leading n by n lower triangular part of A contains the lower
  74 *>          triangular part of the matrix A, and the strictly upper
  75 *>          triangular part of A is not referenced.
  76 *>
  77 *>          On exit, if INFO = 0, the factor U or L from the Cholesky
  78 *>          factorization A = U**H *U  or A = L*L**H.
  79 *> \endverbatim
  80 *>
  81 *> \param[in] LDA
  82 *> \verbatim
  83 *>          LDA is INTEGER
  84 *>          The leading dimension of the array A.  LDA >= max(1,N).
  85 *> \endverbatim
  86 *>
  87 *> \param[out] INFO
  88 *> \verbatim
  89 *>          INFO is INTEGER
  90 *>          = 0: successful exit
  91 *>          < 0: if INFO = -k, the k-th argument had an illegal value
  92 *>          > 0: if INFO = k, the leading minor of order k is not
  93 *>               positive definite, and the factorization could not be
  94 *>               completed.
  95 *> \endverbatim
  96 *
  97 *  Authors:
  98 *  ========
  99 *
 100 *> \author Univ. of Tennessee
 101 *> \author Univ. of California Berkeley
 102 *> \author Univ. of Colorado Denver
 103 *> \author NAG Ltd.
 104 *
 105 *> \date September 2012
 106 *
 107 *> \ingroup complexPOcomputational
 108 *
 109 *  =====================================================================
 110       SUBROUTINE CPOTF2( UPLO, N, A, LDA, INFO )
 111 *
 112 *  -- LAPACK computational routine (version 3.4.2) --
 113 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 114 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 115 *     September 2012
 116 *
 117 *     .. Scalar Arguments ..
 118       CHARACTER          UPLO
 119       INTEGER            INFO, LDA, N
 120 *     ..
 121 *     .. Array Arguments ..
 122       COMPLEX            A( LDA, * )
 123 *     ..
 124 *
 125 *  =====================================================================
 126 *
 127 *     .. Parameters ..
 128       REAL               ONE, ZERO
 129       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
 130       COMPLEX            CONE
 131       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
 132 *     ..
 133 *     .. Local Scalars ..
 134       LOGICAL            UPPER
 135       INTEGER            J
 136       REAL               AJJ
 137 *     ..
 138 *     .. External Functions ..
 139       LOGICAL            LSAME, SISNAN
 140       COMPLEX            CDOTC
 141       EXTERNAL           LSAME, CDOTC, SISNAN
 142 *     ..
 143 *     .. External Subroutines ..
 144       EXTERNAL           CGEMV, CLACGV, CSSCAL, XERBLA
 145 *     ..
 146 *     .. Intrinsic Functions ..
 147       INTRINSIC          MAX, REAL, SQRT
 148 *     ..
 149 *     .. Executable Statements ..
 150 *
 151 *     Test the input parameters.
 152 *
 153       INFO = 0
 154       UPPER = LSAME( UPLO, 'U' )
 155       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 156          INFO = -1
 157       ELSE IF( N.LT.0 ) THEN
 158          INFO = -2
 159       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 160          INFO = -4
 161       END IF
 162       IF( INFO.NE.0 ) THEN
 163          CALL XERBLA( 'CPOTF2', -INFO )
 164          RETURN
 165       END IF
 166 *
 167 *     Quick return if possible
 168 *
 169       IF( N.EQ.0 )
 170      $   RETURN
 171 *
 172       IF( UPPER ) THEN
 173 *
 174 *        Compute the Cholesky factorization A = U**H *U.
 175 *
 176          DO 10 J = 1, N
 177 *
 178 *           Compute U(J,J) and test for non-positive-definiteness.
 179 *
 180             AJJ = REAL( A( J, J ) ) - CDOTC( J-1, A( 1, J ), 1,
 181      $            A( 1, J ), 1 )
 182             IF( AJJ.LE.ZERO.OR.SISNAN( AJJ ) ) THEN
 183                A( J, J ) = AJJ
 184                GO TO 30
 185             END IF
 186             AJJ = SQRT( AJJ )
 187             A( J, J ) = AJJ
 188 *
 189 *           Compute elements J+1:N of row J.
 190 *
 191             IF( J.LT.N ) THEN
 192                CALL CLACGV( J-1, A( 1, J ), 1 )
 193                CALL CGEMV( 'Transpose', J-1, N-J, -CONE, A( 1, J+1 ),
 194      $                     LDA, A( 1, J ), 1, CONE, A( J, J+1 ), LDA )
 195                CALL CLACGV( J-1, A( 1, J ), 1 )
 196                CALL CSSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
 197             END IF
 198    10    CONTINUE
 199       ELSE
 200 *
 201 *        Compute the Cholesky factorization A = L*L**H.
 202 *
 203          DO 20 J = 1, N
 204 *
 205 *           Compute L(J,J) and test for non-positive-definiteness.
 206 *
 207             AJJ = REAL( A( J, J ) ) - CDOTC( J-1, A( J, 1 ), LDA,
 208      $            A( J, 1 ), LDA )
 209             IF( AJJ.LE.ZERO.OR.SISNAN( AJJ ) ) THEN
 210                A( J, J ) = AJJ
 211                GO TO 30
 212             END IF
 213             AJJ = SQRT( AJJ )
 214             A( J, J ) = AJJ
 215 *
 216 *           Compute elements J+1:N of column J.
 217 *
 218             IF( J.LT.N ) THEN
 219                CALL CLACGV( J-1, A( J, 1 ), LDA )
 220                CALL CGEMV( 'No transpose', N-J, J-1, -CONE, A( J+1, 1 ),
 221      $                     LDA, A( J, 1 ), LDA, CONE, A( J+1, J ), 1 )
 222                CALL CLACGV( J-1, A( J, 1 ), LDA )
 223                CALL CSSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
 224             END IF
 225    20    CONTINUE
 226       END IF
 227       GO TO 40
 228 *
 229    30 CONTINUE
 230       INFO = J
 231 *
 232    40 CONTINUE
 233       RETURN
 234 *
 235 *     End of CPOTF2
 236 *
 237       END