LAPACK 3.3.0

zgbtf2.f

Go to the documentation of this file.
00001       SUBROUTINE ZGBTF2( 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       COMPLEX*16         AB( LDAB, * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  ZGBTF2 computes an LU factorization of a complex m-by-n band matrix
00020 *  A 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) COMPLEX*16 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       COMPLEX*16         ONE, ZERO
00091       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
00092      $                   ZERO = ( 0.0D+0, 0.0D+0 ) )
00093 *     ..
00094 *     .. Local Scalars ..
00095       INTEGER            I, J, JP, JU, KM, KV
00096 *     ..
00097 *     .. External Functions ..
00098       INTEGER            IZAMAX
00099       EXTERNAL           IZAMAX
00100 *     ..
00101 *     .. External Subroutines ..
00102       EXTERNAL           XERBLA, ZGERU, ZSCAL, ZSWAP
00103 *     ..
00104 *     .. Intrinsic Functions ..
00105       INTRINSIC          MAX, MIN
00106 *     ..
00107 *     .. Executable Statements ..
00108 *
00109 *     KV is the number of superdiagonals in the factor U, allowing for
00110 *     fill-in.
00111 *
00112       KV = KU + KL
00113 *
00114 *     Test the input parameters.
00115 *
00116       INFO = 0
00117       IF( M.LT.0 ) THEN
00118          INFO = -1
00119       ELSE IF( N.LT.0 ) THEN
00120          INFO = -2
00121       ELSE IF( KL.LT.0 ) THEN
00122          INFO = -3
00123       ELSE IF( KU.LT.0 ) THEN
00124          INFO = -4
00125       ELSE IF( LDAB.LT.KL+KV+1 ) THEN
00126          INFO = -6
00127       END IF
00128       IF( INFO.NE.0 ) THEN
00129          CALL XERBLA( 'ZGBTF2', -INFO )
00130          RETURN
00131       END IF
00132 *
00133 *     Quick return if possible
00134 *
00135       IF( M.EQ.0 .OR. N.EQ.0 )
00136      $   RETURN
00137 *
00138 *     Gaussian elimination with partial pivoting
00139 *
00140 *     Set fill-in elements in columns KU+2 to KV to zero.
00141 *
00142       DO 20 J = KU + 2, MIN( KV, N )
00143          DO 10 I = KV - J + 2, KL
00144             AB( I, J ) = ZERO
00145    10    CONTINUE
00146    20 CONTINUE
00147 *
00148 *     JU is the index of the last column affected by the current stage
00149 *     of the factorization.
00150 *
00151       JU = 1
00152 *
00153       DO 40 J = 1, MIN( M, N )
00154 *
00155 *        Set fill-in elements in column J+KV to zero.
00156 *
00157          IF( J+KV.LE.N ) THEN
00158             DO 30 I = 1, KL
00159                AB( I, J+KV ) = ZERO
00160    30       CONTINUE
00161          END IF
00162 *
00163 *        Find pivot and test for singularity. KM is the number of
00164 *        subdiagonal elements in the current column.
00165 *
00166          KM = MIN( KL, M-J )
00167          JP = IZAMAX( KM+1, AB( KV+1, J ), 1 )
00168          IPIV( J ) = JP + J - 1
00169          IF( AB( KV+JP, J ).NE.ZERO ) THEN
00170             JU = MAX( JU, MIN( J+KU+JP-1, N ) )
00171 *
00172 *           Apply interchange to columns J to JU.
00173 *
00174             IF( JP.NE.1 )
00175      $         CALL ZSWAP( JU-J+1, AB( KV+JP, J ), LDAB-1,
00176      $                     AB( KV+1, J ), LDAB-1 )
00177             IF( KM.GT.0 ) THEN
00178 *
00179 *              Compute multipliers.
00180 *
00181                CALL ZSCAL( 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 ZGERU( 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 ZGBTF2
00202 *
00203       END
 All Files Functions