ScaLAPACK  2.0.2 ScaLAPACK: Scalable Linear Algebra PACKage
ddbtrf.f
Go to the documentation of this file.
```00001       SUBROUTINE DDBTRF( M, N, KL, KU, AB, LDAB, INFO )
00002 *
00003 *  -- ScaLAPACK auxiliary routine (version 2.0) --
00004 *     Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver
00005 *
00006 *     Written by Andrew J. Cleary, University of Tennessee.
00007 *     August, 1996.
00008 *     Modified from DGBTRF:
00009 *  -- LAPACK routine (preliminary version) --
00010 *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
00011 *     Courant Institute, Argonne National Lab, and Rice University
00012 *     August 6, 1991
00013 *
00014 *     .. Scalar Arguments ..
00015       INTEGER            INFO, KL, KU, LDAB, M, N
00016 *     ..
00017 *     .. Array Arguments ..
00018       DOUBLE PRECISION   AB( LDAB, * )
00019 *     ..
00020 *
00021 *  Purpose
00022 *  =======
00023 *
00024 *  Ddbtrf computes an LU factorization of a real m-by-n band matrix A
00025 *  without using partial pivoting or row interchanges.
00026 *
00027 *  This is the blocked version of the algorithm, calling Level 3 BLAS.
00028 *
00029 *  Arguments
00030 *  =========
00031 *
00032 *  M       (input) INTEGER
00033 *          The number of rows of the matrix A.  M >= 0.
00034 *
00035 *  N       (input) INTEGER
00036 *          The number of columns of the matrix A.  N >= 0.
00037 *
00038 *  KL      (input) INTEGER
00039 *          The number of subdiagonals within the band of A.  KL >= 0.
00040 *
00041 *  KU      (input) INTEGER
00042 *          The number of superdiagonals within the band of A.  KU >= 0.
00043 *
00044 *  AB      (input/output) REAL array, dimension (LDAB,N)
00045 *          On entry, the matrix A in band storage, in rows KL+1 to
00046 *          2*KL+KU+1; rows 1 to KL of the array need not be set.
00047 *          The j-th column of A is stored in the j-th column of the
00048 *          array AB as follows:
00049 *          AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)
00050 *
00051 *          On exit, details of the factorization: U is stored as an
00052 *          upper triangular band matrix with KL+KU superdiagonals in
00053 *          rows 1 to KL+KU+1, and the multipliers used during the
00054 *          factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
00055 *          See below for further details.
00056 *
00057 *  LDAB    (input) INTEGER
00058 *          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
00059 *
00060 *  INFO    (output) INTEGER
00061 *          = 0: successful exit
00062 *          < 0: if INFO = -i, the i-th argument had an illegal value
00063 *          > 0: if INFO = +i, U(i,i) is exactly zero. The factorization
00064 *               has been completed, but the factor U is exactly
00065 *               singular, and division by zero will occur if it is used
00066 *               to solve a system of equations.
00067 *
00068 *  Further Details
00069 *  ===============
00070 *
00071 *  The band storage scheme is illustrated by the following example, when
00072 *  M = N = 6, KL = 2, KU = 1:
00073 *
00074 *  On entry:                       On exit:
00075 *
00076 *      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
00077 *     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
00078 *     a21  a32  a43  a54  a65   *      m21  m32  m43  m54  m65   *
00079 *     a31  a42  a53  a64   *    *      m31  m42  m53  m64   *    *
00080 *
00081 *  Array elements marked * are not used by the routine.
00082 *
00083 *  =====================================================================
00084 *
00085 *     .. Parameters ..
00086       DOUBLE PRECISION   ONE, ZERO
00087       PARAMETER          ( ONE = 1.0D+0 )
00088       PARAMETER          ( ZERO = 0.0D+0 )
00089       INTEGER            NBMAX, LDWORK
00090       PARAMETER          ( NBMAX = 64, LDWORK = NBMAX+1 )
00091 *     ..
00092 *     .. Local Scalars ..
00093       INTEGER            I, I2, I3, II, J, J2, J3, JB, JJ, JM, JP,
00094      \$                   JU, KM, KV, NB, NW
00095 *     ..
00096 *     .. Local Arrays ..
00097       DOUBLE PRECISION     WORK13( LDWORK, NBMAX ),
00098      \$                   WORK31( LDWORK, NBMAX )
00099 *     ..
00100 *     .. External Functions ..
00101       INTEGER            ILAENV, ISAMAX
00102       EXTERNAL           ILAENV, ISAMAX
00103 *     ..
00104 *     .. External Subroutines ..
00105       EXTERNAL           DCOPY, DDBTF2, DGEMM, DGER, DSCAL,
00106      \$                   DSWAP, DTRSM, XERBLA
00107 *     ..
00108 *     .. Intrinsic Functions ..
00109       INTRINSIC          MAX, MIN
00110 *     ..
00111 *     .. Executable Statements ..
00112 *
00113 *     KV is the number of superdiagonals in the factor U
00114 *
00115       KV = KU
00116 *
00117 *     Test the input parameters.
00118 *
00119       INFO = 0
00120       IF( M.LT.0 ) THEN
00121          INFO = -1
00122       ELSE IF( N.LT.0 ) THEN
00123          INFO = -2
00124       ELSE IF( KL.LT.0 ) THEN
00125          INFO = -3
00126       ELSE IF( KU.LT.0 ) THEN
00127          INFO = -4
00128       ELSE IF( LDAB.LT.MIN( MIN( KL+KV+1,M ),N ) ) THEN
00129          INFO = -6
00130       END IF
00131       IF( INFO.NE.0 ) THEN
00132          CALL XERBLA( 'DDBTRF', -INFO )
00133          RETURN
00134       END IF
00135 *
00136 *     Quick return if possible
00137 *
00138       IF( M.EQ.0 .OR. N.EQ.0 )
00139      \$     RETURN
00140 *
00141 *     Determine the block size for this environment
00142 *
00143       NB = ILAENV( 1, 'DDBTRF', ' ', M, N, KL, KU )
00144 *
00145 *     The block size must not exceed the limit set by the size of the
00146 *     local arrays WORK13 and WORK31.
00147 *
00148       NB = MIN( NB, NBMAX )
00149 *
00150       IF( NB.LE.1 .OR. NB.GT.KL ) THEN
00151 *
00152 *        Use unblocked code
00153 *
00154          CALL DDBTF2( M, N, KL, KU, AB, LDAB, INFO )
00155       ELSE
00156 *
00157 *        Use blocked code
00158 *
00159 *        Zero the superdiagonal elements of the work array WORK13
00160 *
00161          DO 20 J = 1, NB
00162             DO 10 I = 1, J - 1
00163                WORK13( I, J ) = ZERO
00164    10       CONTINUE
00165    20    CONTINUE
00166 *
00167 *        Zero the subdiagonal elements of the work array WORK31
00168 *
00169          DO 40 J = 1, NB
00170             DO 30 I = J + 1, NB
00171                WORK31( I, J ) = ZERO
00172    30       CONTINUE
00173    40    CONTINUE
00174 *
00175 *        JU is the index of the last column affected by the current
00176 *        stage of the factorization
00177 *
00178          JU = 1
00179 *
00180          DO 180 J = 1, MIN( M, N ), NB
00181             JB = MIN( NB, MIN( M, N )-J+1 )
00182 *
00183 *           The active part of the matrix is partitioned
00184 *
00185 *              A11   A12   A13
00186 *              A21   A22   A23
00187 *              A31   A32   A33
00188 *
00189 *           Here A11, A21 and A31 denote the current block of JB columns
00190 *           which is about to be factorized. The number of rows in the
00191 *           partitioning are JB, I2, I3 respectively, and the numbers
00192 *           of columns are JB, J2, J3. The superdiagonal elements of A13
00193 *           and the subdiagonal elements of A31 lie outside the band.
00194 *
00195             I2 = MIN( KL-JB, M-J-JB+1 )
00196             I3 = MIN( JB, M-J-KL+1 )
00197 *
00198 *           J2 and J3 are computed after JU has been updated.
00199 *
00200 *           Factorize the current block of JB columns
00201 *
00202             DO 80 JJ = J, J + JB - 1
00203 *
00204 *              Find pivot and test for singularity. KM is the number of
00205 *              subdiagonal elements in the current column.
00206 *
00207                KM = MIN( KL, M-JJ )
00208                JP = 1
00209                IF( AB( KV+JP, JJ ).NE.ZERO ) THEN
00210                   JU = MAX( JU, MIN( JJ+KU+JP-1, N ) )
00211 *
00212 *                 Compute multipliers
00213 *
00214                   CALL DSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, JJ ),
00215      \$                 1 )
00216 *
00217 *                 Update trailing submatrix within the band and within
00218 *                 the current block. JM is the index of the last column
00219 *                 which needs to be updated.
00220 *
00221                   JM = MIN( JU, J+JB-1 )
00222                   IF( JM.GT.JJ ) THEN
00223                      CALL DGER( KM, JM-JJ, -ONE, AB( KV+2, JJ ), 1,
00224      \$                          AB( KV, JJ+1 ), LDAB-1,
00225      \$                          AB( KV+1, JJ+1 ), LDAB-1 )
00226                   END IF
00227                END IF
00228 *
00229 *              Copy current column of A31 into the work array WORK31
00230 *
00231                NW = MIN( JJ-J+1, I3 )
00232                IF( NW.GT.0 )
00233      \$            CALL DCOPY( NW, AB( KV+KL+1-JJ+J, JJ ), 1,
00234      \$                        WORK31( 1, JJ-J+1 ), 1 )
00235    80       CONTINUE
00236             IF( J+JB.LE.N ) THEN
00237 *
00238 *              Apply the row interchanges to the other blocks.
00239 *
00240                J2 = MIN( JU-J+1, KV ) - JB
00241                J3 = MAX( 0, JU-J-KV+1 )
00242 *
00243 *              Update the relevant part of the trailing submatrix
00244 *
00245                IF( J2.GT.0 ) THEN
00246 *
00247 *                 Update A12
00248 *
00249                   CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit',
00250      \$                        JB, J2, ONE, AB( KV+1, J ), LDAB-1,
00251      \$                        AB( KV+1-JB, J+JB ), LDAB-1 )
00252 *
00253                   IF( I2.GT.0 ) THEN
00254 *
00255 *                    Update A22
00256 *
00257                      CALL DGEMM( 'No transpose', 'No transpose', I2, J2,
00258      \$                           JB, -ONE, AB( KV+1+JB, J ), LDAB-1,
00259      \$                           AB( KV+1-JB, J+JB ), LDAB-1, ONE,
00260      \$                           AB( KV+1, J+JB ), LDAB-1 )
00261                   END IF
00262 *
00263                   IF( I3.GT.0 ) THEN
00264 *
00265 *                    Update A32
00266 *
00267                      CALL DGEMM( 'No transpose', 'No transpose', I3, J2,
00268      \$                           JB, -ONE, WORK31, LDWORK,
00269      \$                           AB( KV+1-JB, J+JB ), LDAB-1, ONE,
00270      \$                           AB( KV+KL+1-JB, J+JB ), LDAB-1 )
00271                   END IF
00272                END IF
00273 *
00274                IF( J3.GT.0 ) THEN
00275 *
00276 *                 Copy the lower triangle of A13 into the work array
00277 *                 WORK13
00278 *
00279                   DO 130 JJ = 1, J3
00280                      DO 120 II = JJ, JB
00281                         WORK13( II, JJ ) = AB( II-JJ+1, JJ+J+KV-1 )
00282   120                CONTINUE
00283   130             CONTINUE
00284 *
00285 *                 Update A13 in the work array
00286 *
00287                   CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit',
00288      \$                        JB, J3, ONE, AB( KV+1, J ), LDAB-1,
00289      \$                        WORK13, LDWORK )
00290 *
00291                   IF( I2.GT.0 ) THEN
00292 *
00293 *                    Update A23
00294 *
00295                      CALL DGEMM( 'No transpose', 'No transpose', I2, J3,
00296      \$                           JB, -ONE, AB( KV+1+JB, J ), LDAB-1,
00297      \$                           WORK13, LDWORK, ONE, AB( 1+JB, J+KV ),
00298      \$                           LDAB-1 )
00299                   END IF
00300 *
00301                   IF( I3.GT.0 ) THEN
00302 *
00303 *                    Update A33
00304 *
00305                      CALL DGEMM( 'No transpose', 'No transpose', I3, J3,
00306      \$                         JB, -ONE, WORK31, LDWORK, WORK13,
00307      \$                         LDWORK, ONE, AB( 1+KL, J+KV ), LDAB-1 )
00308                   END IF
00309 *
00310 *                 Copy the lower triangle of A13 back into place
00311 *
00312                   DO 150 JJ = 1, J3
00313                      DO 140 II = JJ, JB
00314                         AB( II-JJ+1, JJ+J+KV-1 ) = WORK13( II, JJ )
00315   140                CONTINUE
00316   150             CONTINUE
00317                END IF
00318             ELSE
00319             END IF
00320 *
00321 *           copy the upper triangle of A31 back into place
00322 *
00323             DO 170 JJ = J + JB - 1, J, -1
00324 *
00325 *              Copy the current column of A31 back into place
00326 *
00327                NW = MIN( I3, JJ-J+1 )
00328                IF( NW.GT.0 )
00329      \$            CALL DCOPY( NW, WORK31( 1, JJ-J+1 ), 1,
00330      \$                        AB( KV+KL+1-JJ+J, JJ ), 1 )
00331   170       CONTINUE
00332   180    CONTINUE
00333       END IF
00334 *
00335       RETURN
00336 *
00337 *     End of DDBTRF
00338 *
00339       END
```