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