SRC/zpotrf2.f

   1 *> \brief \b ZPOTRF2
   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 ZPOTRF2( UPLO, N, A, LDA, INFO )
  12 *
  13 *       .. Scalar Arguments ..
  14 *       CHARACTER          UPLO
  15 *       INTEGER            INFO, LDA, N
  16 *       ..
  17 *       .. Array Arguments ..
  18 *       COMPLEX*16         A( LDA, * )
  19 *       ..
  20 *
  21 *
  22 *> \par Purpose:
  23 *  =============
  24 *>
  25 *> \verbatim
  26 *>
  27 *> ZPOTRF2 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**H * U,  if UPLO = 'U', or
  32 *>    A = L  * L**H,  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 COMPLEX*16 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**H*U or A = L*L**H.
  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 complex16POcomputational
 105 *
 106 *  =====================================================================
 107       RECURSIVE SUBROUTINE ZPOTRF2( 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       COMPLEX*16         A( LDA, * )
 120 *     ..
 121 *
 122 *  =====================================================================
 123 *
 124 *     .. Parameters ..
 125       DOUBLE PRECISION   ONE, ZERO
 126       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 127       COMPLEX*16         CONE
 128       PARAMETER          ( CONE = (1.0D+0, 0.0D+0) )
 129 *     ..
 130 *     .. Local Scalars ..
 131       LOGICAL            UPPER
 132       INTEGER            N1, N2, IINFO
 133       DOUBLE PRECISION   AJJ
 134 *     ..
 135 *     .. External Functions ..
 136       LOGICAL            LSAME, DISNAN
 137       EXTERNAL           LSAME, DISNAN
 138 *     ..
 139 *     .. External Subroutines ..
 140       EXTERNAL           ZHERK, ZTRSM, XERBLA
 141 *     ..
 142 *     .. Intrinsic Functions ..
 143       INTRINSIC          MAX, DBLE, SQRT
 144 *     ..
 145 *     .. Executable Statements ..
 146 *
 147 *     Test the input parameters
 148 *
 149       INFO = 0
 150       UPPER = LSAME( UPLO, 'U' )
 151       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 152          INFO = -1
 153       ELSE IF( N.LT.0 ) THEN
 154          INFO = -2
 155       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 156          INFO = -4
 157       END IF
 158       IF( INFO.NE.0 ) THEN
 159          CALL XERBLA( 'ZPOTRF2', -INFO )
 160          RETURN
 161       END IF
 162 *
 163 *     Quick return if possible
 164 *
 165       IF( N.EQ.0 )
 166      $   RETURN
 167 *
 168 *     N=1 case
 169 *
 170       IF( N.EQ.1 ) THEN
 171 *
 172 *        Test for non-positive-definiteness
 173 *
 174          AJJ = DBLE( A( 1, 1 ) )
 175          IF( AJJ.LE.ZERO.OR.DISNAN( AJJ ) ) THEN
 176             INFO = 1
 177             RETURN
 178          END IF
 179 *
 180 *        Factor
 181 *
 182          A( 1, 1 ) = SQRT( AJJ )
 183 *
 184 *     Use recursive code
 185 *
 186       ELSE
 187          N1 = N/2
 188          N2 = N-N1
 189 *
 190 *        Factor A11
 191 *
 192          CALL ZPOTRF2( UPLO, N1, A( 1, 1 ), LDA, IINFO )
 193          IF ( IINFO.NE.0 ) THEN
 194             INFO = IINFO
 195             RETURN
 196          END IF
 197 *
 198 *        Compute the Cholesky factorization A = U**H*U
 199 *
 200          IF( UPPER ) THEN
 201 *
 202 *           Update and scale A12
 203 *
 204             CALL ZTRSM( 'L', 'U', 'C', 'N', N1, N2, CONE,
 205      $                  A( 1, 1 ), LDA, A( 1, N1+1 ), LDA )
 206 *
 207 *           Update and factor A22
 208 *
 209             CALL ZHERK( UPLO, 'C', N2, N1, -ONE, A( 1, N1+1 ), LDA,
 210      $                  ONE, A( N1+1, N1+1 ), LDA )
 211             CALL ZPOTRF2( UPLO, N2, A( N1+1, N1+1 ), LDA, IINFO )
 212             IF ( IINFO.NE.0 ) THEN
 213                INFO = IINFO + N1
 214                RETURN
 215             END IF
 216 *
 217 *        Compute the Cholesky factorization A = L*L**H
 218 *
 219          ELSE
 220 *
 221 *           Update and scale A21
 222 *
 223             CALL ZTRSM( 'R', 'L', 'C', 'N', N2, N1, CONE,
 224      $                  A( 1, 1 ), LDA, A( N1+1, 1 ), LDA )
 225 *
 226 *           Update and factor A22
 227 *
 228             CALL ZHERK( UPLO, 'N', N2, N1, -ONE, A( N1+1, 1 ), LDA,
 229      $                  ONE, A( N1+1, N1+1 ), LDA )
 230             CALL ZPOTRF2( UPLO, N2, A( N1+1, N1+1 ), LDA, IINFO )
 231             IF ( IINFO.NE.0 ) THEN
 232                INFO = IINFO + N1
 233                RETURN
 234             END IF
 235          END IF
 236       END IF
 237       RETURN
 238 *
 239 *     End of ZPOTRF2
 240 *
 241       END