LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cpbtrf.f
Go to the documentation of this file.
00001 *> \brief \b CPBTRF
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CPBTRF + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cpbtrf.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cpbtrf.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cpbtrf.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CPBTRF( UPLO, N, KD, AB, LDAB, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       CHARACTER          UPLO
00025 *       INTEGER            INFO, KD, LDAB, N
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       COMPLEX            AB( LDAB, * )
00029 *       ..
00030 *  
00031 *
00032 *> \par Purpose:
00033 *  =============
00034 *>
00035 *> \verbatim
00036 *>
00037 *> CPBTRF computes the Cholesky factorization of a complex Hermitian
00038 *> positive definite band matrix A.
00039 *>
00040 *> The factorization has the form
00041 *>    A = U**H * U,  if UPLO = 'U', or
00042 *>    A = L  * L**H,  if UPLO = 'L',
00043 *> where U is an upper triangular matrix and L is lower triangular.
00044 *> \endverbatim
00045 *
00046 *  Arguments:
00047 *  ==========
00048 *
00049 *> \param[in] UPLO
00050 *> \verbatim
00051 *>          UPLO is CHARACTER*1
00052 *>          = 'U':  Upper triangle of A is stored;
00053 *>          = 'L':  Lower triangle of A is stored.
00054 *> \endverbatim
00055 *>
00056 *> \param[in] N
00057 *> \verbatim
00058 *>          N is INTEGER
00059 *>          The order of the matrix A.  N >= 0.
00060 *> \endverbatim
00061 *>
00062 *> \param[in] KD
00063 *> \verbatim
00064 *>          KD is INTEGER
00065 *>          The number of superdiagonals of the matrix A if UPLO = 'U',
00066 *>          or the number of subdiagonals if UPLO = 'L'.  KD >= 0.
00067 *> \endverbatim
00068 *>
00069 *> \param[in,out] AB
00070 *> \verbatim
00071 *>          AB is COMPLEX array, dimension (LDAB,N)
00072 *>          On entry, the upper or lower triangle of the Hermitian band
00073 *>          matrix A, stored in the first KD+1 rows of the array.  The
00074 *>          j-th column of A is stored in the j-th column of the array AB
00075 *>          as follows:
00076 *>          if UPLO = 'U', AB(kd+1+i-j,j) = A(i,j) for max(1,j-kd)<=i<=j;
00077 *>          if UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
00078 *>
00079 *>          On exit, if INFO = 0, the triangular factor U or L from the
00080 *>          Cholesky factorization A = U**H*U or A = L*L**H of the band
00081 *>          matrix A, in the same storage format as A.
00082 *> \endverbatim
00083 *>
00084 *> \param[in] LDAB
00085 *> \verbatim
00086 *>          LDAB is INTEGER
00087 *>          The leading dimension of the array AB.  LDAB >= KD+1.
00088 *> \endverbatim
00089 *>
00090 *> \param[out] INFO
00091 *> \verbatim
00092 *>          INFO is INTEGER
00093 *>          = 0:  successful exit
00094 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00095 *>          > 0:  if INFO = i, the leading minor of order i is not
00096 *>                positive definite, and the factorization could not be
00097 *>                completed.
00098 *> \endverbatim
00099 *
00100 *  Authors:
00101 *  ========
00102 *
00103 *> \author Univ. of Tennessee 
00104 *> \author Univ. of California Berkeley 
00105 *> \author Univ. of Colorado Denver 
00106 *> \author NAG Ltd. 
00107 *
00108 *> \date November 2011
00109 *
00110 *> \ingroup complexOTHERcomputational
00111 *
00112 *> \par Further Details:
00113 *  =====================
00114 *>
00115 *> \verbatim
00116 *>
00117 *>  The band storage scheme is illustrated by the following example, when
00118 *>  N = 6, KD = 2, and UPLO = 'U':
00119 *>
00120 *>  On entry:                       On exit:
00121 *>
00122 *>      *    *   a13  a24  a35  a46      *    *   u13  u24  u35  u46
00123 *>      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
00124 *>     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
00125 *>
00126 *>  Similarly, if UPLO = 'L' the format of A is as follows:
00127 *>
00128 *>  On entry:                       On exit:
00129 *>
00130 *>     a11  a22  a33  a44  a55  a66     l11  l22  l33  l44  l55  l66
00131 *>     a21  a32  a43  a54  a65   *      l21  l32  l43  l54  l65   *
00132 *>     a31  a42  a53  a64   *    *      l31  l42  l53  l64   *    *
00133 *>
00134 *>  Array elements marked * are not used by the routine.
00135 *> \endverbatim
00136 *
00137 *> \par Contributors:
00138 *  ==================
00139 *>
00140 *>  Peter Mayes and Giuseppe Radicati, IBM ECSEC, Rome, March 23, 1989
00141 *
00142 *  =====================================================================
00143       SUBROUTINE CPBTRF( UPLO, N, KD, AB, LDAB, INFO )
00144 *
00145 *  -- LAPACK computational routine (version 3.4.0) --
00146 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00147 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00148 *     November 2011
00149 *
00150 *     .. Scalar Arguments ..
00151       CHARACTER          UPLO
00152       INTEGER            INFO, KD, LDAB, N
00153 *     ..
00154 *     .. Array Arguments ..
00155       COMPLEX            AB( LDAB, * )
00156 *     ..
00157 *
00158 *  =====================================================================
00159 *
00160 *     .. Parameters ..
00161       REAL               ONE, ZERO
00162       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00163       COMPLEX            CONE
00164       PARAMETER          ( CONE = ( 1.0E+0, 0.0E+0 ) )
00165       INTEGER            NBMAX, LDWORK
00166       PARAMETER          ( NBMAX = 32, LDWORK = NBMAX+1 )
00167 *     ..
00168 *     .. Local Scalars ..
00169       INTEGER            I, I2, I3, IB, II, J, JJ, NB
00170 *     ..
00171 *     .. Local Arrays ..
00172       COMPLEX            WORK( LDWORK, NBMAX )
00173 *     ..
00174 *     .. External Functions ..
00175       LOGICAL            LSAME
00176       INTEGER            ILAENV
00177       EXTERNAL           LSAME, ILAENV
00178 *     ..
00179 *     .. External Subroutines ..
00180       EXTERNAL           CGEMM, CHERK, CPBTF2, CPOTF2, CTRSM, XERBLA
00181 *     ..
00182 *     .. Intrinsic Functions ..
00183       INTRINSIC          MIN
00184 *     ..
00185 *     .. Executable Statements ..
00186 *
00187 *     Test the input parameters.
00188 *
00189       INFO = 0
00190       IF( ( .NOT.LSAME( UPLO, 'U' ) ) .AND.
00191      $    ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN
00192          INFO = -1
00193       ELSE IF( N.LT.0 ) THEN
00194          INFO = -2
00195       ELSE IF( KD.LT.0 ) THEN
00196          INFO = -3
00197       ELSE IF( LDAB.LT.KD+1 ) THEN
00198          INFO = -5
00199       END IF
00200       IF( INFO.NE.0 ) THEN
00201          CALL XERBLA( 'CPBTRF', -INFO )
00202          RETURN
00203       END IF
00204 *
00205 *     Quick return if possible
00206 *
00207       IF( N.EQ.0 )
00208      $   RETURN
00209 *
00210 *     Determine the block size for this environment
00211 *
00212       NB = ILAENV( 1, 'CPBTRF', UPLO, N, KD, -1, -1 )
00213 *
00214 *     The block size must not exceed the semi-bandwidth KD, and must not
00215 *     exceed the limit set by the size of the local array WORK.
00216 *
00217       NB = MIN( NB, NBMAX )
00218 *
00219       IF( NB.LE.1 .OR. NB.GT.KD ) THEN
00220 *
00221 *        Use unblocked code
00222 *
00223          CALL CPBTF2( UPLO, N, KD, AB, LDAB, INFO )
00224       ELSE
00225 *
00226 *        Use blocked code
00227 *
00228          IF( LSAME( UPLO, 'U' ) ) THEN
00229 *
00230 *           Compute the Cholesky factorization of a Hermitian band
00231 *           matrix, given the upper triangle of the matrix in band
00232 *           storage.
00233 *
00234 *           Zero the upper triangle of the work array.
00235 *
00236             DO 20 J = 1, NB
00237                DO 10 I = 1, J - 1
00238                   WORK( I, J ) = ZERO
00239    10          CONTINUE
00240    20       CONTINUE
00241 *
00242 *           Process the band matrix one diagonal block at a time.
00243 *
00244             DO 70 I = 1, N, NB
00245                IB = MIN( NB, N-I+1 )
00246 *
00247 *              Factorize the diagonal block
00248 *
00249                CALL CPOTF2( UPLO, IB, AB( KD+1, I ), LDAB-1, II )
00250                IF( II.NE.0 ) THEN
00251                   INFO = I + II - 1
00252                   GO TO 150
00253                END IF
00254                IF( I+IB.LE.N ) THEN
00255 *
00256 *                 Update the relevant part of the trailing submatrix.
00257 *                 If A11 denotes the diagonal block which has just been
00258 *                 factorized, then we need to update the remaining
00259 *                 blocks in the diagram:
00260 *
00261 *                    A11   A12   A13
00262 *                          A22   A23
00263 *                                A33
00264 *
00265 *                 The numbers of rows and columns in the partitioning
00266 *                 are IB, I2, I3 respectively. The blocks A12, A22 and
00267 *                 A23 are empty if IB = KD. The upper triangle of A13
00268 *                 lies outside the band.
00269 *
00270                   I2 = MIN( KD-IB, N-I-IB+1 )
00271                   I3 = MIN( IB, N-I-KD+1 )
00272 *
00273                   IF( I2.GT.0 ) THEN
00274 *
00275 *                    Update A12
00276 *
00277                      CALL CTRSM( 'Left', 'Upper', 'Conjugate transpose',
00278      $                           'Non-unit', IB, I2, CONE,
00279      $                           AB( KD+1, I ), LDAB-1,
00280      $                           AB( KD+1-IB, I+IB ), LDAB-1 )
00281 *
00282 *                    Update A22
00283 *
00284                      CALL CHERK( 'Upper', 'Conjugate transpose', I2, IB,
00285      $                           -ONE, AB( KD+1-IB, I+IB ), LDAB-1, ONE,
00286      $                           AB( KD+1, I+IB ), LDAB-1 )
00287                   END IF
00288 *
00289                   IF( I3.GT.0 ) THEN
00290 *
00291 *                    Copy the lower triangle of A13 into the work array.
00292 *
00293                      DO 40 JJ = 1, I3
00294                         DO 30 II = JJ, IB
00295                            WORK( II, JJ ) = AB( II-JJ+1, JJ+I+KD-1 )
00296    30                   CONTINUE
00297    40                CONTINUE
00298 *
00299 *                    Update A13 (in the work array).
00300 *
00301                      CALL CTRSM( 'Left', 'Upper', 'Conjugate transpose',
00302      $                           'Non-unit', IB, I3, CONE,
00303      $                           AB( KD+1, I ), LDAB-1, WORK, LDWORK )
00304 *
00305 *                    Update A23
00306 *
00307                      IF( I2.GT.0 )
00308      $                  CALL CGEMM( 'Conjugate transpose',
00309      $                              'No transpose', I2, I3, IB, -CONE,
00310      $                              AB( KD+1-IB, I+IB ), LDAB-1, WORK,
00311      $                              LDWORK, CONE, AB( 1+IB, I+KD ),
00312      $                              LDAB-1 )
00313 *
00314 *                    Update A33
00315 *
00316                      CALL CHERK( 'Upper', 'Conjugate transpose', I3, IB,
00317      $                           -ONE, WORK, LDWORK, ONE,
00318      $                           AB( KD+1, I+KD ), LDAB-1 )
00319 *
00320 *                    Copy the lower triangle of A13 back into place.
00321 *
00322                      DO 60 JJ = 1, I3
00323                         DO 50 II = JJ, IB
00324                            AB( II-JJ+1, JJ+I+KD-1 ) = WORK( II, JJ )
00325    50                   CONTINUE
00326    60                CONTINUE
00327                   END IF
00328                END IF
00329    70       CONTINUE
00330          ELSE
00331 *
00332 *           Compute the Cholesky factorization of a Hermitian band
00333 *           matrix, given the lower triangle of the matrix in band
00334 *           storage.
00335 *
00336 *           Zero the lower triangle of the work array.
00337 *
00338             DO 90 J = 1, NB
00339                DO 80 I = J + 1, NB
00340                   WORK( I, J ) = ZERO
00341    80          CONTINUE
00342    90       CONTINUE
00343 *
00344 *           Process the band matrix one diagonal block at a time.
00345 *
00346             DO 140 I = 1, N, NB
00347                IB = MIN( NB, N-I+1 )
00348 *
00349 *              Factorize the diagonal block
00350 *
00351                CALL CPOTF2( UPLO, IB, AB( 1, I ), LDAB-1, II )
00352                IF( II.NE.0 ) THEN
00353                   INFO = I + II - 1
00354                   GO TO 150
00355                END IF
00356                IF( I+IB.LE.N ) THEN
00357 *
00358 *                 Update the relevant part of the trailing submatrix.
00359 *                 If A11 denotes the diagonal block which has just been
00360 *                 factorized, then we need to update the remaining
00361 *                 blocks in the diagram:
00362 *
00363 *                    A11
00364 *                    A21   A22
00365 *                    A31   A32   A33
00366 *
00367 *                 The numbers of rows and columns in the partitioning
00368 *                 are IB, I2, I3 respectively. The blocks A21, A22 and
00369 *                 A32 are empty if IB = KD. The lower triangle of A31
00370 *                 lies outside the band.
00371 *
00372                   I2 = MIN( KD-IB, N-I-IB+1 )
00373                   I3 = MIN( IB, N-I-KD+1 )
00374 *
00375                   IF( I2.GT.0 ) THEN
00376 *
00377 *                    Update A21
00378 *
00379                      CALL CTRSM( 'Right', 'Lower',
00380      $                           'Conjugate transpose', 'Non-unit', I2,
00381      $                           IB, CONE, AB( 1, I ), LDAB-1,
00382      $                           AB( 1+IB, I ), LDAB-1 )
00383 *
00384 *                    Update A22
00385 *
00386                      CALL CHERK( 'Lower', 'No transpose', I2, IB, -ONE,
00387      $                           AB( 1+IB, I ), LDAB-1, ONE,
00388      $                           AB( 1, I+IB ), LDAB-1 )
00389                   END IF
00390 *
00391                   IF( I3.GT.0 ) THEN
00392 *
00393 *                    Copy the upper triangle of A31 into the work array.
00394 *
00395                      DO 110 JJ = 1, IB
00396                         DO 100 II = 1, MIN( JJ, I3 )
00397                            WORK( II, JJ ) = AB( KD+1-JJ+II, JJ+I-1 )
00398   100                   CONTINUE
00399   110                CONTINUE
00400 *
00401 *                    Update A31 (in the work array).
00402 *
00403                      CALL CTRSM( 'Right', 'Lower',
00404      $                           'Conjugate transpose', 'Non-unit', I3,
00405      $                           IB, CONE, AB( 1, I ), LDAB-1, WORK,
00406      $                           LDWORK )
00407 *
00408 *                    Update A32
00409 *
00410                      IF( I2.GT.0 )
00411      $                  CALL CGEMM( 'No transpose',
00412      $                              'Conjugate transpose', I3, I2, IB,
00413      $                              -CONE, WORK, LDWORK, AB( 1+IB, I ),
00414      $                              LDAB-1, CONE, AB( 1+KD-IB, I+IB ),
00415      $                              LDAB-1 )
00416 *
00417 *                    Update A33
00418 *
00419                      CALL CHERK( 'Lower', 'No transpose', I3, IB, -ONE,
00420      $                           WORK, LDWORK, ONE, AB( 1, I+KD ),
00421      $                           LDAB-1 )
00422 *
00423 *                    Copy the upper triangle of A31 back into place.
00424 *
00425                      DO 130 JJ = 1, IB
00426                         DO 120 II = 1, MIN( JJ, I3 )
00427                            AB( KD+1-JJ+II, JJ+I-1 ) = WORK( II, JJ )
00428   120                   CONTINUE
00429   130                CONTINUE
00430                   END IF
00431                END IF
00432   140       CONTINUE
00433          END IF
00434       END IF
00435       RETURN
00436 *
00437   150 CONTINUE
00438       RETURN
00439 *
00440 *     End of CPBTRF
00441 *
00442       END
 All Files Functions