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