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