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