LAPACK 3.3.0

dgbtf2.f

Go to the documentation of this file.
00001       SUBROUTINE DGBTF2( 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 *  DGBTF2 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 unblocked version of the algorithm, calling Level 2 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
00085 *  interchanges.
00086 *
00087 *  =====================================================================
00088 *
00089 *     .. Parameters ..
00090       DOUBLE PRECISION   ONE, ZERO
00091       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00092 *     ..
00093 *     .. Local Scalars ..
00094       INTEGER            I, J, JP, JU, KM, KV
00095 *     ..
00096 *     .. External Functions ..
00097       INTEGER            IDAMAX
00098       EXTERNAL           IDAMAX
00099 *     ..
00100 *     .. External Subroutines ..
00101       EXTERNAL           DGER, DSCAL, DSWAP, XERBLA
00102 *     ..
00103 *     .. Intrinsic Functions ..
00104       INTRINSIC          MAX, MIN
00105 *     ..
00106 *     .. Executable Statements ..
00107 *
00108 *     KV is the number of superdiagonals in the factor U, allowing for
00109 *     fill-in.
00110 *
00111       KV = KU + KL
00112 *
00113 *     Test the input parameters.
00114 *
00115       INFO = 0
00116       IF( M.LT.0 ) THEN
00117          INFO = -1
00118       ELSE IF( N.LT.0 ) THEN
00119          INFO = -2
00120       ELSE IF( KL.LT.0 ) THEN
00121          INFO = -3
00122       ELSE IF( KU.LT.0 ) THEN
00123          INFO = -4
00124       ELSE IF( LDAB.LT.KL+KV+1 ) THEN
00125          INFO = -6
00126       END IF
00127       IF( INFO.NE.0 ) THEN
00128          CALL XERBLA( 'DGBTF2', -INFO )
00129          RETURN
00130       END IF
00131 *
00132 *     Quick return if possible
00133 *
00134       IF( M.EQ.0 .OR. N.EQ.0 )
00135      $   RETURN
00136 *
00137 *     Gaussian elimination with partial pivoting
00138 *
00139 *     Set fill-in elements in columns KU+2 to KV to zero.
00140 *
00141       DO 20 J = KU + 2, MIN( KV, N )
00142          DO 10 I = KV - J + 2, KL
00143             AB( I, J ) = ZERO
00144    10    CONTINUE
00145    20 CONTINUE
00146 *
00147 *     JU is the index of the last column affected by the current stage
00148 *     of the factorization.
00149 *
00150       JU = 1
00151 *
00152       DO 40 J = 1, MIN( M, N )
00153 *
00154 *        Set fill-in elements in column J+KV to zero.
00155 *
00156          IF( J+KV.LE.N ) THEN
00157             DO 30 I = 1, KL
00158                AB( I, J+KV ) = ZERO
00159    30       CONTINUE
00160          END IF
00161 *
00162 *        Find pivot and test for singularity. KM is the number of
00163 *        subdiagonal elements in the current column.
00164 *
00165          KM = MIN( KL, M-J )
00166          JP = IDAMAX( KM+1, AB( KV+1, J ), 1 )
00167          IPIV( J ) = JP + J - 1
00168          IF( AB( KV+JP, J ).NE.ZERO ) THEN
00169             JU = MAX( JU, MIN( J+KU+JP-1, N ) )
00170 *
00171 *           Apply interchange to columns J to JU.
00172 *
00173             IF( JP.NE.1 )
00174      $         CALL DSWAP( JU-J+1, AB( KV+JP, J ), LDAB-1,
00175      $                     AB( KV+1, J ), LDAB-1 )
00176 *
00177             IF( KM.GT.0 ) THEN
00178 *
00179 *              Compute multipliers.
00180 *
00181                CALL DSCAL( KM, ONE / AB( KV+1, J ), AB( KV+2, J ), 1 )
00182 *
00183 *              Update trailing submatrix within the band.
00184 *
00185                IF( JU.GT.J )
00186      $            CALL DGER( KM, JU-J, -ONE, AB( KV+2, J ), 1,
00187      $                       AB( KV, J+1 ), LDAB-1, AB( KV+1, J+1 ),
00188      $                       LDAB-1 )
00189             END IF
00190          ELSE
00191 *
00192 *           If pivot is zero, set INFO to the index of the pivot
00193 *           unless a zero pivot has already been found.
00194 *
00195             IF( INFO.EQ.0 )
00196      $         INFO = J
00197          END IF
00198    40 CONTINUE
00199       RETURN
00200 *
00201 *     End of DGBTF2
00202 *
00203       END
 All Files Functions