reference/dpotrff.f

   1       SUBROUTINE DPOTRFF( UPLO, N, A, LDA, INFO )
   2 *
   3 *  -- LAPACK routine (version 3.0) --
   4 *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
   5 *     Courant Institute, Argonne National Lab, and Rice University
   6 *     March 31, 1993
   7 *
   8 *     .. Scalar Arguments ..
   9       CHARACTER          UPLO
  10       INTEGER            INFO, LDA, N
  11 *     ..
  12 *     .. Array Arguments ..
  13       DOUBLE PRECISION   A( LDA, * )
  14 *     ..
  15 *
  16 *  Purpose
  17 *  =======
  18 *
  19 *  DPOTRF computes the Cholesky factorization of a real symmetric
  20 *  positive definite matrix A.
  21 *
  22 *  The factorization has the form
  23 *     A = U**T * U,  if UPLO = 'U', or
  24 *     A = L  * L**T,  if UPLO = 'L',
  25 *  where U is an upper triangular matrix and L is lower triangular.
  26 *
  27 *  This is the block version of the algorithm, calling Level 3 BLAS.
  28 *
  29 *  Arguments
  30 *  =========
  31 *
  32 *  UPLO    (input) CHARACTER*1
  33 *          = 'U':  Upper triangle of A is stored;
  34 *          = 'L':  Lower triangle of A is stored.
  35 *
  36 *  N       (input) INTEGER
  37 *          The order of the matrix A.  N >= 0.
  38 *
  39 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  40 *          On entry, the symmetric matrix A.  If UPLO = 'U', the leading
  41 *          N-by-N upper triangular part of A contains the upper
  42 *          triangular part of the matrix A, and the strictly lower
  43 *          triangular part of A is not referenced.  If UPLO = 'L', the
  44 *          leading N-by-N lower triangular part of A contains the lower
  45 *          triangular part of the matrix A, and the strictly upper
  46 *          triangular part of A is not referenced.
  47 *
  48 *          On exit, if INFO = 0, the factor U or L from the Cholesky
  49 *          factorization A = U**T*U or A = L*L**T.
  50 *
  51 *  LDA     (input) INTEGER
  52 *          The leading dimension of the array A.  LDA >= max(1,N).
  53 *
  54 *  INFO    (output) INTEGER
  55 *          = 0:  successful exit
  56 *          < 0:  if INFO = -i, the i-th argument had an illegal value
  57 *          > 0:  if INFO = i, the leading minor of order i is not
  58 *                positive definite, and the factorization could not be
  59 *                completed.
  60 *
  61 *  =====================================================================
  62 *
  63 *     .. Parameters ..
  64       DOUBLE PRECISION   ONE
  65       PARAMETER          ( ONE = 1.0D+0 )
  66 *     ..
  67 *     .. Local Scalars ..
  68       LOGICAL            UPPER
  69       INTEGER            J, JB, NB
  70 *     ..
  71 *     .. External Functions ..
  72       LOGICAL            LSAME
  73       EXTERNAL           LSAME
  74 *     ..
  75 *     .. External Subroutines ..
  76       EXTERNAL           DGEMM, DPOTF2, DSYRK, DTRSM, XERBLA
  77 *     ..
  78 *     .. Intrinsic Functions ..
  79       INTRINSIC          MAX, MIN
  80 *     ..
  81 *     .. Executable Statements ..
  82 *
  83 *     Test the input parameters.
  84 *
  85       INFO = 0
  86       UPPER = LSAME( UPLO, 'U' )
  87       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
  88          INFO = -1
  89       ELSE IF( N.LT.0 ) THEN
  90          INFO = -2
  91       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
  92          INFO = -4
  93       END IF
  94       IF( INFO.NE.0 ) THEN
  95          CALL XERBLA( 'DPOTRF', -INFO )
  96          RETURN
  97       END IF
  98 *
  99 *     Quick return if possible
 100 *
 101       IF( N.EQ.0 )
 102      $   RETURN
 103 *
 104 *     Determine the block size for this environment.
 105 *
 106       NB = 224
 107
 108       IF( NB.LE.1 .OR. NB.GE.N ) THEN
 109 *
 110 *        Use unblocked code.
 111 *
 112          CALL DPOTF2( UPLO, N, A, LDA, INFO )
 113       ELSE
 114 *
 115 *        Use blocked code.
 116 *
 117          IF( UPPER ) THEN
 118 *
 119 *           Compute the Cholesky factorization A = U'*U.
 120 *
 121             DO 10 J = 1, N, NB
 122 *
 123 *              Update and factorize the current diagonal block and test
 124 *              for non-positive-definiteness.
 125 *
 126                JB = MIN( NB, N-J+1 )
 127                CALL DSYRK( 'Upper', 'Transpose', JB, J-1, -ONE,
 128      $                     A( 1, J ), LDA, ONE, A( J, J ), LDA )
 129                CALL DPOTF2( 'Upper', JB, A( J, J ), LDA, INFO )
 130                IF( INFO.NE.0 )
 131      $            GO TO 30
 132                IF( J+JB.LE.N ) THEN
 133 *
 134 *                 Compute the current block row.
 135 *
 136                   CALL DGEMM( 'Transpose', 'No transpose', JB, N-J-JB+1,
 137      $                        J-1, -ONE, A( 1, J ), LDA, A( 1, J+JB ),
 138      $                        LDA, ONE, A( J, J+JB ), LDA )
 139                   CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit',
 140      $                        JB, N-J-JB+1, ONE, A( J, J ), LDA,
 141      $                        A( J, J+JB ), LDA )
 142                END IF
 143    10       CONTINUE
 144 *
 145          ELSE
 146 *
 147 *           Compute the Cholesky factorization A = L*L'.
 148 *
 149             DO 20 J = 1, N, NB
 150 *
 151 *              Update and factorize the current diagonal block and test
 152 *              for non-positive-definiteness.
 153 *
 154                JB = MIN( NB, N-J+1 )
 155                CALL DSYRK( 'Lower', 'No transpose', JB, J-1, -ONE,
 156      $                     A( J, 1 ), LDA, ONE, A( J, J ), LDA )
 157                CALL DPOTF2( 'Lower', JB, A( J, J ), LDA, INFO )
 158                IF( INFO.NE.0 )
 159      $            GO TO 30
 160                IF( J+JB.LE.N ) THEN
 161 *
 162 *                 Compute the current block column.
 163 *
 164                   CALL DGEMM( 'No transpose', 'Transpose', N-J-JB+1, JB,
 165      $                        J-1, -ONE, A( J+JB, 1 ), LDA, A( J, 1 ),
 166      $                        LDA, ONE, A( J+JB, J ), LDA )
 167                   CALL DTRSM( 'Right', 'Lower', 'Transpose', 'Non-unit',
 168      $                        N-J-JB+1, JB, ONE, A( J, J ), LDA,
 169      $                        A( J+JB, J ), LDA )
 170                END IF
 171    20       CONTINUE
 172          END IF
 173       END IF
 174       GO TO 40
 175 *
 176    30 CONTINUE
 177       INFO = INFO + J - 1
 178 *
 179    40 CONTINUE
 180       RETURN
 181 *
 182 *     End of DPOTRF
 183 *
 184       END