LAPACK: SRC/dgbtrf.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DGBTRF( M, N, KL, KU, AB, LDAB, IPIV, 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       INTEGER            INFO, KL, KU, LDAB, M, N
00010 *     ..
00011 *     .. Array Arguments ..
00012       INTEGER            IPIV( * )
00013       DOUBLE PRECISION   AB( LDAB, * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DGBTRF computes an LU factorization of a real m-by-n band matrix A
00020 *  using partial pivoting with row interchanges.
00021 *
00022 *  This is the blocked version of the algorithm, calling Level 3 BLAS.
00023 *
00024 *  Arguments
00025 *  =========
00026 *
00027 *  M       (input) INTEGER
00028 *          The number of rows of the matrix A.  M >= 0.
00029 *
00030 *  N       (input) INTEGER
00031 *          The number of columns of the matrix A.  N >= 0.
00032 *
00033 *  KL      (input) INTEGER
00034 *          The number of subdiagonals within the band of A.  KL >= 0.
00035 *
00036 *  KU      (input) INTEGER
00037 *          The number of superdiagonals within the band of A.  KU >= 0.
00038 *
00039 *  AB      (input/output) DOUBLE PRECISION array, dimension (LDAB,N)
00040 *          On entry, the matrix A in band storage, in rows KL+1 to
00041 *          2*KL+KU+1; rows 1 to KL of the array need not be set.
00042 *          The j-th column of A is stored in the j-th column of the
00043 *          array AB as follows:
00044 *          AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl)
00045 *
00046 *          On exit, details of the factorization: U is stored as an
00047 *          upper triangular band matrix with KL+KU superdiagonals in
00048 *          rows 1 to KL+KU+1, and the multipliers used during the
00049 *          factorization are stored in rows KL+KU+2 to 2*KL+KU+1.
00050 *          See below for further details.
00051 *
00052 *  LDAB    (input) INTEGER
00053 *          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
00054 *
00055 *  IPIV    (output) INTEGER array, dimension (min(M,N))
00056 *          The pivot indices; for 1 <= i <= min(M,N), row i of the
00057 *          matrix was interchanged with row IPIV(i).
00058 *
00059 *  INFO    (output) INTEGER
00060 *          = 0: successful exit
00061 *          < 0: if INFO = -i, the i-th argument had an illegal value
00062 *          > 0: if INFO = +i, U(i,i) is exactly zero. The factorization
00063 *               has been completed, but the factor U is exactly
00064 *               singular, and division by zero will occur if it is used
00065 *               to solve a system of equations.
00066 *
00067 *  Further Details
00068 *  ===============
00069 *
00070 *  The band storage scheme is illustrated by the following example, when
00071 *  M = N = 6, KL = 2, KU = 1:
00072 *
00073 *  On entry:                       On exit:
00074 *
00075 *      *    *    *    +    +    +       *    *    *   u14  u25  u36
00076 *      *    *    +    +    +    +       *    *   u13  u24  u35  u46
00077 *      *   a12  a23  a34  a45  a56      *   u12  u23  u34  u45  u56
00078 *     a11  a22  a33  a44  a55  a66     u11  u22  u33  u44  u55  u66
00079 *     a21  a32  a43  a54  a65   *      m21  m32  m43  m54  m65   *
00080 *     a31  a42  a53  a64   *    *      m31  m42  m53  m64   *    *
00081 *
00082 *  Array elements marked * are not used by the routine; elements marked
00083 *  + need not be set on entry, but are required by the routine to store
00084 *  elements of U because of fill-in resulting from the row interchanges.
00085 *
00086 *  =====================================================================
00087 *
00088 *     .. Parameters ..
00089       DOUBLE PRECISION   ONE, ZERO
00090       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00091       INTEGER            NBMAX, LDWORK
00092       PARAMETER          ( NBMAX = 64, LDWORK = NBMAX+1 )
00093 *     ..
00094 *     .. Local Scalars ..
00095       INTEGER            I, I2, I3, II, IP, J, J2, J3, JB, JJ, JM, JP,
00096      $                   JU, K2, KM, KV, NB, NW
00097       DOUBLE PRECISION   TEMP
00098 *     ..
00099 *     .. Local Arrays ..
00100       DOUBLE PRECISION   WORK13( LDWORK, NBMAX ),
00101      $                   WORK31( LDWORK, NBMAX )
00102 *     ..
00103 *     .. External Functions ..
00104       INTEGER            IDAMAX, ILAENV
00105       EXTERNAL           IDAMAX, ILAENV
00106 *     ..
00107 *     .. External Subroutines ..
00108       EXTERNAL           DCOPY, DGBTF2, DGEMM, DGER, DLASWP, DSCAL,
00109      $                   DSWAP, DTRSM, XERBLA
00110 *     ..
00111 *     .. Intrinsic Functions ..
00112       INTRINSIC          MAX, MIN
00113 *     ..
00114 *     .. Executable Statements ..
00115 *
00116 *     KV is the number of superdiagonals in the factor U, allowing for
00117 *     fill-in
00118 *
00119       KV = KU + KL
00120 *
00121 *     Test the input parameters.
00122 *
00123       INFO = 0
00124       IF( M.LT.0 ) THEN
00125          INFO = -1
00126       ELSE IF( N.LT.0 ) THEN
00127          INFO = -2
00128       ELSE IF( KL.LT.0 ) THEN
00129          INFO = -3
00130       ELSE IF( KU.LT.0 ) THEN
00131          INFO = -4
00132       ELSE IF( LDAB.LT.KL+KV+1 ) THEN
00133          INFO = -6
00134       END IF
00135       IF( INFO.NE.0 ) THEN
00136          CALL XERBLA( 'DGBTRF', -INFO )
00137          RETURN
00138       END IF
00139 *
00140 *     Quick return if possible
00141 *
00142       IF( M.EQ.0 .OR. N.EQ.0 )
00143      $   RETURN
00144 *
00145 *     Determine the block size for this environment
00146 *
00147       NB = ILAENV( 1, 'DGBTRF', ' ', M, N, KL, KU )
00148 *
00149 *     The block size must not exceed the limit set by the size of the
00150 *     local arrays WORK13 and WORK31.
00151 *
00152       NB = MIN( NB, NBMAX )
00153 *
00154       IF( NB.LE.1 .OR. NB.GT.KL ) THEN
00155 *
00156 *        Use unblocked code
00157 *
00158          CALL DGBTF2( M, N, KL, KU, AB, LDAB, IPIV, INFO )
00159       ELSE
00160 *
00161 *        Use blocked code
00162 *
00163 *        Zero the superdiagonal elements of the work array WORK13
00164 *
00165          DO 20 J = 1, NB
00166             DO 10 I = 1, J - 1
00167                WORK13( I, J ) = ZERO
00168    10       CONTINUE
00169    20    CONTINUE
00170 *
00171 *        Zero the subdiagonal elements of the work array WORK31
00172 *
00173          DO 40 J = 1, NB
00174             DO 30 I = J + 1, NB
00175                WORK31( I, J ) = ZERO
00176    30       CONTINUE
00177    40    CONTINUE
00178 *
00179 *        Gaussian elimination with partial pivoting
00180 *
00181 *        Set fill-in elements in columns KU+2 to KV to zero
00182 *
00183          DO 60 J = KU + 2, MIN( KV, N )
00184             DO 50 I = KV - J + 2, KL
00185                AB( I, J ) = ZERO
00186    50       CONTINUE
00187    60    CONTINUE
00188 *
00189 *        JU is the index of the last column affected by the current
00190 *        stage of the factorization
00191 *
00192          JU = 1
00193 *
00194          DO 180 J = 1, MIN( M, N ), NB
00195             JB = MIN( NB, MIN( M, N )-J+1 )
00196 *
00197 *           The active part of the matrix is partitioned
00198 *
00199 *              A11   A12   A13
00200 *              A21   A22   A23
00201 *              A31   A32   A33
00202 *
00203 *           Here A11, A21 and A31 denote the current block of JB columns
00204 *           which is about to be factorized. The number of rows in the
00205 *           partitioning are JB, I2, I3 respectively, and the numbers
00206 *           of columns are JB, J2, J3. The superdiagonal elements of A13
00207 *           and the subdiagonal elements of A31 lie outside the band.
00208 *
00209             I2 = MIN( KL-JB, M-J-JB+1 )
00210             I3 = MIN( JB, M-J-KL+1 )
00211 *
00212 *           J2 and J3 are computed after JU has been updated.
00213 *
00214 *           Factorize the current block of JB columns
00215 *
00216             DO 80 JJ = J, J + JB - 1
00217 *
00218 *              Set fill-in elements in column JJ+KV to zero
00219 *
00220                IF( JJ+KV.LE.N ) THEN
00221                   DO 70 I = 1, KL
00222                      AB( I, JJ+KV ) = ZERO
00223    70             CONTINUE
00224                END IF
00225 *
00226 *              Find pivot and test for singularity. KM is the number of
00227 *              subdiagonal elements in the current column.
00228 *
00229                KM = MIN( KL, M-JJ )
00230                JP = IDAMAX( KM+1, AB( KV+1, JJ ), 1 )
00231                IPIV( JJ ) = JP + JJ - J
00232                IF( AB( KV+JP, JJ ).NE.ZERO ) THEN
00233                   JU = MAX( JU, MIN( JJ+KU+JP-1, N ) )
00234                   IF( JP.NE.1 ) THEN
00235 *
00236 *                    Apply interchange to columns J to J+JB-1
00237 *
00238                      IF( JP+JJ-1.LT.J+KL ) THEN
00239 *
00240                         CALL DSWAP( JB, AB( KV+1+JJ-J, J ), LDAB-1,
00241      $                              AB( KV+JP+JJ-J, J ), LDAB-1 )
00242                      ELSE
00243 *
00244 *                       The interchange affects columns J to JJ-1 of A31
00245 *                       which are stored in the work array WORK31
00246 *
00247                         CALL DSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1,
00248      $                              WORK31( JP+JJ-J-KL, 1 ), LDWORK )
00249                         CALL DSWAP( J+JB-JJ, AB( KV+1, JJ ), LDAB-1,
00250      $                              AB( KV+JP, JJ ), LDAB-1 )
00251                      END IF
00252                   END IF
00253 *
00254 *                 Compute multipliers
00255 *
00256                   CALL DSCAL( KM, ONE / AB( KV+1, JJ ), AB( KV+2, JJ ),
00257      $                        1 )
00258 *
00259 *                 Update trailing submatrix within the band and within
00260 *                 the current block. JM is the index of the last column
00261 *                 which needs to be updated.
00262 *
00263                   JM = MIN( JU, J+JB-1 )
00264                   IF( JM.GT.JJ )
00265      $               CALL DGER( KM, JM-JJ, -ONE, AB( KV+2, JJ ), 1,
00266      $                          AB( KV, JJ+1 ), LDAB-1,
00267      $                          AB( KV+1, JJ+1 ), LDAB-1 )
00268                ELSE
00269 *
00270 *                 If pivot is zero, set INFO to the index of the pivot
00271 *                 unless a zero pivot has already been found.
00272 *
00273                   IF( INFO.EQ.0 )
00274      $               INFO = JJ
00275                END IF
00276 *
00277 *              Copy current column of A31 into the work array WORK31
00278 *
00279                NW = MIN( JJ-J+1, I3 )
00280                IF( NW.GT.0 )
00281      $            CALL DCOPY( NW, AB( KV+KL+1-JJ+J, JJ ), 1,
00282      $                        WORK31( 1, JJ-J+1 ), 1 )
00283    80       CONTINUE
00284             IF( J+JB.LE.N ) THEN
00285 *
00286 *              Apply the row interchanges to the other blocks.
00287 *
00288                J2 = MIN( JU-J+1, KV ) - JB
00289                J3 = MAX( 0, JU-J-KV+1 )
00290 *
00291 *              Use DLASWP to apply the row interchanges to A12, A22, and
00292 *              A32.
00293 *
00294                CALL DLASWP( J2, AB( KV+1-JB, J+JB ), LDAB-1, 1, JB,
00295      $                      IPIV( J ), 1 )
00296 *
00297 *              Adjust the pivot indices.
00298 *
00299                DO 90 I = J, J + JB - 1
00300                   IPIV( I ) = IPIV( I ) + J - 1
00301    90          CONTINUE
00302 *
00303 *              Apply the row interchanges to A13, A23, and A33
00304 *              columnwise.
00305 *
00306                K2 = J - 1 + JB + J2
00307                DO 110 I = 1, J3
00308                   JJ = K2 + I
00309                   DO 100 II = J + I - 1, J + JB - 1
00310                      IP = IPIV( II )
00311                      IF( IP.NE.II ) THEN
00312                         TEMP = AB( KV+1+II-JJ, JJ )
00313                         AB( KV+1+II-JJ, JJ ) = AB( KV+1+IP-JJ, JJ )
00314                         AB( KV+1+IP-JJ, JJ ) = TEMP
00315                      END IF
00316   100             CONTINUE
00317   110          CONTINUE
00318 *
00319 *              Update the relevant part of the trailing submatrix
00320 *
00321                IF( J2.GT.0 ) THEN
00322 *
00323 *                 Update A12
00324 *
00325                   CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit',
00326      $                        JB, J2, ONE, AB( KV+1, J ), LDAB-1,
00327      $                        AB( KV+1-JB, J+JB ), LDAB-1 )
00328 *
00329                   IF( I2.GT.0 ) THEN
00330 *
00331 *                    Update A22
00332 *
00333                      CALL DGEMM( 'No transpose', 'No transpose', I2, J2,
00334      $                           JB, -ONE, AB( KV+1+JB, J ), LDAB-1,
00335      $                           AB( KV+1-JB, J+JB ), LDAB-1, ONE,
00336      $                           AB( KV+1, J+JB ), LDAB-1 )
00337                   END IF
00338 *
00339                   IF( I3.GT.0 ) THEN
00340 *
00341 *                    Update A32
00342 *
00343                      CALL DGEMM( 'No transpose', 'No transpose', I3, J2,
00344      $                           JB, -ONE, WORK31, LDWORK,
00345      $                           AB( KV+1-JB, J+JB ), LDAB-1, ONE,
00346      $                           AB( KV+KL+1-JB, J+JB ), LDAB-1 )
00347                   END IF
00348                END IF
00349 *
00350                IF( J3.GT.0 ) THEN
00351 *
00352 *                 Copy the lower triangle of A13 into the work array
00353 *                 WORK13
00354 *
00355                   DO 130 JJ = 1, J3
00356                      DO 120 II = JJ, JB
00357                         WORK13( II, JJ ) = AB( II-JJ+1, JJ+J+KV-1 )
00358   120                CONTINUE
00359   130             CONTINUE
00360 *
00361 *                 Update A13 in the work array
00362 *
00363                   CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Unit',
00364      $                        JB, J3, ONE, AB( KV+1, J ), LDAB-1,
00365      $                        WORK13, LDWORK )
00366 *
00367                   IF( I2.GT.0 ) THEN
00368 *
00369 *                    Update A23
00370 *
00371                      CALL DGEMM( 'No transpose', 'No transpose', I2, J3,
00372      $                           JB, -ONE, AB( KV+1+JB, J ), LDAB-1,
00373      $                           WORK13, LDWORK, ONE, AB( 1+JB, J+KV ),
00374      $                           LDAB-1 )
00375                   END IF
00376 *
00377                   IF( I3.GT.0 ) THEN
00378 *
00379 *                    Update A33
00380 *
00381                      CALL DGEMM( 'No transpose', 'No transpose', I3, J3,
00382      $                           JB, -ONE, WORK31, LDWORK, WORK13,
00383      $                           LDWORK, ONE, AB( 1+KL, J+KV ), LDAB-1 )
00384                   END IF
00385 *
00386 *                 Copy the lower triangle of A13 back into place
00387 *
00388                   DO 150 JJ = 1, J3
00389                      DO 140 II = JJ, JB
00390                         AB( II-JJ+1, JJ+J+KV-1 ) = WORK13( II, JJ )
00391   140                CONTINUE
00392   150             CONTINUE
00393                END IF
00394             ELSE
00395 *
00396 *              Adjust the pivot indices.
00397 *
00398                DO 160 I = J, J + JB - 1
00399                   IPIV( I ) = IPIV( I ) + J - 1
00400   160          CONTINUE
00401             END IF
00402 *
00403 *           Partially undo the interchanges in the current block to
00404 *           restore the upper triangular form of A31 and copy the upper
00405 *           triangle of A31 back into place
00406 *
00407             DO 170 JJ = J + JB - 1, J, -1
00408                JP = IPIV( JJ ) - JJ + 1
00409                IF( JP.NE.1 ) THEN
00410 *
00411 *                 Apply interchange to columns J to JJ-1
00412 *
00413                   IF( JP+JJ-1.LT.J+KL ) THEN
00414 *
00415 *                    The interchange does not affect A31
00416 *
00417                      CALL DSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1,
00418      $                           AB( KV+JP+JJ-J, J ), LDAB-1 )
00419                   ELSE
00420 *
00421 *                    The interchange does affect A31
00422 *
00423                      CALL DSWAP( JJ-J, AB( KV+1+JJ-J, J ), LDAB-1,
00424      $                           WORK31( JP+JJ-J-KL, 1 ), LDWORK )
00425                   END IF
00426                END IF
00427 *
00428 *              Copy the current column of A31 back into place
00429 *
00430                NW = MIN( I3, JJ-J+1 )
00431                IF( NW.GT.0 )
00432      $            CALL DCOPY( NW, WORK31( 1, JJ-J+1 ), 1,
00433      $                        AB( KV+KL+1-JJ+J, JJ ), 1 )
00434   170       CONTINUE
00435   180    CONTINUE
00436       END IF
00437 *
00438       RETURN
00439 *
00440 *     End of DGBTRF
00441 *
00442       END