SRC/cpbtrf.f

   1 *> \brief \b CPBTRF
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CPBTRF + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cpbtrf.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cpbtrf.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cpbtrf.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CPBTRF( 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            AB( LDAB, * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> CPBTRF 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 and L is lower triangular.
  44 *> \endverbatim
  45 *
  46 *  Arguments:
  47 *  ==========
  48 *
  49 *> \param[in] UPLO
  50 *> \verbatim
  51 *>          UPLO is CHARACTER*1
  52 *>          = 'U':  Upper triangle of A is stored;
  53 *>          = 'L':  Lower triangle of A is stored.
  54 *> \endverbatim
  55 *>
  56 *> \param[in] N
  57 *> \verbatim
  58 *>          N is INTEGER
  59 *>          The order of the matrix A.  N >= 0.
  60 *> \endverbatim
  61 *>
  62 *> \param[in] KD
  63 *> \verbatim
  64 *>          KD is INTEGER
  65 *>          The number of superdiagonals of the matrix A if UPLO = 'U',
  66 *>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
  67 *> \endverbatim
  68 *>
  69 *> \param[in,out] AB
  70 *> \verbatim
  71 *>          AB is COMPLEX array, dimension (LDAB,N)
  72 *>          On entry, the upper or lower triangle of the Hermitian band
  73 *>          matrix A, stored in the first KD+1 rows of the array.  The
  74 *>          j-th column of A is stored in the j-th column of the array AB
  75 *>          as follows:
  76 *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
  77 *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
  78 *>
  79 *>          On exit, if INFO = 0, the triangular factor U or L from the
  80 *>          Cholesky factorization A = U**H*U or A = L*L**H of the band
  81 *>          matrix A, in the same storage format as A.
  82 *> \endverbatim
  83 *>
  84 *> \param[in] LDAB
  85 *> \verbatim
  86 *>          LDAB is INTEGER
  87 *>          The leading dimension of the array AB.  LDAB >= KD+1.
  88 *> \endverbatim
  89 *>
  90 *> \param[out] INFO
  91 *> \verbatim
  92 *>          INFO is INTEGER
  93 *>          = 0:  successful exit
  94 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  95 *>          > 0:  if INFO = i, the leading minor of order i is not
  96 *>                positive definite, and the factorization could not be
  97 *>                completed.
  98 *> \endverbatim
  99 *
 100 *  Authors:
 101 *  ========
 102 *
 103 *> \author Univ. of Tennessee
 104 *> \author Univ. of California Berkeley
 105 *> \author Univ. of Colorado Denver
 106 *> \author NAG Ltd.
 107 *
 108 *> \date November 2011
 109 *
 110 *> \ingroup complexOTHERcomputational
 111 *
 112 *> \par Further Details:
 113 *  =====================
 114 *>
 115 *> \verbatim
 116 *>
 117 *>  The band storage scheme is illustrated by the following example, when
 118 *>  N = 6, KD = 2, and UPLO = 'U':
 119 *>
 120 *>  On entry:                       On exit:
 121 *>
 122 *>      *    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46
 123 *>      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
 124 *>     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
 125 *>
 126 *>  Similarly, if UPLO = 'L' the format of A is as follows:
 127 *>
 128 *>  On entry:                       On exit:
 129 *>
 130 *>     a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66
 131 *>     a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *
 132 *>     a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *
 133 *>
 134 *>  Array elements marked * are not used by the routine.
 135 *> \endverbatim
 136 *
 137 *> \par Contributors:
 138 *  ==================
 139 *>
 140 *>  Peter Mayes and Giuseppe Radicati, IBM ECSEC, Rome, March 23, 1989
 141 *
 142 *  =====================================================================
 143       SUBROUTINE CPBTRF( UPLO, N, KD, AB, LDAB, INFO )
 144 *
 145 *  -- LAPACK computational routine (version 3.4.0) --
 146 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 147 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 148 *     November 2011
 149 *
 150 *     .. Scalar Arguments ..
 151       CHARACTER          UPLO
 152       INTEGER            INFO, KD, LDAB, N
 153 *     ..
 154 *     .. Array Arguments ..
 155       COMPLEX            AB( LDAB, * )
 156 *     ..
 157 *
 158 *  =====================================================================
 159 *
 160 *     .. Parameters ..
 161       REAL               ONE, ZERO
 162       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
 163       COMPLEX            CONE
 164       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
 165       INTEGER            NBMAX, LDWORK
 166       PARAMETER          ( NBMAX = 32, LDWORK = NBMAX+1 )
 167 *     ..
 168 *     .. Local Scalars ..
 169       INTEGER            I, I2, I3, IB, II, J, JJ, NB
 170 *     ..
 171 *     .. Local Arrays ..
 172       COMPLEX            WORK( LDWORK, NBMAX )
 173 *     ..
 174 *     .. External Functions ..
 175       LOGICAL            LSAME
 176       INTEGER            ILAENV
 177       EXTERNAL           LSAME, ILAENV
 178 *     ..
 179 *     .. External Subroutines ..
 180       EXTERNAL           CGEMM, CHERK, CPBTF2, CPOTF2, CTRSM, XERBLA
 181 *     ..
 182 *     .. Intrinsic Functions ..
 183       INTRINSIC          MIN
 184 *     ..
 185 *     .. Executable Statements ..
 186 *
 187 *     Test the input parameters.
 188 *
 189       INFO = 0
 190       IF( ( .NOT.LSAME( UPLO, 'U' ) ) .AND.
 191      $    ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN
 192          INFO = -1
 193       ELSE IF( N.LT.0 ) THEN
 194          INFO = -2
 195       ELSE IF( KD.LT.0 ) THEN
 196          INFO = -3
 197       ELSE IF( LDAB.LT.KD+1 ) THEN
 198          INFO = -5
 199       END IF
 200       IF( INFO.NE.0 ) THEN
 201          CALL XERBLA( 'CPBTRF', -INFO )
 202          RETURN
 203       END IF
 204 *
 205 *     Quick return if possible
 206 *
 207       IF( N.EQ.0 )
 208      $   RETURN
 209 *
 210 *     Determine the block size for this environment
 211 *
 212       NB = ILAENV( 1, 'CPBTRF', UPLO, N, KD, -1, -1 )
 213 *
 214 *     The block size must not exceed the semi-bandwidth KD, and must not
 215 *     exceed the limit set by the size of the local array WORK.
 216 *
 217       NB = MIN( NB, NBMAX )
 218 *
 219       IF( NB.LE.1 .OR. NB.GT.KD ) THEN
 220 *
 221 *        Use unblocked code
 222 *
 223          CALL CPBTF2( UPLO, N, KD, AB, LDAB, INFO )
 224       ELSE
 225 *
 226 *        Use blocked code
 227 *
 228          IF( LSAME( UPLO, 'U' ) ) THEN
 229 *
 230 *           Compute the Cholesky factorization of a Hermitian band
 231 *           matrix, given the upper triangle of the matrix in band
 232 *           storage.
 233 *
 234 *           Zero the upper triangle of the work array.
 235 *
 236             DO 20 J = 1, NB
 237                DO 10 I = 1, J - 1
 238                   WORK( I, J ) = ZERO
 239    10          CONTINUE
 240    20       CONTINUE
 241 *
 242 *           Process the band matrix one diagonal block at a time.
 243 *
 244             DO 70 I = 1, N, NB
 245                IB = MIN( NB, N-I+1 )
 246 *
 247 *              Factorize the diagonal block
 248 *
 249                CALL CPOTF2( UPLO, IB, AB( KD+1, I ), LDAB-1, II )
 250                IF( II.NE.0 ) THEN
 251                   INFO = I + II - 1
 252                   GO TO 150
 253                END IF
 254                IF( I+IB.LE.N ) THEN
 255 *
 256 *                 Update the relevant part of the trailing submatrix.
 257 *                 If A11 denotes the diagonal block which has just been
 258 *                 factorized, then we need to update the remaining
 259 *                 blocks in the diagram:
 260 *
 261 *                    A11   A12   A13
 262 *                          A22   A23
 263 *                                A33
 264 *
 265 *                 The numbers of rows and columns in the partitioning
 266 *                 are IB, I2, I3 respectively. The blocks A12, A22 and
 267 *                 A23 are empty if IB = KD. The upper triangle of A13
 268 *                 lies outside the band.
 269 *
 270                   I2 = MIN( KD-IB, N-I-IB+1 )
 271                   I3 = MIN( IB, N-I-KD+1 )
 272 *
 273                   IF( I2.GT.0 ) THEN
 274 *
 275 *                    Update A12
 276 *
 277                      CALL CTRSM( 'Left', 'Upper', 'Conjugate transpose',
 278      $                           'Non-unit', IB, I2, CONE,
 279      $                           AB( KD+1, I ), LDAB-1,
 280      $                           AB( KD+1-IB, I+IB ), LDAB-1 )
 281 *
 282 *                    Update A22
 283 *
 284                      CALL CHERK( 'Upper', 'Conjugate transpose', I2, IB,
 285      $                           -ONE, AB( KD+1-IB, I+IB ), LDAB-1, ONE,
 286      $                           AB( KD+1, I+IB ), LDAB-1 )
 287                   END IF
 288 *
 289                   IF( I3.GT.0 ) THEN
 290 *
 291 *                    Copy the lower triangle of A13 into the work array.
 292 *
 293                      DO 40 JJ = 1, I3
 294                         DO 30 II = JJ, IB
 295                            WORK( II, JJ ) = AB( II-JJ+1, JJ+I+KD-1 )
 296    30                   CONTINUE
 297    40                CONTINUE
 298 *
 299 *                    Update A13 (in the work array).
 300 *
 301                      CALL CTRSM( 'Left', 'Upper', 'Conjugate transpose',
 302      $                           'Non-unit', IB, I3, CONE,
 303      $                           AB( KD+1, I ), LDAB-1, WORK, LDWORK )
 304 *
 305 *                    Update A23
 306 *
 307                      IF( I2.GT.0 )
 308      $                  CALL CGEMM( 'Conjugate transpose',
 309      $                              'No transpose', I2, I3, IB, -CONE,
 310      $                              AB( KD+1-IB, I+IB ), LDAB-1, WORK,
 311      $                              LDWORK, CONE, AB( 1+IB, I+KD ),
 312      $                              LDAB-1 )
 313 *
 314 *                    Update A33
 315 *
 316                      CALL CHERK( 'Upper', 'Conjugate transpose', I3, IB,
 317      $                           -ONE, WORK, LDWORK, ONE,
 318      $                           AB( KD+1, I+KD ), LDAB-1 )
 319 *
 320 *                    Copy the lower triangle of A13 back into place.
 321 *
 322                      DO 60 JJ = 1, I3
 323                         DO 50 II = JJ, IB
 324                            AB( II-JJ+1, JJ+I+KD-1 ) = WORK( II, JJ )
 325    50                   CONTINUE
 326    60                CONTINUE
 327                   END IF
 328                END IF
 329    70       CONTINUE
 330          ELSE
 331 *
 332 *           Compute the Cholesky factorization of a Hermitian band
 333 *           matrix, given the lower triangle of the matrix in band
 334 *           storage.
 335 *
 336 *           Zero the lower triangle of the work array.
 337 *
 338             DO 90 J = 1, NB
 339                DO 80 I = J + 1, NB
 340                   WORK( I, J ) = ZERO
 341    80          CONTINUE
 342    90       CONTINUE
 343 *
 344 *           Process the band matrix one diagonal block at a time.
 345 *
 346             DO 140 I = 1, N, NB
 347                IB = MIN( NB, N-I+1 )
 348 *
 349 *              Factorize the diagonal block
 350 *
 351                CALL CPOTF2( UPLO, IB, AB( 1, I ), LDAB-1, II )
 352                IF( II.NE.0 ) THEN
 353                   INFO = I + II - 1
 354                   GO TO 150
 355                END IF
 356                IF( I+IB.LE.N ) THEN
 357 *
 358 *                 Update the relevant part of the trailing submatrix.
 359 *                 If A11 denotes the diagonal block which has just been
 360 *                 factorized, then we need to update the remaining
 361 *                 blocks in the diagram:
 362 *
 363 *                    A11
 364 *                    A21   A22
 365 *                    A31   A32   A33
 366 *
 367 *                 The numbers of rows and columns in the partitioning
 368 *                 are IB, I2, I3 respectively. The blocks A21, A22 and
 369 *                 A32 are empty if IB = KD. The lower triangle of A31
 370 *                 lies outside the band.
 371 *
 372                   I2 = MIN( KD-IB, N-I-IB+1 )
 373                   I3 = MIN( IB, N-I-KD+1 )
 374 *
 375                   IF( I2.GT.0 ) THEN
 376 *
 377 *                    Update A21
 378 *
 379                      CALL CTRSM( 'Right', 'Lower',
 380      $                           'Conjugate transpose', 'Non-unit', I2,
 381      $                           IB, CONE, AB( 1, I ), LDAB-1,
 382      $                           AB( 1+IB, I ), LDAB-1 )
 383 *
 384 *                    Update A22
 385 *
 386                      CALL CHERK( 'Lower', 'No transpose', I2, IB, -ONE,
 387      $                           AB( 1+IB, I ), LDAB-1, ONE,
 388      $                           AB( 1, I+IB ), LDAB-1 )
 389                   END IF
 390 *
 391                   IF( I3.GT.0 ) THEN
 392 *
 393 *                    Copy the upper triangle of A31 into the work array.
 394 *
 395                      DO 110 JJ = 1, IB
 396                         DO 100 II = 1, MIN( JJ, I3 )
 397                            WORK( II, JJ ) = AB( KD+1-JJ+II, JJ+I-1 )
 398   100                   CONTINUE
 399   110                CONTINUE
 400 *
 401 *                    Update A31 (in the work array).
 402 *
 403                      CALL CTRSM( 'Right', 'Lower',
 404      $                           'Conjugate transpose', 'Non-unit', I3,
 405      $                           IB, CONE, AB( 1, I ), LDAB-1, WORK,
 406      $                           LDWORK )
 407 *
 408 *                    Update A32
 409 *
 410                      IF( I2.GT.0 )
 411      $                  CALL CGEMM( 'No transpose',
 412      $                              'Conjugate transpose', I3, I2, IB,
 413      $                              -CONE, WORK, LDWORK, AB( 1+IB, I ),
 414      $                              LDAB-1, CONE, AB( 1+KD-IB, I+IB ),
 415      $                              LDAB-1 )
 416 *
 417 *                    Update A33
 418 *
 419                      CALL CHERK( 'Lower', 'No transpose', I3, IB, -ONE,
 420      $                           WORK, LDWORK, ONE, AB( 1, I+KD ),
 421      $                           LDAB-1 )
 422 *
 423 *                    Copy the upper triangle of A31 back into place.
 424 *
 425                      DO 130 JJ = 1, IB
 426                         DO 120 II = 1, MIN( JJ, I3 )
 427                            AB( KD+1-JJ+II, JJ+I-1 ) = WORK( II, JJ )
 428   120                   CONTINUE
 429   130                CONTINUE
 430                   END IF
 431                END IF
 432   140       CONTINUE
 433          END IF
 434       END IF
 435       RETURN
 436 *
 437   150 CONTINUE
 438       RETURN
 439 *
 440 *     End of CPBTRF
 441 *
 442       END