![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b CGBEQUB 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download CGBEQUB + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgbequb.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgbequb.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgbequb.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE CGBEQUB( M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND, 00022 * AMAX, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER INFO, KL, KU, LDAB, M, N 00026 * REAL AMAX, COLCND, ROWCND 00027 * .. 00028 * .. Array Arguments .. 00029 * REAL C( * ), R( * ) 00030 * COMPLEX AB( LDAB, * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> CGBEQUB computes row and column scalings intended to equilibrate an 00040 *> M-by-N matrix A and reduce its condition number. R returns the row 00041 *> scale factors and C the column scale factors, chosen to try to make 00042 *> the largest element in each row and column of the matrix B with 00043 *> elements B(i,j)=R(i)*A(i,j)*C(j) have an absolute value of at most 00044 *> the radix. 00045 *> 00046 *> R(i) and C(j) are restricted to be a power of the radix between 00047 *> SMLNUM = smallest safe number and BIGNUM = largest safe number. Use 00048 *> of these scaling factors is not guaranteed to reduce the condition 00049 *> number of A but works well in practice. 00050 *> 00051 *> This routine differs from CGEEQU by restricting the scaling factors 00052 *> to a power of the radix. Baring over- and underflow, scaling by 00053 *> these factors introduces no additional rounding errors. However, the 00054 *> scaled entries' magnitured are no longer approximately 1 but lie 00055 *> between sqrt(radix) and 1/sqrt(radix). 00056 *> \endverbatim 00057 * 00058 * Arguments: 00059 * ========== 00060 * 00061 *> \param[in] M 00062 *> \verbatim 00063 *> M is INTEGER 00064 *> The number of rows of the matrix A. M >= 0. 00065 *> \endverbatim 00066 *> 00067 *> \param[in] N 00068 *> \verbatim 00069 *> N is INTEGER 00070 *> The number of columns of the matrix A. N >= 0. 00071 *> \endverbatim 00072 *> 00073 *> \param[in] KL 00074 *> \verbatim 00075 *> KL is INTEGER 00076 *> The number of subdiagonals within the band of A. KL >= 0. 00077 *> \endverbatim 00078 *> 00079 *> \param[in] KU 00080 *> \verbatim 00081 *> KU is INTEGER 00082 *> The number of superdiagonals within the band of A. KU >= 0. 00083 *> \endverbatim 00084 *> 00085 *> \param[in] AB 00086 *> \verbatim 00087 *> AB is 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 *> \endverbatim 00093 *> 00094 *> \param[in] LDAB 00095 *> \verbatim 00096 *> LDAB is INTEGER 00097 *> The leading dimension of the array A. LDAB >= max(1,M). 00098 *> \endverbatim 00099 *> 00100 *> \param[out] R 00101 *> \verbatim 00102 *> R is REAL array, dimension (M) 00103 *> If INFO = 0 or INFO > M, R contains the row scale factors 00104 *> for A. 00105 *> \endverbatim 00106 *> 00107 *> \param[out] C 00108 *> \verbatim 00109 *> C is REAL array, dimension (N) 00110 *> If INFO = 0, C contains the column scale factors for A. 00111 *> \endverbatim 00112 *> 00113 *> \param[out] ROWCND 00114 *> \verbatim 00115 *> ROWCND is REAL 00116 *> If INFO = 0 or INFO > M, ROWCND contains the ratio of the 00117 *> smallest R(i) to the largest R(i). If ROWCND >= 0.1 and 00118 *> AMAX is neither too large nor too small, it is not worth 00119 *> scaling by R. 00120 *> \endverbatim 00121 *> 00122 *> \param[out] COLCND 00123 *> \verbatim 00124 *> COLCND is REAL 00125 *> If INFO = 0, COLCND contains the ratio of the smallest 00126 *> C(i) to the largest C(i). If COLCND >= 0.1, it is not 00127 *> worth scaling by C. 00128 *> \endverbatim 00129 *> 00130 *> \param[out] AMAX 00131 *> \verbatim 00132 *> AMAX is REAL 00133 *> Absolute value of largest matrix element. If AMAX is very 00134 *> close to overflow or very close to underflow, the matrix 00135 *> should be scaled. 00136 *> \endverbatim 00137 *> 00138 *> \param[out] INFO 00139 *> \verbatim 00140 *> INFO is INTEGER 00141 *> = 0: successful exit 00142 *> < 0: if INFO = -i, the i-th argument had an illegal value 00143 *> > 0: if INFO = i, and i is 00144 *> <= M: the i-th row of A is exactly zero 00145 *> > M: the (i-M)-th column of A is exactly zero 00146 *> \endverbatim 00147 * 00148 * Authors: 00149 * ======== 00150 * 00151 *> \author Univ. of Tennessee 00152 *> \author Univ. of California Berkeley 00153 *> \author Univ. of Colorado Denver 00154 *> \author NAG Ltd. 00155 * 00156 *> \date November 2011 00157 * 00158 *> \ingroup complexGBcomputational 00159 * 00160 * ===================================================================== 00161 SUBROUTINE CGBEQUB( M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND, 00162 $ AMAX, INFO ) 00163 * 00164 * -- LAPACK computational routine (version 3.4.0) -- 00165 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00166 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00167 * November 2011 00168 * 00169 * .. Scalar Arguments .. 00170 INTEGER INFO, KL, KU, LDAB, M, N 00171 REAL AMAX, COLCND, ROWCND 00172 * .. 00173 * .. Array Arguments .. 00174 REAL C( * ), R( * ) 00175 COMPLEX AB( LDAB, * ) 00176 * .. 00177 * 00178 * ===================================================================== 00179 * 00180 * .. Parameters .. 00181 REAL ONE, ZERO 00182 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00183 * .. 00184 * .. Local Scalars .. 00185 INTEGER I, J, KD 00186 REAL BIGNUM, RCMAX, RCMIN, SMLNUM, RADIX, 00187 $ LOGRDX 00188 COMPLEX ZDUM 00189 * .. 00190 * .. External Functions .. 00191 REAL SLAMCH 00192 EXTERNAL SLAMCH 00193 * .. 00194 * .. External Subroutines .. 00195 EXTERNAL XERBLA 00196 * .. 00197 * .. Intrinsic Functions .. 00198 INTRINSIC ABS, MAX, MIN, LOG, REAL, AIMAG 00199 * .. 00200 * .. Statement Functions .. 00201 REAL CABS1 00202 * .. 00203 * .. Statement Function definitions .. 00204 CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) ) 00205 * .. 00206 * .. Executable Statements .. 00207 * 00208 * Test the input parameters. 00209 * 00210 INFO = 0 00211 IF( M.LT.0 ) THEN 00212 INFO = -1 00213 ELSE IF( N.LT.0 ) THEN 00214 INFO = -2 00215 ELSE IF( KL.LT.0 ) THEN 00216 INFO = -3 00217 ELSE IF( KU.LT.0 ) THEN 00218 INFO = -4 00219 ELSE IF( LDAB.LT.KL+KU+1 ) THEN 00220 INFO = -6 00221 END IF 00222 IF( INFO.NE.0 ) THEN 00223 CALL XERBLA( 'CGBEQUB', -INFO ) 00224 RETURN 00225 END IF 00226 * 00227 * Quick return if possible. 00228 * 00229 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00230 ROWCND = ONE 00231 COLCND = ONE 00232 AMAX = ZERO 00233 RETURN 00234 END IF 00235 * 00236 * Get machine constants. Assume SMLNUM is a power of the radix. 00237 * 00238 SMLNUM = SLAMCH( 'S' ) 00239 BIGNUM = ONE / SMLNUM 00240 RADIX = SLAMCH( 'B' ) 00241 LOGRDX = LOG(RADIX) 00242 * 00243 * Compute row scale factors. 00244 * 00245 DO 10 I = 1, M 00246 R( I ) = ZERO 00247 10 CONTINUE 00248 * 00249 * Find the maximum element in each row. 00250 * 00251 KD = KU + 1 00252 DO 30 J = 1, N 00253 DO 20 I = MAX( J-KU, 1 ), MIN( J+KL, M ) 00254 R( I ) = MAX( R( I ), CABS1( AB( KD+I-J, J ) ) ) 00255 20 CONTINUE 00256 30 CONTINUE 00257 DO I = 1, M 00258 IF( R( I ).GT.ZERO ) THEN 00259 R( I ) = RADIX**INT( LOG( R( I ) ) / LOGRDX ) 00260 END IF 00261 END DO 00262 * 00263 * Find the maximum and minimum scale factors. 00264 * 00265 RCMIN = BIGNUM 00266 RCMAX = ZERO 00267 DO 40 I = 1, M 00268 RCMAX = MAX( RCMAX, R( I ) ) 00269 RCMIN = MIN( RCMIN, R( I ) ) 00270 40 CONTINUE 00271 AMAX = RCMAX 00272 * 00273 IF( RCMIN.EQ.ZERO ) THEN 00274 * 00275 * Find the first zero scale factor and return an error code. 00276 * 00277 DO 50 I = 1, M 00278 IF( R( I ).EQ.ZERO ) THEN 00279 INFO = I 00280 RETURN 00281 END IF 00282 50 CONTINUE 00283 ELSE 00284 * 00285 * Invert the scale factors. 00286 * 00287 DO 60 I = 1, M 00288 R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM ) 00289 60 CONTINUE 00290 * 00291 * Compute ROWCND = min(R(I)) / max(R(I)). 00292 * 00293 ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM ) 00294 END IF 00295 * 00296 * Compute column scale factors. 00297 * 00298 DO 70 J = 1, N 00299 C( J ) = ZERO 00300 70 CONTINUE 00301 * 00302 * Find the maximum element in each column, 00303 * assuming the row scaling computed above. 00304 * 00305 DO 90 J = 1, N 00306 DO 80 I = MAX( J-KU, 1 ), MIN( J+KL, M ) 00307 C( J ) = MAX( C( J ), CABS1( AB( KD+I-J, J ) )*R( I ) ) 00308 80 CONTINUE 00309 IF( C( J ).GT.ZERO ) THEN 00310 C( J ) = RADIX**INT( LOG( C( J ) ) / LOGRDX ) 00311 END IF 00312 90 CONTINUE 00313 * 00314 * Find the maximum and minimum scale factors. 00315 * 00316 RCMIN = BIGNUM 00317 RCMAX = ZERO 00318 DO 100 J = 1, N 00319 RCMIN = MIN( RCMIN, C( J ) ) 00320 RCMAX = MAX( RCMAX, C( J ) ) 00321 100 CONTINUE 00322 * 00323 IF( RCMIN.EQ.ZERO ) THEN 00324 * 00325 * Find the first zero scale factor and return an error code. 00326 * 00327 DO 110 J = 1, N 00328 IF( C( J ).EQ.ZERO ) THEN 00329 INFO = M + J 00330 RETURN 00331 END IF 00332 110 CONTINUE 00333 ELSE 00334 * 00335 * Invert the scale factors. 00336 * 00337 DO 120 J = 1, N 00338 C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM ) 00339 120 CONTINUE 00340 * 00341 * Compute COLCND = min(C(J)) / max(C(J)). 00342 * 00343 COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM ) 00344 END IF 00345 * 00346 RETURN 00347 * 00348 * End of CGBEQUB 00349 * 00350 END