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