SRC/spotrf2.f

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