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