00001 SUBROUTINE SGEEQUB( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX, 00002 $ INFO ) 00003 * 00004 * -- LAPACK routine (version 3.2) -- 00005 * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- 00006 * -- Jason Riedy of Univ. of California Berkeley. -- 00007 * -- November 2008 -- 00008 * 00009 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00010 * -- Univ. of California Berkeley and NAG Ltd. -- 00011 * 00012 IMPLICIT NONE 00013 * .. 00014 * .. Scalar Arguments .. 00015 INTEGER INFO, LDA, M, N 00016 REAL AMAX, COLCND, ROWCND 00017 * .. 00018 * .. Array Arguments .. 00019 REAL A( LDA, * ), C( * ), R( * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * SGEEQUB computes row and column scalings intended to equilibrate an 00026 * M-by-N matrix A and reduce its condition number. R returns the row 00027 * scale factors and C the column scale factors, chosen to try to make 00028 * the largest element in each row and column of the matrix B with 00029 * elements B(i,j)=R(i)*A(i,j)*C(j) have an absolute value of at most 00030 * the radix. 00031 * 00032 * R(i) and C(j) are restricted to be a power of the radix between 00033 * SMLNUM = smallest safe number and BIGNUM = largest safe number. Use 00034 * of these scaling factors is not guaranteed to reduce the condition 00035 * number of A but works well in practice. 00036 * 00037 * This routine differs from SGEEQU by restricting the scaling factors 00038 * to a power of the radix. Baring over- and underflow, scaling by 00039 * these factors introduces no additional rounding errors. However, the 00040 * scaled entries' magnitured are no longer approximately 1 but lie 00041 * between sqrt(radix) and 1/sqrt(radix). 00042 * 00043 * Arguments 00044 * ========= 00045 * 00046 * M (input) INTEGER 00047 * The number of rows of the matrix A. M >= 0. 00048 * 00049 * N (input) INTEGER 00050 * The number of columns of the matrix A. N >= 0. 00051 * 00052 * A (input) REAL array, dimension (LDA,N) 00053 * The M-by-N matrix whose equilibration factors are 00054 * to be computed. 00055 * 00056 * LDA (input) INTEGER 00057 * The leading dimension of the array A. LDA >= max(1,M). 00058 * 00059 * R (output) REAL array, dimension (M) 00060 * If INFO = 0 or INFO > M, R contains the row scale factors 00061 * for A. 00062 * 00063 * C (output) REAL array, dimension (N) 00064 * If INFO = 0, C contains the column scale factors for A. 00065 * 00066 * ROWCND (output) REAL 00067 * If INFO = 0 or INFO > M, ROWCND contains the ratio of the 00068 * smallest R(i) to the largest R(i). If ROWCND >= 0.1 and 00069 * AMAX is neither too large nor too small, it is not worth 00070 * scaling by R. 00071 * 00072 * COLCND (output) REAL 00073 * If INFO = 0, COLCND contains the ratio of the smallest 00074 * C(i) to the largest C(i). If COLCND >= 0.1, it is not 00075 * worth scaling by C. 00076 * 00077 * AMAX (output) REAL 00078 * Absolute value of largest matrix element. If AMAX is very 00079 * close to overflow or very close to underflow, the matrix 00080 * should be scaled. 00081 * 00082 * INFO (output) INTEGER 00083 * = 0: successful exit 00084 * < 0: if INFO = -i, the i-th argument had an illegal value 00085 * > 0: if INFO = i, and i is 00086 * <= M: the i-th row of A is exactly zero 00087 * > M: the (i-M)-th column of A is exactly zero 00088 * 00089 * ===================================================================== 00090 * 00091 * .. Parameters .. 00092 REAL ONE, ZERO 00093 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00094 * .. 00095 * .. Local Scalars .. 00096 INTEGER I, J 00097 REAL BIGNUM, RCMAX, RCMIN, SMLNUM, RADIX, LOGRDX 00098 * .. 00099 * .. External Functions .. 00100 REAL SLAMCH 00101 EXTERNAL SLAMCH 00102 * .. 00103 * .. External Subroutines .. 00104 EXTERNAL XERBLA 00105 * .. 00106 * .. Intrinsic Functions .. 00107 INTRINSIC ABS, MAX, MIN, LOG 00108 * .. 00109 * .. Executable Statements .. 00110 * 00111 * Test the input parameters. 00112 * 00113 INFO = 0 00114 IF( M.LT.0 ) THEN 00115 INFO = -1 00116 ELSE IF( N.LT.0 ) THEN 00117 INFO = -2 00118 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 00119 INFO = -4 00120 END IF 00121 IF( INFO.NE.0 ) THEN 00122 CALL XERBLA( 'SGEEQUB', -INFO ) 00123 RETURN 00124 END IF 00125 * 00126 * Quick return if possible. 00127 * 00128 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00129 ROWCND = ONE 00130 COLCND = ONE 00131 AMAX = ZERO 00132 RETURN 00133 END IF 00134 * 00135 * Get machine constants. Assume SMLNUM is a power of the radix. 00136 * 00137 SMLNUM = SLAMCH( 'S' ) 00138 BIGNUM = ONE / SMLNUM 00139 RADIX = SLAMCH( 'B' ) 00140 LOGRDX = LOG( RADIX ) 00141 * 00142 * Compute row scale factors. 00143 * 00144 DO 10 I = 1, M 00145 R( I ) = ZERO 00146 10 CONTINUE 00147 * 00148 * Find the maximum element in each row. 00149 * 00150 DO 30 J = 1, N 00151 DO 20 I = 1, M 00152 R( I ) = MAX( R( I ), ABS( A( I, J ) ) ) 00153 20 CONTINUE 00154 30 CONTINUE 00155 DO I = 1, M 00156 IF( R( I ).GT.ZERO ) THEN 00157 R( I ) = RADIX**INT( LOG( R( I ) ) / LOGRDX ) 00158 END IF 00159 END DO 00160 * 00161 * Find the maximum and minimum scale factors. 00162 * 00163 RCMIN = BIGNUM 00164 RCMAX = ZERO 00165 DO 40 I = 1, M 00166 RCMAX = MAX( RCMAX, R( I ) ) 00167 RCMIN = MIN( RCMIN, R( I ) ) 00168 40 CONTINUE 00169 AMAX = RCMAX 00170 * 00171 IF( RCMIN.EQ.ZERO ) THEN 00172 * 00173 * Find the first zero scale factor and return an error code. 00174 * 00175 DO 50 I = 1, M 00176 IF( R( I ).EQ.ZERO ) THEN 00177 INFO = I 00178 RETURN 00179 END IF 00180 50 CONTINUE 00181 ELSE 00182 * 00183 * Invert the scale factors. 00184 * 00185 DO 60 I = 1, M 00186 R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM ) 00187 60 CONTINUE 00188 * 00189 * Compute ROWCND = min(R(I)) / max(R(I)). 00190 * 00191 ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM ) 00192 END IF 00193 * 00194 * Compute column scale factors 00195 * 00196 DO 70 J = 1, N 00197 C( J ) = ZERO 00198 70 CONTINUE 00199 * 00200 * Find the maximum element in each column, 00201 * assuming the row scaling computed above. 00202 * 00203 DO 90 J = 1, N 00204 DO 80 I = 1, M 00205 C( J ) = MAX( C( J ), ABS( A( I, J ) )*R( I ) ) 00206 80 CONTINUE 00207 IF( C( J ).GT.ZERO ) THEN 00208 C( J ) = RADIX**INT( LOG( C( J ) ) / LOGRDX ) 00209 END IF 00210 90 CONTINUE 00211 * 00212 * Find the maximum and minimum scale factors. 00213 * 00214 RCMIN = BIGNUM 00215 RCMAX = ZERO 00216 DO 100 J = 1, N 00217 RCMIN = MIN( RCMIN, C( J ) ) 00218 RCMAX = MAX( RCMAX, C( J ) ) 00219 100 CONTINUE 00220 * 00221 IF( RCMIN.EQ.ZERO ) THEN 00222 * 00223 * Find the first zero scale factor and return an error code. 00224 * 00225 DO 110 J = 1, N 00226 IF( C( J ).EQ.ZERO ) THEN 00227 INFO = M + J 00228 RETURN 00229 END IF 00230 110 CONTINUE 00231 ELSE 00232 * 00233 * Invert the scale factors. 00234 * 00235 DO 120 J = 1, N 00236 C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM ) 00237 120 CONTINUE 00238 * 00239 * Compute COLCND = min(C(J)) / max(C(J)). 00240 * 00241 COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM ) 00242 END IF 00243 * 00244 RETURN 00245 * 00246 * End of SGEEQUB 00247 * 00248 END
1.7.2