00001 /* zgbequ.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 /* Subroutine */ int zgbequ_(integer *m, integer *n, integer *kl, integer *ku, 00017 doublecomplex *ab, integer *ldab, doublereal *r__, doublereal *c__, 00018 doublereal *rowcnd, doublereal *colcnd, doublereal *amax, integer * 00019 info) 00020 { 00021 /* System generated locals */ 00022 integer ab_dim1, ab_offset, i__1, i__2, i__3, i__4; 00023 doublereal d__1, d__2, d__3, d__4; 00024 00025 /* Builtin functions */ 00026 double d_imag(doublecomplex *); 00027 00028 /* Local variables */ 00029 integer i__, j, kd; 00030 doublereal rcmin, rcmax; 00031 extern doublereal dlamch_(char *); 00032 extern /* Subroutine */ int xerbla_(char *, integer *); 00033 doublereal bignum, smlnum; 00034 00035 00036 /* -- LAPACK routine (version 3.2) -- */ 00037 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 00038 /* November 2006 */ 00039 00040 /* .. Scalar Arguments .. */ 00041 /* .. */ 00042 /* .. Array Arguments .. */ 00043 /* .. */ 00044 00045 /* Purpose */ 00046 /* ======= */ 00047 00048 /* ZGBEQU computes row and column scalings intended to equilibrate an */ 00049 /* M-by-N band matrix A and reduce its condition number. R returns the */ 00050 /* row scale factors and C the column scale factors, chosen to try to */ 00051 /* make the largest element in each row and column of the matrix B with */ 00052 /* elements B(i,j)=R(i)*A(i,j)*C(j) have absolute value 1. */ 00053 00054 /* R(i) and C(j) are restricted to be between SMLNUM = smallest safe */ 00055 /* number and BIGNUM = largest safe number. Use of these scaling */ 00056 /* factors is not guaranteed to reduce the condition number of A but */ 00057 /* works well in practice. */ 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) COMPLEX*16 array, dimension (LDAB,N) */ 00075 /* The band matrix A, stored in rows 1 to KL+KU+1. The j-th */ 00076 /* column of A is stored in the j-th column of the array AB as */ 00077 /* follows: */ 00078 /* AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl). */ 00079 00080 /* LDAB (input) INTEGER */ 00081 /* The leading dimension of the array AB. LDAB >= KL+KU+1. */ 00082 00083 /* R (output) DOUBLE PRECISION array, dimension (M) */ 00084 /* If INFO = 0, or INFO > M, R contains the row scale factors */ 00085 /* for A. */ 00086 00087 /* C (output) DOUBLE PRECISION array, dimension (N) */ 00088 /* If INFO = 0, C contains the column scale factors for A. */ 00089 00090 /* ROWCND (output) DOUBLE PRECISION */ 00091 /* If INFO = 0 or INFO > M, ROWCND contains the ratio of the */ 00092 /* smallest R(i) to the largest R(i). If ROWCND >= 0.1 and */ 00093 /* AMAX is neither too large nor too small, it is not worth */ 00094 /* scaling by R. */ 00095 00096 /* COLCND (output) DOUBLE PRECISION */ 00097 /* If INFO = 0, COLCND contains the ratio of the smallest */ 00098 /* C(i) to the largest C(i). If COLCND >= 0.1, it is not */ 00099 /* worth scaling by C. */ 00100 00101 /* AMAX (output) DOUBLE PRECISION */ 00102 /* Absolute value of largest matrix element. If AMAX is very */ 00103 /* close to overflow or very close to underflow, the matrix */ 00104 /* should be scaled. */ 00105 00106 /* INFO (output) INTEGER */ 00107 /* = 0: successful exit */ 00108 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 00109 /* > 0: if INFO = i, and i is */ 00110 /* <= M: the i-th row of A is exactly zero */ 00111 /* > M: the (i-M)-th column of A is exactly zero */ 00112 00113 /* ===================================================================== */ 00114 00115 /* .. Parameters .. */ 00116 /* .. */ 00117 /* .. Local Scalars .. */ 00118 /* .. */ 00119 /* .. External Functions .. */ 00120 /* .. */ 00121 /* .. External Subroutines .. */ 00122 /* .. */ 00123 /* .. Intrinsic Functions .. */ 00124 /* .. */ 00125 /* .. Statement Functions .. */ 00126 /* .. */ 00127 /* .. Statement Function definitions .. */ 00128 /* .. */ 00129 /* .. Executable Statements .. */ 00130 00131 /* Test the input parameters */ 00132 00133 /* Parameter adjustments */ 00134 ab_dim1 = *ldab; 00135 ab_offset = 1 + ab_dim1; 00136 ab -= ab_offset; 00137 --r__; 00138 --c__; 00139 00140 /* Function Body */ 00141 *info = 0; 00142 if (*m < 0) { 00143 *info = -1; 00144 } else if (*n < 0) { 00145 *info = -2; 00146 } else if (*kl < 0) { 00147 *info = -3; 00148 } else if (*ku < 0) { 00149 *info = -4; 00150 } else if (*ldab < *kl + *ku + 1) { 00151 *info = -6; 00152 } 00153 if (*info != 0) { 00154 i__1 = -(*info); 00155 xerbla_("ZGBEQU", &i__1); 00156 return 0; 00157 } 00158 00159 /* Quick return if possible */ 00160 00161 if (*m == 0 || *n == 0) { 00162 *rowcnd = 1.; 00163 *colcnd = 1.; 00164 *amax = 0.; 00165 return 0; 00166 } 00167 00168 /* Get machine constants. */ 00169 00170 smlnum = dlamch_("S"); 00171 bignum = 1. / smlnum; 00172 00173 /* Compute row scale factors. */ 00174 00175 i__1 = *m; 00176 for (i__ = 1; i__ <= i__1; ++i__) { 00177 r__[i__] = 0.; 00178 /* L10: */ 00179 } 00180 00181 /* Find the maximum element in each row. */ 00182 00183 kd = *ku + 1; 00184 i__1 = *n; 00185 for (j = 1; j <= i__1; ++j) { 00186 /* Computing MAX */ 00187 i__2 = j - *ku; 00188 /* Computing MIN */ 00189 i__4 = j + *kl; 00190 i__3 = min(i__4,*m); 00191 for (i__ = max(i__2,1); i__ <= i__3; ++i__) { 00192 /* Computing MAX */ 00193 i__2 = kd + i__ - j + j * ab_dim1; 00194 d__3 = r__[i__], d__4 = (d__1 = ab[i__2].r, abs(d__1)) + (d__2 = 00195 d_imag(&ab[kd + i__ - j + j * ab_dim1]), abs(d__2)); 00196 r__[i__] = max(d__3,d__4); 00197 /* L20: */ 00198 } 00199 /* L30: */ 00200 } 00201 00202 /* Find the maximum and minimum scale factors. */ 00203 00204 rcmin = bignum; 00205 rcmax = 0.; 00206 i__1 = *m; 00207 for (i__ = 1; i__ <= i__1; ++i__) { 00208 /* Computing MAX */ 00209 d__1 = rcmax, d__2 = r__[i__]; 00210 rcmax = max(d__1,d__2); 00211 /* Computing MIN */ 00212 d__1 = rcmin, d__2 = r__[i__]; 00213 rcmin = min(d__1,d__2); 00214 /* L40: */ 00215 } 00216 *amax = rcmax; 00217 00218 if (rcmin == 0.) { 00219 00220 /* Find the first zero scale factor and return an error code. */ 00221 00222 i__1 = *m; 00223 for (i__ = 1; i__ <= i__1; ++i__) { 00224 if (r__[i__] == 0.) { 00225 *info = i__; 00226 return 0; 00227 } 00228 /* L50: */ 00229 } 00230 } else { 00231 00232 /* Invert the scale factors. */ 00233 00234 i__1 = *m; 00235 for (i__ = 1; i__ <= i__1; ++i__) { 00236 /* Computing MIN */ 00237 /* Computing MAX */ 00238 d__2 = r__[i__]; 00239 d__1 = max(d__2,smlnum); 00240 r__[i__] = 1. / min(d__1,bignum); 00241 /* L60: */ 00242 } 00243 00244 /* Compute ROWCND = min(R(I)) / max(R(I)) */ 00245 00246 *rowcnd = max(rcmin,smlnum) / min(rcmax,bignum); 00247 } 00248 00249 /* Compute column scale factors */ 00250 00251 i__1 = *n; 00252 for (j = 1; j <= i__1; ++j) { 00253 c__[j] = 0.; 00254 /* L70: */ 00255 } 00256 00257 /* Find the maximum element in each column, */ 00258 /* assuming the row scaling computed above. */ 00259 00260 kd = *ku + 1; 00261 i__1 = *n; 00262 for (j = 1; j <= i__1; ++j) { 00263 /* Computing MAX */ 00264 i__3 = j - *ku; 00265 /* Computing MIN */ 00266 i__4 = j + *kl; 00267 i__2 = min(i__4,*m); 00268 for (i__ = max(i__3,1); i__ <= i__2; ++i__) { 00269 /* Computing MAX */ 00270 i__3 = kd + i__ - j + j * ab_dim1; 00271 d__3 = c__[j], d__4 = ((d__1 = ab[i__3].r, abs(d__1)) + (d__2 = 00272 d_imag(&ab[kd + i__ - j + j * ab_dim1]), abs(d__2))) * 00273 r__[i__]; 00274 c__[j] = max(d__3,d__4); 00275 /* L80: */ 00276 } 00277 /* L90: */ 00278 } 00279 00280 /* Find the maximum and minimum scale factors. */ 00281 00282 rcmin = bignum; 00283 rcmax = 0.; 00284 i__1 = *n; 00285 for (j = 1; j <= i__1; ++j) { 00286 /* Computing MIN */ 00287 d__1 = rcmin, d__2 = c__[j]; 00288 rcmin = min(d__1,d__2); 00289 /* Computing MAX */ 00290 d__1 = rcmax, d__2 = c__[j]; 00291 rcmax = max(d__1,d__2); 00292 /* L100: */ 00293 } 00294 00295 if (rcmin == 0.) { 00296 00297 /* Find the first zero scale factor and return an error code. */ 00298 00299 i__1 = *n; 00300 for (j = 1; j <= i__1; ++j) { 00301 if (c__[j] == 0.) { 00302 *info = *m + j; 00303 return 0; 00304 } 00305 /* L110: */ 00306 } 00307 } else { 00308 00309 /* Invert the scale factors. */ 00310 00311 i__1 = *n; 00312 for (j = 1; j <= i__1; ++j) { 00313 /* Computing MIN */ 00314 /* Computing MAX */ 00315 d__2 = c__[j]; 00316 d__1 = max(d__2,smlnum); 00317 c__[j] = 1. / min(d__1,bignum); 00318 /* L120: */ 00319 } 00320 00321 /* Compute COLCND = min(C(J)) / max(C(J)) */ 00322 00323 *colcnd = max(rcmin,smlnum) / min(rcmax,bignum); 00324 } 00325 00326 return 0; 00327 00328 /* End of ZGBEQU */ 00329 00330 } /* zgbequ_ */