![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DGEEQUB 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DGEEQUB + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeequb.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeequb.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeequb.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DGEEQUB( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, 00022 * INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER INFO, LDA, M, N 00026 * DOUBLE PRECISION AMAX, COLCND, ROWCND 00027 * .. 00028 * .. Array Arguments .. 00029 * DOUBLE PRECISION A( LDA, * ), C( * ), R( * ) 00030 * .. 00031 * 00032 * 00033 *> \par Purpose: 00034 * ============= 00035 *> 00036 *> \verbatim 00037 *> 00038 *> DGEEQUB 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] A 00073 *> \verbatim 00074 *> A is DOUBLE PRECISION array, dimension (LDA,N) 00075 *> The M-by-N matrix whose equilibration factors are 00076 *> to be computed. 00077 *> \endverbatim 00078 *> 00079 *> \param[in] LDA 00080 *> \verbatim 00081 *> LDA is INTEGER 00082 *> The leading dimension of the array A. LDA >= max(1,M). 00083 *> \endverbatim 00084 *> 00085 *> \param[out] R 00086 *> \verbatim 00087 *> R is DOUBLE PRECISION array, dimension (M) 00088 *> If INFO = 0 or INFO > M, R contains the row scale factors 00089 *> for A. 00090 *> \endverbatim 00091 *> 00092 *> \param[out] C 00093 *> \verbatim 00094 *> C is DOUBLE PRECISION array, dimension (N) 00095 *> If INFO = 0, C contains the column scale factors for A. 00096 *> \endverbatim 00097 *> 00098 *> \param[out] ROWCND 00099 *> \verbatim 00100 *> ROWCND is DOUBLE PRECISION 00101 *> If INFO = 0 or INFO > M, ROWCND contains the ratio of the 00102 *> smallest R(i) to the largest R(i). If ROWCND >= 0.1 and 00103 *> AMAX is neither too large nor too small, it is not worth 00104 *> scaling by R. 00105 *> \endverbatim 00106 *> 00107 *> \param[out] COLCND 00108 *> \verbatim 00109 *> COLCND is 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 *> \endverbatim 00114 *> 00115 *> \param[out] AMAX 00116 *> \verbatim 00117 *> AMAX is DOUBLE PRECISION 00118 *> Absolute value of largest matrix element. If AMAX is very 00119 *> close to overflow or very close to underflow, the matrix 00120 *> should be scaled. 00121 *> \endverbatim 00122 *> 00123 *> \param[out] INFO 00124 *> \verbatim 00125 *> INFO is INTEGER 00126 *> = 0: successful exit 00127 *> < 0: if INFO = -i, the i-th argument had an illegal value 00128 *> > 0: if INFO = i, and i is 00129 *> <= M: the i-th row of A is exactly zero 00130 *> > M: the (i-M)-th column of A is exactly zero 00131 *> \endverbatim 00132 * 00133 * Authors: 00134 * ======== 00135 * 00136 *> \author Univ. of Tennessee 00137 *> \author Univ. of California Berkeley 00138 *> \author Univ. of Colorado Denver 00139 *> \author NAG Ltd. 00140 * 00141 *> \date November 2011 00142 * 00143 *> \ingroup doubleGEcomputational 00144 * 00145 * ===================================================================== 00146 SUBROUTINE DGEEQUB( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, 00147 $ INFO ) 00148 * 00149 * -- LAPACK computational routine (version 3.4.0) -- 00150 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00151 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00152 * November 2011 00153 * 00154 * .. Scalar Arguments .. 00155 INTEGER INFO, LDA, M, N 00156 DOUBLE PRECISION AMAX, COLCND, ROWCND 00157 * .. 00158 * .. Array Arguments .. 00159 DOUBLE PRECISION A( LDA, * ), C( * ), R( * ) 00160 * .. 00161 * 00162 * ===================================================================== 00163 * 00164 * .. Parameters .. 00165 DOUBLE PRECISION ONE, ZERO 00166 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00167 * .. 00168 * .. Local Scalars .. 00169 INTEGER I, J 00170 DOUBLE PRECISION BIGNUM, RCMAX, RCMIN, SMLNUM, RADIX, LOGRDX 00171 * .. 00172 * .. External Functions .. 00173 DOUBLE PRECISION DLAMCH 00174 EXTERNAL DLAMCH 00175 * .. 00176 * .. External Subroutines .. 00177 EXTERNAL XERBLA 00178 * .. 00179 * .. Intrinsic Functions .. 00180 INTRINSIC ABS, MAX, MIN, LOG 00181 * .. 00182 * .. Executable Statements .. 00183 * 00184 * Test the input parameters. 00185 * 00186 INFO = 0 00187 IF( M.LT.0 ) THEN 00188 INFO = -1 00189 ELSE IF( N.LT.0 ) THEN 00190 INFO = -2 00191 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00192 INFO = -4 00193 END IF 00194 IF( INFO.NE.0 ) THEN 00195 CALL XERBLA( 'DGEEQUB', -INFO ) 00196 RETURN 00197 END IF 00198 * 00199 * Quick return if possible. 00200 * 00201 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00202 ROWCND = ONE 00203 COLCND = ONE 00204 AMAX = ZERO 00205 RETURN 00206 END IF 00207 * 00208 * Get machine constants. Assume SMLNUM is a power of the radix. 00209 * 00210 SMLNUM = DLAMCH( 'S' ) 00211 BIGNUM = ONE / SMLNUM 00212 RADIX = DLAMCH( 'B' ) 00213 LOGRDX = LOG( RADIX ) 00214 * 00215 * Compute row scale factors. 00216 * 00217 DO 10 I = 1, M 00218 R( I ) = ZERO 00219 10 CONTINUE 00220 * 00221 * Find the maximum element in each row. 00222 * 00223 DO 30 J = 1, N 00224 DO 20 I = 1, M 00225 R( I ) = MAX( R( I ), ABS( A( I, J ) ) ) 00226 20 CONTINUE 00227 30 CONTINUE 00228 DO I = 1, M 00229 IF( R( I ).GT.ZERO ) THEN 00230 R( I ) = RADIX**INT( LOG( R( I ) ) / LOGRDX ) 00231 END IF 00232 END DO 00233 * 00234 * Find the maximum and minimum scale factors. 00235 * 00236 RCMIN = BIGNUM 00237 RCMAX = ZERO 00238 DO 40 I = 1, M 00239 RCMAX = MAX( RCMAX, R( I ) ) 00240 RCMIN = MIN( RCMIN, R( I ) ) 00241 40 CONTINUE 00242 AMAX = RCMAX 00243 * 00244 IF( RCMIN.EQ.ZERO ) THEN 00245 * 00246 * Find the first zero scale factor and return an error code. 00247 * 00248 DO 50 I = 1, M 00249 IF( R( I ).EQ.ZERO ) THEN 00250 INFO = I 00251 RETURN 00252 END IF 00253 50 CONTINUE 00254 ELSE 00255 * 00256 * Invert the scale factors. 00257 * 00258 DO 60 I = 1, M 00259 R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM ) 00260 60 CONTINUE 00261 * 00262 * Compute ROWCND = min(R(I)) / max(R(I)). 00263 * 00264 ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM ) 00265 END IF 00266 * 00267 * Compute column scale factors 00268 * 00269 DO 70 J = 1, N 00270 C( J ) = ZERO 00271 70 CONTINUE 00272 * 00273 * Find the maximum element in each column, 00274 * assuming the row scaling computed above. 00275 * 00276 DO 90 J = 1, N 00277 DO 80 I = 1, M 00278 C( J ) = MAX( C( J ), ABS( A( I, J ) )*R( I ) ) 00279 80 CONTINUE 00280 IF( C( J ).GT.ZERO ) THEN 00281 C( J ) = RADIX**INT( LOG( C( J ) ) / LOGRDX ) 00282 END IF 00283 90 CONTINUE 00284 * 00285 * Find the maximum and minimum scale factors. 00286 * 00287 RCMIN = BIGNUM 00288 RCMAX = ZERO 00289 DO 100 J = 1, N 00290 RCMIN = MIN( RCMIN, C( J ) ) 00291 RCMAX = MAX( RCMAX, C( J ) ) 00292 100 CONTINUE 00293 * 00294 IF( RCMIN.EQ.ZERO ) THEN 00295 * 00296 * Find the first zero scale factor and return an error code. 00297 * 00298 DO 110 J = 1, N 00299 IF( C( J ).EQ.ZERO ) THEN 00300 INFO = M + J 00301 RETURN 00302 END IF 00303 110 CONTINUE 00304 ELSE 00305 * 00306 * Invert the scale factors. 00307 * 00308 DO 120 J = 1, N 00309 C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM ) 00310 120 CONTINUE 00311 * 00312 * Compute COLCND = min(C(J)) / max(C(J)). 00313 * 00314 COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM ) 00315 END IF 00316 * 00317 RETURN 00318 * 00319 * End of DGEEQUB 00320 * 00321 END