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