LAPACK: SRC/dgbtrs.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DGBTRS( TRANS, N, KL, KU, NRHS, AB, LDAB, IPIV, B, LDB,
00002      $                   INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          TRANS
00011       INTEGER            INFO, KL, KU, LDAB, LDB, N, NRHS
00012 *     ..
00013 *     .. Array Arguments ..
00014       INTEGER            IPIV( * )
00015       DOUBLE PRECISION   AB( LDAB, * ), B( LDB, * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  DGBTRS solves a system of linear equations
00022 *     A * X = B  or  A' * X = B
00023 *  with a general band matrix A using the LU factorization computed
00024 *  by DGBTRF.
00025 *
00026 *  Arguments
00027 *  =========
00028 *
00029 *  TRANS   (input) CHARACTER*1
00030 *          Specifies the form of the system of equations.
00031 *          = 'N':  A * X = B  (No transpose)
00032 *          = 'T':  A'* X = B  (Transpose)
00033 *          = 'C':  A'* X = B  (Conjugate transpose = Transpose)
00034 *
00035 *  N       (input) INTEGER
00036 *          The order 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 *  NRHS    (input) INTEGER
00045 *          The number of right hand sides, i.e., the number of columns
00046 *          of the matrix B.  NRHS >= 0.
00047 *
00048 *  AB      (input) DOUBLE PRECISION array, dimension (LDAB,N)
00049 *          Details of the LU factorization of the band matrix A, as
00050 *          computed by DGBTRF.  U is stored as an upper triangular band
00051 *          matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, and
00052 *          the multipliers used during the factorization are stored in
00053 *          rows KL+KU+2 to 2*KL+KU+1.
00054 *
00055 *  LDAB    (input) INTEGER
00056 *          The leading dimension of the array AB.  LDAB >= 2*KL+KU+1.
00057 *
00058 *  IPIV    (input) INTEGER array, dimension (N)
00059 *          The pivot indices; for 1 <= i <= N, row i of the matrix was
00060 *          interchanged with row IPIV(i).
00061 *
00062 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
00063 *          On entry, the right hand side matrix B.
00064 *          On exit, the solution matrix X.
00065 *
00066 *  LDB     (input) INTEGER
00067 *          The leading dimension of the array B.  LDB >= max(1,N).
00068 *
00069 *  INFO    (output) INTEGER
00070 *          = 0:  successful exit
00071 *          < 0: if INFO = -i, the i-th argument had an illegal value
00072 *
00073 *  =====================================================================
00074 *
00075 *     .. Parameters ..
00076       DOUBLE PRECISION   ONE
00077       PARAMETER          ( ONE = 1.0D+0 )
00078 *     ..
00079 *     .. Local Scalars ..
00080       LOGICAL            LNOTI, NOTRAN
00081       INTEGER            I, J, KD, L, LM
00082 *     ..
00083 *     .. External Functions ..
00084       LOGICAL            LSAME
00085       EXTERNAL           LSAME
00086 *     ..
00087 *     .. External Subroutines ..
00088       EXTERNAL           DGEMV, DGER, DSWAP, DTBSV, XERBLA
00089 *     ..
00090 *     .. Intrinsic Functions ..
00091       INTRINSIC          MAX, MIN
00092 *     ..
00093 *     .. Executable Statements ..
00094 *
00095 *     Test the input parameters.
00096 *
00097       INFO = 0
00098       NOTRAN = LSAME( TRANS, 'N' )
00099       IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
00100      $    LSAME( TRANS, 'C' ) ) THEN
00101          INFO = -1
00102       ELSE IF( N.LT.0 ) THEN
00103          INFO = -2
00104       ELSE IF( KL.LT.0 ) THEN
00105          INFO = -3
00106       ELSE IF( KU.LT.0 ) THEN
00107          INFO = -4
00108       ELSE IF( NRHS.LT.0 ) THEN
00109          INFO = -5
00110       ELSE IF( LDAB.LT.( 2*KL+KU+1 ) ) THEN
00111          INFO = -7
00112       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00113          INFO = -10
00114       END IF
00115       IF( INFO.NE.0 ) THEN
00116          CALL XERBLA( 'DGBTRS', -INFO )
00117          RETURN
00118       END IF
00119 *
00120 *     Quick return if possible
00121 *
00122       IF( N.EQ.0 .OR. NRHS.EQ.0 )
00123      $   RETURN
00124 *
00125       KD = KU + KL + 1
00126       LNOTI = KL.GT.0
00127 *
00128       IF( NOTRAN ) THEN
00129 *
00130 *        Solve  A*X = B.
00131 *
00132 *        Solve L*X = B, overwriting B with X.
00133 *
00134 *        L is represented as a product of permutations and unit lower
00135 *        triangular matrices L = P(1) * L(1) * ... * P(n-1) * L(n-1),
00136 *        where each transformation L(i) is a rank-one modification of
00137 *        the identity matrix.
00138 *
00139          IF( LNOTI ) THEN
00140             DO 10 J = 1, N - 1
00141                LM = MIN( KL, N-J )
00142                L = IPIV( J )
00143                IF( L.NE.J )
00144      $            CALL DSWAP( NRHS, B( L, 1 ), LDB, B( J, 1 ), LDB )
00145                CALL DGER( LM, NRHS, -ONE, AB( KD+1, J ), 1, B( J, 1 ),
00146      $                    LDB, B( J+1, 1 ), LDB )
00147    10       CONTINUE
00148          END IF
00149 *
00150          DO 20 I = 1, NRHS
00151 *
00152 *           Solve U*X = B, overwriting B with X.
00153 *
00154             CALL DTBSV( 'Upper', 'No transpose', 'Non-unit', N, KL+KU,
00155      $                  AB, LDAB, B( 1, I ), 1 )
00156    20    CONTINUE
00157 *
00158       ELSE
00159 *
00160 *        Solve A'*X = B.
00161 *
00162          DO 30 I = 1, NRHS
00163 *
00164 *           Solve U'*X = B, overwriting B with X.
00165 *
00166             CALL DTBSV( 'Upper', 'Transpose', 'Non-unit', N, KL+KU, AB,
00167      $                  LDAB, B( 1, I ), 1 )
00168    30    CONTINUE
00169 *
00170 *        Solve L'*X = B, overwriting B with X.
00171 *
00172          IF( LNOTI ) THEN
00173             DO 40 J = N - 1, 1, -1
00174                LM = MIN( KL, N-J )
00175                CALL DGEMV( 'Transpose', LM, NRHS, -ONE, B( J+1, 1 ),
00176      $                     LDB, AB( KD+1, J ), 1, ONE, B( J, 1 ), LDB )
00177                L = IPIV( J )
00178                IF( L.NE.J )
00179      $            CALL DSWAP( NRHS, B( L, 1 ), LDB, B( J, 1 ), LDB )
00180    40       CONTINUE
00181          END IF
00182       END IF
00183       RETURN
00184 *
00185 *     End of DGBTRS
00186 *
00187       END