1 SUBROUTINE DGBEQUB( M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND,
2 $ AMAX, INFO )
3 *
4 * -- LAPACK routine (version 3.2) --
5 * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
6 * -- Jason Riedy of Univ. of California Berkeley. --
7 * -- November 2008 --
8 *
9 * -- LAPACK is a software package provided by Univ. of Tennessee, --
10 * -- Univ. of California Berkeley and NAG Ltd. --
11 *
12 IMPLICIT NONE
13 * ..
14 * .. Scalar Arguments ..
15 INTEGER INFO, KL, KU, LDAB, M, N
16 DOUBLE PRECISION AMAX, COLCND, ROWCND
17 * ..
18 * .. Array Arguments ..
19 DOUBLE PRECISION AB( LDAB, * ), C( * ), R( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * DGBEQUB computes row and column scalings intended to equilibrate an
26 * M-by-N matrix A and reduce its condition number. R returns the row
27 * scale factors and C the column scale factors, chosen to try to make
28 * the largest element in each row and column of the matrix B with
29 * elements B(i,j)=R(i)*A(i,j)*C(j) have an absolute value of at most
30 * the radix.
31 *
32 * R(i) and C(j) are restricted to be a power of the radix between
33 * SMLNUM = smallest safe number and BIGNUM = largest safe number. Use
34 * of these scaling factors is not guaranteed to reduce the condition
35 * number of A but works well in practice.
36 *
37 * This routine differs from DGEEQU by restricting the scaling factors
38 * to a power of the radix. Baring over- and underflow, scaling by
39 * these factors introduces no additional rounding errors. However, the
40 * scaled entries' magnitured are no longer approximately 1 but lie
41 * between sqrt(radix) and 1/sqrt(radix).
42 *
43 * Arguments
44 * =========
45 *
46 * M (input) INTEGER
47 * The number of rows of the matrix A. M >= 0.
48 *
49 * N (input) INTEGER
50 * The number of columns of the matrix A. N >= 0.
51 *
52 * KL (input) INTEGER
53 * The number of subdiagonals within the band of A. KL >= 0.
54 *
55 * KU (input) INTEGER
56 * The number of superdiagonals within the band of A. KU >= 0.
57 *
58 * AB (input) DOUBLE PRECISION array, dimension (LDAB,N)
59 * On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
60 * The j-th column of A is stored in the j-th column of the
61 * array AB as follows:
62 * AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
63 *
64 * LDAB (input) INTEGER
65 * The leading dimension of the array A. LDAB >= max(1,M).
66 *
67 * R (output) DOUBLE PRECISION array, dimension (M)
68 * If INFO = 0 or INFO > M, R contains the row scale factors
69 * for A.
70 *
71 * C (output) DOUBLE PRECISION array, dimension (N)
72 * If INFO = 0, C contains the column scale factors for A.
73 *
74 * ROWCND (output) DOUBLE PRECISION
75 * If INFO = 0 or INFO > M, ROWCND contains the ratio of the
76 * smallest R(i) to the largest R(i). If ROWCND >= 0.1 and
77 * AMAX is neither too large nor too small, it is not worth
78 * scaling by R.
79 *
80 * COLCND (output) DOUBLE PRECISION
81 * If INFO = 0, COLCND contains the ratio of the smallest
82 * C(i) to the largest C(i). If COLCND >= 0.1, it is not
83 * worth scaling by C.
84 *
85 * AMAX (output) DOUBLE PRECISION
86 * Absolute value of largest matrix element. If AMAX is very
87 * close to overflow or very close to underflow, the matrix
88 * should be scaled.
89 *
90 * INFO (output) INTEGER
91 * = 0: successful exit
92 * < 0: if INFO = -i, the i-th argument had an illegal value
93 * > 0: if INFO = i, and i is
94 * <= M: the i-th row of A is exactly zero
95 * > M: the (i-M)-th column of A is exactly zero
96 *
97 * =====================================================================
98 *
99 * .. Parameters ..
100 DOUBLE PRECISION ONE, ZERO
101 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
102 * ..
103 * .. Local Scalars ..
104 INTEGER I, J, KD
105 DOUBLE PRECISION BIGNUM, RCMAX, RCMIN, SMLNUM, RADIX, LOGRDX
106 * ..
107 * .. External Functions ..
108 DOUBLE PRECISION DLAMCH
109 EXTERNAL DLAMCH
110 * ..
111 * .. External Subroutines ..
112 EXTERNAL XERBLA
113 * ..
114 * .. Intrinsic Functions ..
115 INTRINSIC ABS, MAX, MIN, LOG
116 * ..
117 * .. Executable Statements ..
118 *
119 * Test the input parameters.
120 *
121 INFO = 0
122 IF( M.LT.0 ) THEN
123 INFO = -1
124 ELSE IF( N.LT.0 ) THEN
125 INFO = -2
126 ELSE IF( KL.LT.0 ) THEN
127 INFO = -3
128 ELSE IF( KU.LT.0 ) THEN
129 INFO = -4
130 ELSE IF( LDAB.LT.KL+KU+1 ) THEN
131 INFO = -6
132 END IF
133 IF( INFO.NE.0 ) THEN
134 CALL XERBLA( 'DGBEQUB', -INFO )
135 RETURN
136 END IF
137 *
138 * Quick return if possible.
139 *
140 IF( M.EQ.0 .OR. N.EQ.0 ) THEN
141 ROWCND = ONE
142 COLCND = ONE
143 AMAX = ZERO
144 RETURN
145 END IF
146 *
147 * Get machine constants. Assume SMLNUM is a power of the radix.
148 *
149 SMLNUM = DLAMCH( 'S' )
150 BIGNUM = ONE / SMLNUM
151 RADIX = DLAMCH( 'B' )
152 LOGRDX = LOG(RADIX)
153 *
154 * Compute row scale factors.
155 *
156 DO 10 I = 1, M
157 R( I ) = ZERO
158 10 CONTINUE
159 *
160 * Find the maximum element in each row.
161 *
162 KD = KU + 1
163 DO 30 J = 1, N
164 DO 20 I = MAX( J-KU, 1 ), MIN( J+KL, M )
165 R( I ) = MAX( R( I ), ABS( AB( KD+I-J, J ) ) )
166 20 CONTINUE
167 30 CONTINUE
168 DO I = 1, M
169 IF( R( I ).GT.ZERO ) THEN
170 R( I ) = RADIX**INT( LOG( R( I ) ) / LOGRDX )
171 END IF
172 END DO
173 *
174 * Find the maximum and minimum scale factors.
175 *
176 RCMIN = BIGNUM
177 RCMAX = ZERO
178 DO 40 I = 1, M
179 RCMAX = MAX( RCMAX, R( I ) )
180 RCMIN = MIN( RCMIN, R( I ) )
181 40 CONTINUE
182 AMAX = RCMAX
183 *
184 IF( RCMIN.EQ.ZERO ) THEN
185 *
186 * Find the first zero scale factor and return an error code.
187 *
188 DO 50 I = 1, M
189 IF( R( I ).EQ.ZERO ) THEN
190 INFO = I
191 RETURN
192 END IF
193 50 CONTINUE
194 ELSE
195 *
196 * Invert the scale factors.
197 *
198 DO 60 I = 1, M
199 R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM )
200 60 CONTINUE
201 *
202 * Compute ROWCND = min(R(I)) / max(R(I)).
203 *
204 ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
205 END IF
206 *
207 * Compute column scale factors.
208 *
209 DO 70 J = 1, N
210 C( J ) = ZERO
211 70 CONTINUE
212 *
213 * Find the maximum element in each column,
214 * assuming the row scaling computed above.
215 *
216 DO 90 J = 1, N
217 DO 80 I = MAX( J-KU, 1 ), MIN( J+KL, M )
218 C( J ) = MAX( C( J ), ABS( AB( KD+I-J, J ) )*R( I ) )
219 80 CONTINUE
220 IF( C( J ).GT.ZERO ) THEN
221 C( J ) = RADIX**INT( LOG( C( J ) ) / LOGRDX )
222 END IF
223 90 CONTINUE
224 *
225 * Find the maximum and minimum scale factors.
226 *
227 RCMIN = BIGNUM
228 RCMAX = ZERO
229 DO 100 J = 1, N
230 RCMIN = MIN( RCMIN, C( J ) )
231 RCMAX = MAX( RCMAX, C( J ) )
232 100 CONTINUE
233 *
234 IF( RCMIN.EQ.ZERO ) THEN
235 *
236 * Find the first zero scale factor and return an error code.
237 *
238 DO 110 J = 1, N
239 IF( C( J ).EQ.ZERO ) THEN
240 INFO = M + J
241 RETURN
242 END IF
243 110 CONTINUE
244 ELSE
245 *
246 * Invert the scale factors.
247 *
248 DO 120 J = 1, N
249 C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM )
250 120 CONTINUE
251 *
252 * Compute COLCND = min(C(J)) / max(C(J)).
253 *
254 COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
255 END IF
256 *
257 RETURN
258 *
259 * End of DGBEQUB
260 *
261 END
2 $ AMAX, INFO )
3 *
4 * -- LAPACK routine (version 3.2) --
5 * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
6 * -- Jason Riedy of Univ. of California Berkeley. --
7 * -- November 2008 --
8 *
9 * -- LAPACK is a software package provided by Univ. of Tennessee, --
10 * -- Univ. of California Berkeley and NAG Ltd. --
11 *
12 IMPLICIT NONE
13 * ..
14 * .. Scalar Arguments ..
15 INTEGER INFO, KL, KU, LDAB, M, N
16 DOUBLE PRECISION AMAX, COLCND, ROWCND
17 * ..
18 * .. Array Arguments ..
19 DOUBLE PRECISION AB( LDAB, * ), C( * ), R( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * DGBEQUB computes row and column scalings intended to equilibrate an
26 * M-by-N matrix A and reduce its condition number. R returns the row
27 * scale factors and C the column scale factors, chosen to try to make
28 * the largest element in each row and column of the matrix B with
29 * elements B(i,j)=R(i)*A(i,j)*C(j) have an absolute value of at most
30 * the radix.
31 *
32 * R(i) and C(j) are restricted to be a power of the radix between
33 * SMLNUM = smallest safe number and BIGNUM = largest safe number. Use
34 * of these scaling factors is not guaranteed to reduce the condition
35 * number of A but works well in practice.
36 *
37 * This routine differs from DGEEQU by restricting the scaling factors
38 * to a power of the radix. Baring over- and underflow, scaling by
39 * these factors introduces no additional rounding errors. However, the
40 * scaled entries' magnitured are no longer approximately 1 but lie
41 * between sqrt(radix) and 1/sqrt(radix).
42 *
43 * Arguments
44 * =========
45 *
46 * M (input) INTEGER
47 * The number of rows of the matrix A. M >= 0.
48 *
49 * N (input) INTEGER
50 * The number of columns of the matrix A. N >= 0.
51 *
52 * KL (input) INTEGER
53 * The number of subdiagonals within the band of A. KL >= 0.
54 *
55 * KU (input) INTEGER
56 * The number of superdiagonals within the band of A. KU >= 0.
57 *
58 * AB (input) DOUBLE PRECISION array, dimension (LDAB,N)
59 * On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
60 * The j-th column of A is stored in the j-th column of the
61 * array AB as follows:
62 * AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
63 *
64 * LDAB (input) INTEGER
65 * The leading dimension of the array A. LDAB >= max(1,M).
66 *
67 * R (output) DOUBLE PRECISION array, dimension (M)
68 * If INFO = 0 or INFO > M, R contains the row scale factors
69 * for A.
70 *
71 * C (output) DOUBLE PRECISION array, dimension (N)
72 * If INFO = 0, C contains the column scale factors for A.
73 *
74 * ROWCND (output) DOUBLE PRECISION
75 * If INFO = 0 or INFO > M, ROWCND contains the ratio of the
76 * smallest R(i) to the largest R(i). If ROWCND >= 0.1 and
77 * AMAX is neither too large nor too small, it is not worth
78 * scaling by R.
79 *
80 * COLCND (output) DOUBLE PRECISION
81 * If INFO = 0, COLCND contains the ratio of the smallest
82 * C(i) to the largest C(i). If COLCND >= 0.1, it is not
83 * worth scaling by C.
84 *
85 * AMAX (output) DOUBLE PRECISION
86 * Absolute value of largest matrix element. If AMAX is very
87 * close to overflow or very close to underflow, the matrix
88 * should be scaled.
89 *
90 * INFO (output) INTEGER
91 * = 0: successful exit
92 * < 0: if INFO = -i, the i-th argument had an illegal value
93 * > 0: if INFO = i, and i is
94 * <= M: the i-th row of A is exactly zero
95 * > M: the (i-M)-th column of A is exactly zero
96 *
97 * =====================================================================
98 *
99 * .. Parameters ..
100 DOUBLE PRECISION ONE, ZERO
101 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
102 * ..
103 * .. Local Scalars ..
104 INTEGER I, J, KD
105 DOUBLE PRECISION BIGNUM, RCMAX, RCMIN, SMLNUM, RADIX, LOGRDX
106 * ..
107 * .. External Functions ..
108 DOUBLE PRECISION DLAMCH
109 EXTERNAL DLAMCH
110 * ..
111 * .. External Subroutines ..
112 EXTERNAL XERBLA
113 * ..
114 * .. Intrinsic Functions ..
115 INTRINSIC ABS, MAX, MIN, LOG
116 * ..
117 * .. Executable Statements ..
118 *
119 * Test the input parameters.
120 *
121 INFO = 0
122 IF( M.LT.0 ) THEN
123 INFO = -1
124 ELSE IF( N.LT.0 ) THEN
125 INFO = -2
126 ELSE IF( KL.LT.0 ) THEN
127 INFO = -3
128 ELSE IF( KU.LT.0 ) THEN
129 INFO = -4
130 ELSE IF( LDAB.LT.KL+KU+1 ) THEN
131 INFO = -6
132 END IF
133 IF( INFO.NE.0 ) THEN
134 CALL XERBLA( 'DGBEQUB', -INFO )
135 RETURN
136 END IF
137 *
138 * Quick return if possible.
139 *
140 IF( M.EQ.0 .OR. N.EQ.0 ) THEN
141 ROWCND = ONE
142 COLCND = ONE
143 AMAX = ZERO
144 RETURN
145 END IF
146 *
147 * Get machine constants. Assume SMLNUM is a power of the radix.
148 *
149 SMLNUM = DLAMCH( 'S' )
150 BIGNUM = ONE / SMLNUM
151 RADIX = DLAMCH( 'B' )
152 LOGRDX = LOG(RADIX)
153 *
154 * Compute row scale factors.
155 *
156 DO 10 I = 1, M
157 R( I ) = ZERO
158 10 CONTINUE
159 *
160 * Find the maximum element in each row.
161 *
162 KD = KU + 1
163 DO 30 J = 1, N
164 DO 20 I = MAX( J-KU, 1 ), MIN( J+KL, M )
165 R( I ) = MAX( R( I ), ABS( AB( KD+I-J, J ) ) )
166 20 CONTINUE
167 30 CONTINUE
168 DO I = 1, M
169 IF( R( I ).GT.ZERO ) THEN
170 R( I ) = RADIX**INT( LOG( R( I ) ) / LOGRDX )
171 END IF
172 END DO
173 *
174 * Find the maximum and minimum scale factors.
175 *
176 RCMIN = BIGNUM
177 RCMAX = ZERO
178 DO 40 I = 1, M
179 RCMAX = MAX( RCMAX, R( I ) )
180 RCMIN = MIN( RCMIN, R( I ) )
181 40 CONTINUE
182 AMAX = RCMAX
183 *
184 IF( RCMIN.EQ.ZERO ) THEN
185 *
186 * Find the first zero scale factor and return an error code.
187 *
188 DO 50 I = 1, M
189 IF( R( I ).EQ.ZERO ) THEN
190 INFO = I
191 RETURN
192 END IF
193 50 CONTINUE
194 ELSE
195 *
196 * Invert the scale factors.
197 *
198 DO 60 I = 1, M
199 R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM )
200 60 CONTINUE
201 *
202 * Compute ROWCND = min(R(I)) / max(R(I)).
203 *
204 ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
205 END IF
206 *
207 * Compute column scale factors.
208 *
209 DO 70 J = 1, N
210 C( J ) = ZERO
211 70 CONTINUE
212 *
213 * Find the maximum element in each column,
214 * assuming the row scaling computed above.
215 *
216 DO 90 J = 1, N
217 DO 80 I = MAX( J-KU, 1 ), MIN( J+KL, M )
218 C( J ) = MAX( C( J ), ABS( AB( KD+I-J, J ) )*R( I ) )
219 80 CONTINUE
220 IF( C( J ).GT.ZERO ) THEN
221 C( J ) = RADIX**INT( LOG( C( J ) ) / LOGRDX )
222 END IF
223 90 CONTINUE
224 *
225 * Find the maximum and minimum scale factors.
226 *
227 RCMIN = BIGNUM
228 RCMAX = ZERO
229 DO 100 J = 1, N
230 RCMIN = MIN( RCMIN, C( J ) )
231 RCMAX = MAX( RCMAX, C( J ) )
232 100 CONTINUE
233 *
234 IF( RCMIN.EQ.ZERO ) THEN
235 *
236 * Find the first zero scale factor and return an error code.
237 *
238 DO 110 J = 1, N
239 IF( C( J ).EQ.ZERO ) THEN
240 INFO = M + J
241 RETURN
242 END IF
243 110 CONTINUE
244 ELSE
245 *
246 * Invert the scale factors.
247 *
248 DO 120 J = 1, N
249 C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM )
250 120 CONTINUE
251 *
252 * Compute COLCND = min(C(J)) / max(C(J)).
253 *
254 COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
255 END IF
256 *
257 RETURN
258 *
259 * End of DGBEQUB
260 *
261 END