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