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