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