00001 SUBROUTINE CGBEQUB( M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND, 00002 $ AMAX, 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, KL, KU, LDAB, M, N 00016 REAL AMAX, COLCND, ROWCND 00017 * .. 00018 * .. Array Arguments .. 00019 REAL C( * ), R( * ) 00020 COMPLEX AB( LDAB, * ) 00021 * .. 00022 * 00023 * Purpose 00024 * ======= 00025 * 00026 * CGBEQUB 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 * KL (input) INTEGER 00054 * The number of subdiagonals within the band of A. KL >= 0. 00055 * 00056 * KU (input) INTEGER 00057 * The number of superdiagonals within the band of A. KU >= 0. 00058 * 00059 * AB (input) DOUBLE PRECISION array, dimension (LDAB,N) 00060 * On entry, the matrix A in band storage, in rows 1 to KL+KU+1. 00061 * The j-th column of A is stored in the j-th column of the 00062 * array AB as follows: 00063 * AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl) 00064 * 00065 * LDAB (input) INTEGER 00066 * The leading dimension of the array A. LDAB >= max(1,M). 00067 * 00068 * R (output) REAL array, dimension (M) 00069 * If INFO = 0 or INFO > M, R contains the row scale factors 00070 * for A. 00071 * 00072 * C (output) REAL array, dimension (N) 00073 * If INFO = 0, C contains the column scale factors for A. 00074 * 00075 * ROWCND (output) REAL 00076 * If INFO = 0 or INFO > M, ROWCND contains the ratio of the 00077 * smallest R(i) to the largest R(i). If ROWCND >= 0.1 and 00078 * AMAX is neither too large nor too small, it is not worth 00079 * scaling by R. 00080 * 00081 * COLCND (output) REAL 00082 * If INFO = 0, COLCND contains the ratio of the smallest 00083 * C(i) to the largest C(i). If COLCND >= 0.1, it is not 00084 * worth scaling by C. 00085 * 00086 * AMAX (output) REAL 00087 * Absolute value of largest matrix element. If AMAX is very 00088 * close to overflow or very close to underflow, the matrix 00089 * should be scaled. 00090 * 00091 * INFO (output) INTEGER 00092 * = 0: successful exit 00093 * < 0: if INFO = -i, the i-th argument had an illegal value 00094 * > 0: if INFO = i, and i is 00095 * <= M: the i-th row of A is exactly zero 00096 * > M: the (i-M)-th column of A is exactly zero 00097 * 00098 * ===================================================================== 00099 * 00100 * .. Parameters .. 00101 REAL ONE, ZERO 00102 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00103 * .. 00104 * .. Local Scalars .. 00105 INTEGER I, J, KD 00106 REAL BIGNUM, RCMAX, RCMIN, SMLNUM, RADIX, 00107 $ LOGRDX 00108 COMPLEX ZDUM 00109 * .. 00110 * .. External Functions .. 00111 REAL SLAMCH 00112 EXTERNAL SLAMCH 00113 * .. 00114 * .. External Subroutines .. 00115 EXTERNAL XERBLA 00116 * .. 00117 * .. Intrinsic Functions .. 00118 INTRINSIC ABS, MAX, MIN, LOG, REAL, AIMAG 00119 * .. 00120 * .. Statement Functions .. 00121 REAL CABS1 00122 * .. 00123 * .. Statement Function definitions .. 00124 CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) ) 00125 * .. 00126 * .. Executable Statements .. 00127 * 00128 * Test the input parameters. 00129 * 00130 INFO = 0 00131 IF( M.LT.0 ) THEN 00132 INFO = -1 00133 ELSE IF( N.LT.0 ) THEN 00134 INFO = -2 00135 ELSE IF( KL.LT.0 ) THEN 00136 INFO = -3 00137 ELSE IF( KU.LT.0 ) THEN 00138 INFO = -4 00139 ELSE IF( LDAB.LT.KL+KU+1 ) THEN 00140 INFO = -6 00141 END IF 00142 IF( INFO.NE.0 ) THEN 00143 CALL XERBLA( 'CGBEQUB', -INFO ) 00144 RETURN 00145 END IF 00146 * 00147 * Quick return if possible. 00148 * 00149 IF( M.EQ.0 .OR. N.EQ.0 ) THEN 00150 ROWCND = ONE 00151 COLCND = ONE 00152 AMAX = ZERO 00153 RETURN 00154 END IF 00155 * 00156 * Get machine constants. Assume SMLNUM is a power of the radix. 00157 * 00158 SMLNUM = SLAMCH( 'S' ) 00159 BIGNUM = ONE / SMLNUM 00160 RADIX = SLAMCH( 'B' ) 00161 LOGRDX = LOG(RADIX) 00162 * 00163 * Compute row scale factors. 00164 * 00165 DO 10 I = 1, M 00166 R( I ) = ZERO 00167 10 CONTINUE 00168 * 00169 * Find the maximum element in each row. 00170 * 00171 KD = KU + 1 00172 DO 30 J = 1, N 00173 DO 20 I = MAX( J-KU, 1 ), MIN( J+KL, M ) 00174 R( I ) = MAX( R( I ), CABS1( AB( KD+I-J, J ) ) ) 00175 20 CONTINUE 00176 30 CONTINUE 00177 DO I = 1, M 00178 IF( R( I ).GT.ZERO ) THEN 00179 R( I ) = RADIX**INT( LOG( R( I ) ) / LOGRDX ) 00180 END IF 00181 END DO 00182 * 00183 * Find the maximum and minimum scale factors. 00184 * 00185 RCMIN = BIGNUM 00186 RCMAX = ZERO 00187 DO 40 I = 1, M 00188 RCMAX = MAX( RCMAX, R( I ) ) 00189 RCMIN = MIN( RCMIN, R( I ) ) 00190 40 CONTINUE 00191 AMAX = RCMAX 00192 * 00193 IF( RCMIN.EQ.ZERO ) THEN 00194 * 00195 * Find the first zero scale factor and return an error code. 00196 * 00197 DO 50 I = 1, M 00198 IF( R( I ).EQ.ZERO ) THEN 00199 INFO = I 00200 RETURN 00201 END IF 00202 50 CONTINUE 00203 ELSE 00204 * 00205 * Invert the scale factors. 00206 * 00207 DO 60 I = 1, M 00208 R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM ) 00209 60 CONTINUE 00210 * 00211 * Compute ROWCND = min(R(I)) / max(R(I)). 00212 * 00213 ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM ) 00214 END IF 00215 * 00216 * Compute column scale factors. 00217 * 00218 DO 70 J = 1, N 00219 C( J ) = ZERO 00220 70 CONTINUE 00221 * 00222 * Find the maximum element in each column, 00223 * assuming the row scaling computed above. 00224 * 00225 DO 90 J = 1, N 00226 DO 80 I = MAX( J-KU, 1 ), MIN( J+KL, M ) 00227 C( J ) = MAX( C( J ), CABS1( AB( KD+I-J, J ) )*R( I ) ) 00228 80 CONTINUE 00229 IF( C( J ).GT.ZERO ) THEN 00230 C( J ) = RADIX**INT( LOG( C( J ) ) / LOGRDX ) 00231 END IF 00232 90 CONTINUE 00233 * 00234 * Find the maximum and minimum scale factors. 00235 * 00236 RCMIN = BIGNUM 00237 RCMAX = ZERO 00238 DO 100 J = 1, N 00239 RCMIN = MIN( RCMIN, C( J ) ) 00240 RCMAX = MAX( RCMAX, C( J ) ) 00241 100 CONTINUE 00242 * 00243 IF( RCMIN.EQ.ZERO ) THEN 00244 * 00245 * Find the first zero scale factor and return an error code. 00246 * 00247 DO 110 J = 1, N 00248 IF( C( J ).EQ.ZERO ) THEN 00249 INFO = M + J 00250 RETURN 00251 END IF 00252 110 CONTINUE 00253 ELSE 00254 * 00255 * Invert the scale factors. 00256 * 00257 DO 120 J = 1, N 00258 C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM ) 00259 120 CONTINUE 00260 * 00261 * Compute COLCND = min(C(J)) / max(C(J)). 00262 * 00263 COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM ) 00264 END IF 00265 * 00266 RETURN 00267 * 00268 * End of CGBEQUB 00269 * 00270 END
1.7.2