LAPACK: SRC/zpbtrf.f Source File

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