SRC/zpbtf2.f

   1 *> \brief \b ZPBTF2 computes the Cholesky factorization of a symmetric/Hermitian positive definite band 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 ZPBTF2 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zpbtf2.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zpbtf2.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zpbtf2.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZPBTF2( UPLO, N, KD, AB, LDAB, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          UPLO
  25 *       INTEGER            INFO, KD, LDAB, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       COMPLEX*16         AB( LDAB, * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> ZPBTF2 computes the Cholesky factorization of a complex Hermitian
  38 *> positive definite band 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, U**H is the conjugate transpose
  44 *> of U, and L is lower triangular.
  45 *>
  46 *> This is the unblocked version of the algorithm, calling Level 2 BLAS.
  47 *> \endverbatim
  48 *
  49 *  Arguments:
  50 *  ==========
  51 *
  52 *> \param[in] UPLO
  53 *> \verbatim
  54 *>          UPLO is CHARACTER*1
  55 *>          Specifies whether the upper or lower triangular part of the
  56 *>          Hermitian matrix A is stored:
  57 *>          = 'U':  Upper triangular
  58 *>          = 'L':  Lower triangular
  59 *> \endverbatim
  60 *>
  61 *> \param[in] N
  62 *> \verbatim
  63 *>          N is INTEGER
  64 *>          The order of the matrix A.  N >= 0.
  65 *> \endverbatim
  66 *>
  67 *> \param[in] KD
  68 *> \verbatim
  69 *>          KD is INTEGER
  70 *>          The number of super-diagonals of the matrix A if UPLO = 'U',
  71 *>          or the number of sub-diagonals if UPLO = 'L'.  KD >= 0.
  72 *> \endverbatim
  73 *>
  74 *> \param[in,out] AB
  75 *> \verbatim
  76 *>          AB is COMPLEX*16 array, dimension (LDAB,N)
  77 *>          On entry, the upper or lower triangle of the Hermitian band
  78 *>          matrix A, stored in the first KD+1 rows of the array.  The
  79 *>          j-th column of A is stored in the j-th column of the array AB
  80 *>          as follows:
  81 *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
  82 *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
  83 *>
  84 *>          On exit, if INFO = 0, the triangular factor U or L from the
  85 *>          Cholesky factorization A = U**H *U or A = L*L**H of the band
  86 *>          matrix A, in the same storage format as A.
  87 *> \endverbatim
  88 *>
  89 *> \param[in] LDAB
  90 *> \verbatim
  91 *>          LDAB is INTEGER
  92 *>          The leading dimension of the array AB.  LDAB >= KD+1.
  93 *> \endverbatim
  94 *>
  95 *> \param[out] INFO
  96 *> \verbatim
  97 *>          INFO is INTEGER
  98 *>          = 0: successful exit
  99 *>          < 0: if INFO = -k, the k-th argument had an illegal value
 100 *>          > 0: if INFO = k, the leading minor of order k is not
 101 *>               positive definite, and the factorization could not be
 102 *>               completed.
 103 *> \endverbatim
 104 *
 105 *  Authors:
 106 *  ========
 107 *
 108 *> \author Univ. of Tennessee
 109 *> \author Univ. of California Berkeley
 110 *> \author Univ. of Colorado Denver
 111 *> \author NAG Ltd.
 112 *
 113 *> \date September 2012
 114 *
 115 *> \ingroup complex16OTHERcomputational
 116 *
 117 *> \par Further Details:
 118 *  =====================
 119 *>
 120 *> \verbatim
 121 *>
 122 *>  The band storage scheme is illustrated by the following example, when
 123 *>  N = 6, KD = 2, and UPLO = 'U':
 124 *>
 125 *>  On entry:                       On exit:
 126 *>
 127 *>      *    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46
 128 *>      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
 129 *>     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
 130 *>
 131 *>  Similarly, if UPLO = 'L' the format of A is as follows:
 132 *>
 133 *>  On entry:                       On exit:
 134 *>
 135 *>     a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66
 136 *>     a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *
 137 *>     a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *
 138 *>
 139 *>  Array elements marked * are not used by the routine.
 140 *> \endverbatim
 141 *>
 142 *  =====================================================================
 143       SUBROUTINE ZPBTF2( UPLO, N, KD, AB, LDAB, INFO )
 144 *
 145 *  -- LAPACK computational routine (version 3.4.2) --
 146 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 147 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 148 *     September 2012
 149 *
 150 *     .. Scalar Arguments ..
 151       CHARACTER          UPLO
 152       INTEGER            INFO, KD, LDAB, N
 153 *     ..
 154 *     .. Array Arguments ..
 155       COMPLEX*16         AB( LDAB, * )
 156 *     ..
 157 *
 158 *  =====================================================================
 159 *
 160 *     .. Parameters ..
 161       DOUBLE PRECISION   ONE, ZERO
 162       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 163 *     ..
 164 *     .. Local Scalars ..
 165       LOGICAL            UPPER
 166       INTEGER            J, KLD, KN
 167       DOUBLE PRECISION   AJJ
 168 *     ..
 169 *     .. External Functions ..
 170       LOGICAL            LSAME
 171       EXTERNAL           LSAME
 172 *     ..
 173 *     .. External Subroutines ..
 174       EXTERNAL           XERBLA, ZDSCAL, ZHER, ZLACGV
 175 *     ..
 176 *     .. Intrinsic Functions ..
 177       INTRINSIC          DBLE, MAX, MIN, SQRT
 178 *     ..
 179 *     .. Executable Statements ..
 180 *
 181 *     Test the input parameters.
 182 *
 183       INFO = 0
 184       UPPER = LSAME( UPLO, 'U' )
 185       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
 186          INFO = -1
 187       ELSE IF( N.LT.0 ) THEN
 188          INFO = -2
 189       ELSE IF( KD.LT.0 ) THEN
 190          INFO = -3
 191       ELSE IF( LDAB.LT.KD+1 ) THEN
 192          INFO = -5
 193       END IF
 194       IF( INFO.NE.0 ) THEN
 195          CALL XERBLA( 'ZPBTF2', -INFO )
 196          RETURN
 197       END IF
 198 *
 199 *     Quick return if possible
 200 *
 201       IF( N.EQ.0 )
 202      $   RETURN
 203 *
 204       KLD = MAX( 1, LDAB-1 )
 205 *
 206       IF( UPPER ) THEN
 207 *
 208 *        Compute the Cholesky factorization A = U**H * U.
 209 *
 210          DO 10 J = 1, N
 211 *
 212 *           Compute U(J,J) and test for non-positive-definiteness.
 213 *
 214             AJJ = DBLE( AB( KD+1, J ) )
 215             IF( AJJ.LE.ZERO ) THEN
 216                AB( KD+1, J ) = AJJ
 217                GO TO 30
 218             END IF
 219             AJJ = SQRT( AJJ )
 220             AB( KD+1, J ) = AJJ
 221 *
 222 *           Compute elements J+1:J+KN of row J and update the
 223 *           trailing submatrix within the band.
 224 *
 225             KN = MIN( KD, N-J )
 226             IF( KN.GT.0 ) THEN
 227                CALL ZDSCAL( KN, ONE / AJJ, AB( KD, J+1 ), KLD )
 228                CALL ZLACGV( KN, AB( KD, J+1 ), KLD )
 229                CALL ZHER( 'Upper', KN, -ONE, AB( KD, J+1 ), KLD,
 230      $                    AB( KD+1, J+1 ), KLD )
 231                CALL ZLACGV( KN, AB( KD, J+1 ), KLD )
 232             END IF
 233    10    CONTINUE
 234       ELSE
 235 *
 236 *        Compute the Cholesky factorization A = L*L**H.
 237 *
 238          DO 20 J = 1, N
 239 *
 240 *           Compute L(J,J) and test for non-positive-definiteness.
 241 *
 242             AJJ = DBLE( AB( 1, J ) )
 243             IF( AJJ.LE.ZERO ) THEN
 244                AB( 1, J ) = AJJ
 245                GO TO 30
 246             END IF
 247             AJJ = SQRT( AJJ )
 248             AB( 1, J ) = AJJ
 249 *
 250 *           Compute elements J+1:J+KN of column J and update the
 251 *           trailing submatrix within the band.
 252 *
 253             KN = MIN( KD, N-J )
 254             IF( KN.GT.0 ) THEN
 255                CALL ZDSCAL( KN, ONE / AJJ, AB( 2, J ), 1 )
 256                CALL ZHER( 'Lower', KN, -ONE, AB( 2, J ), 1,
 257      $                    AB( 1, J+1 ), KLD )
 258             END IF
 259    20    CONTINUE
 260       END IF
 261       RETURN
 262 *
 263    30 CONTINUE
 264       INFO = J
 265       RETURN
 266 *
 267 *     End of ZPBTF2
 268 *
 269       END