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