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