SRC/dpotf2.f

   1 *> \brief \b DPOTF2 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 DPOTF2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dpotf2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dpotf2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dpotf2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DPOTF2( UPLO, N, A, LDA, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          UPLO
  25 *       INTEGER            INFO, LDA, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       DOUBLE PRECISION   A( LDA, * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> DPOTF2 computes the Cholesky factorization of a real symmetric
  38 *> positive definite matrix A.
  39 *>
  40 *> The factorization has the form
  41 *>    A = U**T * U ,  if UPLO = 'U', or
  42 *>    A = L  * L**T,  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 *>          symmetric 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 DOUBLE PRECISION array, dimension (LDA,N)
  69 *>          On entry, the symmetric 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**T *U  or A = L*L**T.
  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 doublePOcomputational
 108 *
 109 *  =====================================================================
 110       SUBROUTINE DPOTF2( 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       DOUBLE PRECISION   A( LDA, * )
 123 *     ..
 124 *
 125 *  =====================================================================
 126 *
 127 *     .. Parameters ..
 128       DOUBLE PRECISION   ONE, ZERO
 129       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 130 *     ..
 131 *     .. Local Scalars ..
 132       LOGICAL            UPPER
 133       INTEGER            J
 134       DOUBLE PRECISION   AJJ
 135 *     ..
 136 *     .. External Functions ..
 137       LOGICAL            LSAME, DISNAN
 138       DOUBLE PRECISION   DDOT
 139       EXTERNAL           LSAME, DDOT, DISNAN
 140 *     ..
 141 *     .. External Subroutines ..
 142       EXTERNAL           DGEMV, DSCAL, XERBLA
 143 *     ..
 144 *     .. Intrinsic Functions ..
 145       INTRINSIC          MAX, SQRT
 146 *     ..
 147 *     .. Executable Statements ..
 148 *
 149 *     Test the input parameters.
 150 *
 151       INFO = 0
 152       UPPER = LSAME( UPLO, 'U' )
 153       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 154          INFO = -1
 155       ELSE IF( N.LT.0 ) THEN
 156          INFO = -2
 157       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 158          INFO = -4
 159       END IF
 160       IF( INFO.NE.0 ) THEN
 161          CALL XERBLA( 'DPOTF2', -INFO )
 162          RETURN
 163       END IF
 164 *
 165 *     Quick return if possible
 166 *
 167       IF( N.EQ.0 )
 168      $   RETURN
 169 *
 170       IF( UPPER ) THEN
 171 *
 172 *        Compute the Cholesky factorization A = U**T *U.
 173 *
 174          DO 10 J = 1, N
 175 *
 176 *           Compute U(J,J) and test for non-positive-definiteness.
 177 *
 178             AJJ = A( J, J ) - DDOT( J-1, A( 1, J ), 1, A( 1, J ), 1 )
 179             IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN
 180                A( J, J ) = AJJ
 181                GO TO 30
 182             END IF
 183             AJJ = SQRT( AJJ )
 184             A( J, J ) = AJJ
 185 *
 186 *           Compute elements J+1:N of row J.
 187 *
 188             IF( J.LT.N ) THEN
 189                CALL DGEMV( 'Transpose', J-1, N-J, -ONE, A( 1, J+1 ),
 190      $                     LDA, A( 1, J ), 1, ONE, A( J, J+1 ), LDA )
 191                CALL DSCAL( N-J, ONE / AJJ, A( J, J+1 ), LDA )
 192             END IF
 193    10    CONTINUE
 194       ELSE
 195 *
 196 *        Compute the Cholesky factorization A = L*L**T.
 197 *
 198          DO 20 J = 1, N
 199 *
 200 *           Compute L(J,J) and test for non-positive-definiteness.
 201 *
 202             AJJ = A( J, J ) - DDOT( J-1, A( J, 1 ), LDA, A( J, 1 ),
 203      $            LDA )
 204             IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN
 205                A( J, J ) = AJJ
 206                GO TO 30
 207             END IF
 208             AJJ = SQRT( AJJ )
 209             A( J, J ) = AJJ
 210 *
 211 *           Compute elements J+1:N of column J.
 212 *
 213             IF( J.LT.N ) THEN
 214                CALL DGEMV( 'No transpose', N-J, J-1, -ONE, A( J+1, 1 ),
 215      $                     LDA, A( J, 1 ), LDA, ONE, A( J+1, J ), 1 )
 216                CALL DSCAL( N-J, ONE / AJJ, A( J+1, J ), 1 )
 217             END IF
 218    20    CONTINUE
 219       END IF
 220       GO TO 40
 221 *
 222    30 CONTINUE
 223       INFO = J
 224 *
 225    40 CONTINUE
 226       RETURN
 227 *
 228 *     End of DPOTF2
 229 *
 230       END