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