00001 /* dgbtf2.f -- translated by f2c (version 20061008). 00002 You must link the resulting object file with libf2c: 00003 on Microsoft Windows system, link with libf2c.lib; 00004 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 00005 or, if you install libf2c.a in a standard place, with -lf2c -lm 00006 -- in that order, at the end of the command line, as in 00007 cc *.o -lf2c -lm 00008 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 00009 00010 http://www.netlib.org/f2c/libf2c.zip 00011 */ 00012 00013 #include "f2c.h" 00014 #include "blaswrap.h" 00015 00016 /* Table of constant values */ 00017 00018 static integer c__1 = 1; 00019 static doublereal c_b9 = -1.; 00020 00021 /* Subroutine */ int dgbtf2_(integer *m, integer *n, integer *kl, integer *ku, 00022 doublereal *ab, integer *ldab, integer *ipiv, integer *info) 00023 { 00024 /* System generated locals */ 00025 integer ab_dim1, ab_offset, i__1, i__2, i__3, i__4; 00026 doublereal d__1; 00027 00028 /* Local variables */ 00029 integer i__, j, km, jp, ju, kv; 00030 extern /* Subroutine */ int dger_(integer *, integer *, doublereal *, 00031 doublereal *, integer *, doublereal *, integer *, doublereal *, 00032 integer *), dscal_(integer *, doublereal *, doublereal *, integer 00033 *), dswap_(integer *, doublereal *, integer *, doublereal *, 00034 integer *); 00035 extern integer idamax_(integer *, doublereal *, integer *); 00036 extern /* Subroutine */ int xerbla_(char *, integer *); 00037 00038 00039 /* -- LAPACK routine (version 3.2) -- */ 00040 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00041 /* November 2006 */ 00042 00043 /* .. Scalar Arguments .. */ 00044 /* .. */ 00045 /* .. Array Arguments .. */ 00046 /* .. */ 00047 00048 /* Purpose */ 00049 /* ======= */ 00050 00051 /* DGBTF2 computes an LU factorization of a real m-by-n band matrix A */ 00052 /* using partial pivoting with row interchanges. */ 00053 00054 /* This is the unblocked version of the algorithm, calling Level 2 BLAS. */ 00055 00056 /* Arguments */ 00057 /* ========= */ 00058 00059 /* M (input) INTEGER */ 00060 /* The number of rows of the matrix A. M >= 0. */ 00061 00062 /* N (input) INTEGER */ 00063 /* The number of columns of the matrix A. N >= 0. */ 00064 00065 /* KL (input) INTEGER */ 00066 /* The number of subdiagonals within the band of A. KL >= 0. */ 00067 00068 /* KU (input) INTEGER */ 00069 /* The number of superdiagonals within the band of A. KU >= 0. */ 00070 00071 /* AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N) */ 00072 /* On entry, the matrix A in band storage, in rows KL+1 to */ 00073 /* 2*KL+KU+1; rows 1 to KL of the array need not be set. */ 00074 /* The j-th column of A is stored in the j-th column of the */ 00075 /* array AB as follows: */ 00076 /* AB(kl+ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl) */ 00077 00078 /* On exit, details of the factorization: U is stored as an */ 00079 /* upper triangular band matrix with KL+KU superdiagonals in */ 00080 /* rows 1 to KL+KU+1, and the multipliers used during the */ 00081 /* factorization are stored in rows KL+KU+2 to 2*KL+KU+1. */ 00082 /* See below for further details. */ 00083 00084 /* LDAB (input) INTEGER */ 00085 /* The leading dimension of the array AB. LDAB >= 2*KL+KU+1. */ 00086 00087 /* IPIV (output) INTEGER array, dimension (min(M,N)) */ 00088 /* The pivot indices; for 1 <= i <= min(M,N), row i of the */ 00089 /* matrix was interchanged with row IPIV(i). */ 00090 00091 /* INFO (output) INTEGER */ 00092 /* = 0: successful exit */ 00093 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00094 /* > 0: if INFO = +i, U(i,i) is exactly zero. The factorization */ 00095 /* has been completed, but the factor U is exactly */ 00096 /* singular, and division by zero will occur if it is used */ 00097 /* to solve a system of equations. */ 00098 00099 /* Further Details */ 00100 /* =============== */ 00101 00102 /* The band storage scheme is illustrated by the following example, when */ 00103 /* M = N = 6, KL = 2, KU = 1: */ 00104 00105 /* On entry: On exit: */ 00106 00107 /* * * * + + + * * * u14 u25 u36 */ 00108 /* * * + + + + * * u13 u24 u35 u46 */ 00109 /* * a12 a23 a34 a45 a56 * u12 u23 u34 u45 u56 */ 00110 /* a11 a22 a33 a44 a55 a66 u11 u22 u33 u44 u55 u66 */ 00111 /* a21 a32 a43 a54 a65 * m21 m32 m43 m54 m65 * */ 00112 /* a31 a42 a53 a64 * * m31 m42 m53 m64 * * */ 00113 00114 /* Array elements marked * are not used by the routine; elements marked */ 00115 /* + need not be set on entry, but are required by the routine to store */ 00116 /* elements of U, because of fill-in resulting from the row */ 00117 /* interchanges. */ 00118 00119 /* ===================================================================== */ 00120 00121 /* .. Parameters .. */ 00122 /* .. */ 00123 /* .. Local Scalars .. */ 00124 /* .. */ 00125 /* .. External Functions .. */ 00126 /* .. */ 00127 /* .. External Subroutines .. */ 00128 /* .. */ 00129 /* .. Intrinsic Functions .. */ 00130 /* .. */ 00131 /* .. Executable Statements .. */ 00132 00133 /* KV is the number of superdiagonals in the factor U, allowing for */ 00134 /* fill-in. */ 00135 00136 /* Parameter adjustments */ 00137 ab_dim1 = *ldab; 00138 ab_offset = 1 + ab_dim1; 00139 ab -= ab_offset; 00140 --ipiv; 00141 00142 /* Function Body */ 00143 kv = *ku + *kl; 00144 00145 /* Test the input parameters. */ 00146 00147 *info = 0; 00148 if (*m < 0) { 00149 *info = -1; 00150 } else if (*n < 0) { 00151 *info = -2; 00152 } else if (*kl < 0) { 00153 *info = -3; 00154 } else if (*ku < 0) { 00155 *info = -4; 00156 } else if (*ldab < *kl + kv + 1) { 00157 *info = -6; 00158 } 00159 if (*info != 0) { 00160 i__1 = -(*info); 00161 xerbla_("DGBTF2", &i__1); 00162 return 0; 00163 } 00164 00165 /* Quick return if possible */ 00166 00167 if (*m == 0 || *n == 0) { 00168 return 0; 00169 } 00170 00171 /* Gaussian elimination with partial pivoting */ 00172 00173 /* Set fill-in elements in columns KU+2 to KV to zero. */ 00174 00175 i__1 = min(kv,*n); 00176 for (j = *ku + 2; j <= i__1; ++j) { 00177 i__2 = *kl; 00178 for (i__ = kv - j + 2; i__ <= i__2; ++i__) { 00179 ab[i__ + j * ab_dim1] = 0.; 00180 /* L10: */ 00181 } 00182 /* L20: */ 00183 } 00184 00185 /* JU is the index of the last column affected by the current stage */ 00186 /* of the factorization. */ 00187 00188 ju = 1; 00189 00190 i__1 = min(*m,*n); 00191 for (j = 1; j <= i__1; ++j) { 00192 00193 /* Set fill-in elements in column J+KV to zero. */ 00194 00195 if (j + kv <= *n) { 00196 i__2 = *kl; 00197 for (i__ = 1; i__ <= i__2; ++i__) { 00198 ab[i__ + (j + kv) * ab_dim1] = 0.; 00199 /* L30: */ 00200 } 00201 } 00202 00203 /* Find pivot and test for singularity. KM is the number of */ 00204 /* subdiagonal elements in the current column. */ 00205 00206 /* Computing MIN */ 00207 i__2 = *kl, i__3 = *m - j; 00208 km = min(i__2,i__3); 00209 i__2 = km + 1; 00210 jp = idamax_(&i__2, &ab[kv + 1 + j * ab_dim1], &c__1); 00211 ipiv[j] = jp + j - 1; 00212 if (ab[kv + jp + j * ab_dim1] != 0.) { 00213 /* Computing MAX */ 00214 /* Computing MIN */ 00215 i__4 = j + *ku + jp - 1; 00216 i__2 = ju, i__3 = min(i__4,*n); 00217 ju = max(i__2,i__3); 00218 00219 /* Apply interchange to columns J to JU. */ 00220 00221 if (jp != 1) { 00222 i__2 = ju - j + 1; 00223 i__3 = *ldab - 1; 00224 i__4 = *ldab - 1; 00225 dswap_(&i__2, &ab[kv + jp + j * ab_dim1], &i__3, &ab[kv + 1 + 00226 j * ab_dim1], &i__4); 00227 } 00228 00229 if (km > 0) { 00230 00231 /* Compute multipliers. */ 00232 00233 d__1 = 1. / ab[kv + 1 + j * ab_dim1]; 00234 dscal_(&km, &d__1, &ab[kv + 2 + j * ab_dim1], &c__1); 00235 00236 /* Update trailing submatrix within the band. */ 00237 00238 if (ju > j) { 00239 i__2 = ju - j; 00240 i__3 = *ldab - 1; 00241 i__4 = *ldab - 1; 00242 dger_(&km, &i__2, &c_b9, &ab[kv + 2 + j * ab_dim1], &c__1, 00243 &ab[kv + (j + 1) * ab_dim1], &i__3, &ab[kv + 1 + 00244 (j + 1) * ab_dim1], &i__4); 00245 } 00246 } 00247 } else { 00248 00249 /* If pivot is zero, set INFO to the index of the pivot */ 00250 /* unless a zero pivot has already been found. */ 00251 00252 if (*info == 0) { 00253 *info = j; 00254 } 00255 } 00256 /* L40: */ 00257 } 00258 return 0; 00259 00260 /* End of DGBTF2 */ 00261 00262 } /* dgbtf2_ */