LAPACK 3.3.1 Linear Algebra PACKage

# spbtrf.f

Go to the documentation of this file.
00001       SUBROUTINE SPBTRF( 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       REAL               AB( LDAB, * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  SPBTRF computes the Cholesky factorization of a real symmetric
00020 *  positive definite band matrix A.
00021 *
00022 *  The factorization has the form
00023 *     A = U**T * U,  if UPLO = 'U', or
00024 *     A = L  * L**T,  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) REAL array, dimension (LDAB,N)
00042 *          On entry, the upper or lower triangle of the symmetric 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**T*U or A = L*L**T 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       REAL               ONE, ZERO
00092       PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
00093       INTEGER            NBMAX, LDWORK
00094       PARAMETER          ( NBMAX = 32, LDWORK = NBMAX+1 )
00095 *     ..
00096 *     .. Local Scalars ..
00097       INTEGER            I, I2, I3, IB, II, J, JJ, NB
00098 *     ..
00099 *     .. Local Arrays ..
00100       REAL               WORK( LDWORK, NBMAX )
00101 *     ..
00102 *     .. External Functions ..
00103       LOGICAL            LSAME
00104       INTEGER            ILAENV
00105       EXTERNAL           LSAME, ILAENV
00106 *     ..
00107 *     .. External Subroutines ..
00108       EXTERNAL           SGEMM, SPBTF2, SPOTF2, SSYRK, STRSM, XERBLA
00109 *     ..
00110 *     .. Intrinsic Functions ..
00111       INTRINSIC          MIN
00112 *     ..
00113 *     .. Executable Statements ..
00114 *
00115 *     Test the input parameters.
00116 *
00117       INFO = 0
00118       IF( ( .NOT.LSAME( UPLO, 'U' ) ) .AND.
00119      \$    ( .NOT.LSAME( UPLO, 'L' ) ) ) THEN
00120          INFO = -1
00121       ELSE IF( N.LT.0 ) THEN
00122          INFO = -2
00123       ELSE IF( KD.LT.0 ) THEN
00124          INFO = -3
00125       ELSE IF( LDAB.LT.KD+1 ) THEN
00126          INFO = -5
00127       END IF
00128       IF( INFO.NE.0 ) THEN
00129          CALL XERBLA( 'SPBTRF', -INFO )
00130          RETURN
00131       END IF
00132 *
00133 *     Quick return if possible
00134 *
00135       IF( N.EQ.0 )
00136      \$   RETURN
00137 *
00138 *     Determine the block size for this environment
00139 *
00140       NB = ILAENV( 1, 'SPBTRF', UPLO, N, KD, -1, -1 )
00141 *
00142 *     The block size must not exceed the semi-bandwidth KD, and must not
00143 *     exceed the limit set by the size of the local array WORK.
00144 *
00145       NB = MIN( NB, NBMAX )
00146 *
00147       IF( NB.LE.1 .OR. NB.GT.KD ) THEN
00148 *
00149 *        Use unblocked code
00150 *
00151          CALL SPBTF2( UPLO, N, KD, AB, LDAB, INFO )
00152       ELSE
00153 *
00154 *        Use blocked code
00155 *
00156          IF( LSAME( UPLO, 'U' ) ) THEN
00157 *
00158 *           Compute the Cholesky factorization of a symmetric band
00159 *           matrix, given the upper triangle of the matrix in band
00160 *           storage.
00161 *
00162 *           Zero the upper triangle of the work array.
00163 *
00164             DO 20 J = 1, NB
00165                DO 10 I = 1, J - 1
00166                   WORK( I, J ) = ZERO
00167    10          CONTINUE
00168    20       CONTINUE
00169 *
00170 *           Process the band matrix one diagonal block at a time.
00171 *
00172             DO 70 I = 1, N, NB
00173                IB = MIN( NB, N-I+1 )
00174 *
00175 *              Factorize the diagonal block
00176 *
00177                CALL SPOTF2( UPLO, IB, AB( KD+1, I ), LDAB-1, II )
00178                IF( II.NE.0 ) THEN
00179                   INFO = I + II - 1
00180                   GO TO 150
00181                END IF
00182                IF( I+IB.LE.N ) THEN
00183 *
00184 *                 Update the relevant part of the trailing submatrix.
00185 *                 If A11 denotes the diagonal block which has just been
00186 *                 factorized, then we need to update the remaining
00187 *                 blocks in the diagram:
00188 *
00189 *                    A11   A12   A13
00190 *                          A22   A23
00191 *                                A33
00192 *
00193 *                 The numbers of rows and columns in the partitioning
00194 *                 are IB, I2, I3 respectively. The blocks A12, A22 and
00195 *                 A23 are empty if IB = KD. The upper triangle of A13
00196 *                 lies outside the band.
00197 *
00198                   I2 = MIN( KD-IB, N-I-IB+1 )
00199                   I3 = MIN( IB, N-I-KD+1 )
00200 *
00201                   IF( I2.GT.0 ) THEN
00202 *
00203 *                    Update A12
00204 *
00205                      CALL STRSM( 'Left', 'Upper', 'Transpose',
00206      \$                           'Non-unit', IB, I2, ONE, AB( KD+1, I ),
00207      \$                           LDAB-1, AB( KD+1-IB, I+IB ), LDAB-1 )
00208 *
00209 *                    Update A22
00210 *
00211                      CALL SSYRK( 'Upper', 'Transpose', I2, IB, -ONE,
00212      \$                           AB( KD+1-IB, I+IB ), LDAB-1, ONE,
00213      \$                           AB( KD+1, I+IB ), LDAB-1 )
00214                   END IF
00215 *
00216                   IF( I3.GT.0 ) THEN
00217 *
00218 *                    Copy the lower triangle of A13 into the work array.
00219 *
00220                      DO 40 JJ = 1, I3
00221                         DO 30 II = JJ, IB
00222                            WORK( II, JJ ) = AB( II-JJ+1, JJ+I+KD-1 )
00223    30                   CONTINUE
00224    40                CONTINUE
00225 *
00226 *                    Update A13 (in the work array).
00227 *
00228                      CALL STRSM( 'Left', 'Upper', 'Transpose',
00229      \$                           'Non-unit', IB, I3, ONE, AB( KD+1, I ),
00230      \$                           LDAB-1, WORK, LDWORK )
00231 *
00232 *                    Update A23
00233 *
00234                      IF( I2.GT.0 )
00235      \$                  CALL SGEMM( 'Transpose', 'No Transpose', I2, I3,
00236      \$                              IB, -ONE, AB( KD+1-IB, I+IB ),
00237      \$                              LDAB-1, WORK, LDWORK, ONE,
00238      \$                              AB( 1+IB, I+KD ), LDAB-1 )
00239 *
00240 *                    Update A33
00241 *
00242                      CALL SSYRK( 'Upper', 'Transpose', I3, IB, -ONE,
00243      \$                           WORK, LDWORK, ONE, AB( KD+1, I+KD ),
00244      \$                           LDAB-1 )
00245 *
00246 *                    Copy the lower triangle of A13 back into place.
00247 *
00248                      DO 60 JJ = 1, I3
00249                         DO 50 II = JJ, IB
00250                            AB( II-JJ+1, JJ+I+KD-1 ) = WORK( II, JJ )
00251    50                   CONTINUE
00252    60                CONTINUE
00253                   END IF
00254                END IF
00255    70       CONTINUE
00256          ELSE
00257 *
00258 *           Compute the Cholesky factorization of a symmetric band
00259 *           matrix, given the lower triangle of the matrix in band
00260 *           storage.
00261 *
00262 *           Zero the lower triangle of the work array.
00263 *
00264             DO 90 J = 1, NB
00265                DO 80 I = J + 1, NB
00266                   WORK( I, J ) = ZERO
00267    80          CONTINUE
00268    90       CONTINUE
00269 *
00270 *           Process the band matrix one diagonal block at a time.
00271 *
00272             DO 140 I = 1, N, NB
00273                IB = MIN( NB, N-I+1 )
00274 *
00275 *              Factorize the diagonal block
00276 *
00277                CALL SPOTF2( UPLO, IB, AB( 1, I ), LDAB-1, II )
00278                IF( II.NE.0 ) THEN
00279                   INFO = I + II - 1
00280                   GO TO 150
00281                END IF
00282                IF( I+IB.LE.N ) THEN
00283 *
00284 *                 Update the relevant part of the trailing submatrix.
00285 *                 If A11 denotes the diagonal block which has just been
00286 *                 factorized, then we need to update the remaining
00287 *                 blocks in the diagram:
00288 *
00289 *                    A11
00290 *                    A21   A22
00291 *                    A31   A32   A33
00292 *
00293 *                 The numbers of rows and columns in the partitioning
00294 *                 are IB, I2, I3 respectively. The blocks A21, A22 and
00295 *                 A32 are empty if IB = KD. The lower triangle of A31
00296 *                 lies outside the band.
00297 *
00298                   I2 = MIN( KD-IB, N-I-IB+1 )
00299                   I3 = MIN( IB, N-I-KD+1 )
00300 *
00301                   IF( I2.GT.0 ) THEN
00302 *
00303 *                    Update A21
00304 *
00305                      CALL STRSM( 'Right', 'Lower', 'Transpose',
00306      \$                           'Non-unit', I2, IB, ONE, AB( 1, I ),
00307      \$                           LDAB-1, AB( 1+IB, I ), LDAB-1 )
00308 *
00309 *                    Update A22
00310 *
00311                      CALL SSYRK( 'Lower', 'No Transpose', I2, IB, -ONE,
00312      \$                           AB( 1+IB, I ), LDAB-1, ONE,
00313      \$                           AB( 1, I+IB ), LDAB-1 )
00314                   END IF
00315 *
00316                   IF( I3.GT.0 ) THEN
00317 *
00318 *                    Copy the upper triangle of A31 into the work array.
00319 *
00320                      DO 110 JJ = 1, IB
00321                         DO 100 II = 1, MIN( JJ, I3 )
00322                            WORK( II, JJ ) = AB( KD+1-JJ+II, JJ+I-1 )
00323   100                   CONTINUE
00324   110                CONTINUE
00325 *
00326 *                    Update A31 (in the work array).
00327 *
00328                      CALL STRSM( 'Right', 'Lower', 'Transpose',
00329      \$                           'Non-unit', I3, IB, ONE, AB( 1, I ),
00330      \$                           LDAB-1, WORK, LDWORK )
00331 *
00332 *                    Update A32
00333 *
00334                      IF( I2.GT.0 )
00335      \$                  CALL SGEMM( 'No transpose', 'Transpose', I3, I2,
00336      \$                              IB, -ONE, WORK, LDWORK,
00337      \$                              AB( 1+IB, I ), LDAB-1, ONE,
00338      \$                              AB( 1+KD-IB, I+IB ), LDAB-1 )
00339 *
00340 *                    Update A33
00341 *
00342                      CALL SSYRK( 'Lower', 'No Transpose', I3, IB, -ONE,
00343      \$                           WORK, LDWORK, ONE, AB( 1, I+KD ),
00344      \$                           LDAB-1 )
00345 *
00346 *                    Copy the upper triangle of A31 back into place.
00347 *
00348                      DO 130 JJ = 1, IB
00349                         DO 120 II = 1, MIN( JJ, I3 )
00350                            AB( KD+1-JJ+II, JJ+I-1 ) = WORK( II, JJ )
00351   120                   CONTINUE
00352   130                CONTINUE
00353                   END IF
00354                END IF
00355   140       CONTINUE
00356          END IF
00357       END IF
00358       RETURN
00359 *
00360   150 CONTINUE
00361       RETURN
00362 *
00363 *     End of SPBTRF
00364 *
00365       END