1 SUBROUTINE DGEEQUB( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
2 $ 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, LDA, M, N
16 DOUBLE PRECISION AMAX, COLCND, ROWCND
17 * ..
18 * .. Array Arguments ..
19 DOUBLE PRECISION A( LDA, * ), C( * ), R( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * DGEEQUB 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 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
53 * The M-by-N matrix whose equilibration factors are
54 * to be computed.
55 *
56 * LDA (input) INTEGER
57 * The leading dimension of the array A. LDA >= max(1,M).
58 *
59 * R (output) DOUBLE PRECISION array, dimension (M)
60 * If INFO = 0 or INFO > M, R contains the row scale factors
61 * for A.
62 *
63 * C (output) DOUBLE PRECISION array, dimension (N)
64 * If INFO = 0, C contains the column scale factors for A.
65 *
66 * ROWCND (output) DOUBLE PRECISION
67 * If INFO = 0 or INFO > M, ROWCND contains the ratio of the
68 * smallest R(i) to the largest R(i). If ROWCND >= 0.1 and
69 * AMAX is neither too large nor too small, it is not worth
70 * scaling by R.
71 *
72 * COLCND (output) DOUBLE PRECISION
73 * If INFO = 0, COLCND contains the ratio of the smallest
74 * C(i) to the largest C(i). If COLCND >= 0.1, it is not
75 * worth scaling by C.
76 *
77 * AMAX (output) DOUBLE PRECISION
78 * Absolute value of largest matrix element. If AMAX is very
79 * close to overflow or very close to underflow, the matrix
80 * should be scaled.
81 *
82 * INFO (output) INTEGER
83 * = 0: successful exit
84 * < 0: if INFO = -i, the i-th argument had an illegal value
85 * > 0: if INFO = i, and i is
86 * <= M: the i-th row of A is exactly zero
87 * > M: the (i-M)-th column of A is exactly zero
88 *
89 * =====================================================================
90 *
91 * .. Parameters ..
92 DOUBLE PRECISION ONE, ZERO
93 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
94 * ..
95 * .. Local Scalars ..
96 INTEGER I, J
97 DOUBLE PRECISION BIGNUM, RCMAX, RCMIN, SMLNUM, RADIX, LOGRDX
98 * ..
99 * .. External Functions ..
100 DOUBLE PRECISION DLAMCH
101 EXTERNAL DLAMCH
102 * ..
103 * .. External Subroutines ..
104 EXTERNAL XERBLA
105 * ..
106 * .. Intrinsic Functions ..
107 INTRINSIC ABS, MAX, MIN, LOG
108 * ..
109 * .. Executable Statements ..
110 *
111 * Test the input parameters.
112 *
113 INFO = 0
114 IF( M.LT.0 ) THEN
115 INFO = -1
116 ELSE IF( N.LT.0 ) THEN
117 INFO = -2
118 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
119 INFO = -4
120 END IF
121 IF( INFO.NE.0 ) THEN
122 CALL XERBLA( 'DGEEQUB', -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. Assume SMLNUM is a power of the radix.
136 *
137 SMLNUM = DLAMCH( 'S' )
138 BIGNUM = ONE / SMLNUM
139 RADIX = DLAMCH( 'B' )
140 LOGRDX = LOG( RADIX )
141 *
142 * Compute row scale factors.
143 *
144 DO 10 I = 1, M
145 R( I ) = ZERO
146 10 CONTINUE
147 *
148 * Find the maximum element in each row.
149 *
150 DO 30 J = 1, N
151 DO 20 I = 1, M
152 R( I ) = MAX( R( I ), ABS( A( I, J ) ) )
153 20 CONTINUE
154 30 CONTINUE
155 DO I = 1, M
156 IF( R( I ).GT.ZERO ) THEN
157 R( I ) = RADIX**INT( LOG( R( I ) ) / LOGRDX )
158 END IF
159 END DO
160 *
161 * Find the maximum and minimum scale factors.
162 *
163 RCMIN = BIGNUM
164 RCMAX = ZERO
165 DO 40 I = 1, M
166 RCMAX = MAX( RCMAX, R( I ) )
167 RCMIN = MIN( RCMIN, R( I ) )
168 40 CONTINUE
169 AMAX = RCMAX
170 *
171 IF( RCMIN.EQ.ZERO ) THEN
172 *
173 * Find the first zero scale factor and return an error code.
174 *
175 DO 50 I = 1, M
176 IF( R( I ).EQ.ZERO ) THEN
177 INFO = I
178 RETURN
179 END IF
180 50 CONTINUE
181 ELSE
182 *
183 * Invert the scale factors.
184 *
185 DO 60 I = 1, M
186 R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM )
187 60 CONTINUE
188 *
189 * Compute ROWCND = min(R(I)) / max(R(I)).
190 *
191 ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
192 END IF
193 *
194 * Compute column scale factors
195 *
196 DO 70 J = 1, N
197 C( J ) = ZERO
198 70 CONTINUE
199 *
200 * Find the maximum element in each column,
201 * assuming the row scaling computed above.
202 *
203 DO 90 J = 1, N
204 DO 80 I = 1, M
205 C( J ) = MAX( C( J ), ABS( A( I, J ) )*R( I ) )
206 80 CONTINUE
207 IF( C( J ).GT.ZERO ) THEN
208 C( J ) = RADIX**INT( LOG( C( J ) ) / LOGRDX )
209 END IF
210 90 CONTINUE
211 *
212 * Find the maximum and minimum scale factors.
213 *
214 RCMIN = BIGNUM
215 RCMAX = ZERO
216 DO 100 J = 1, N
217 RCMIN = MIN( RCMIN, C( J ) )
218 RCMAX = MAX( RCMAX, C( J ) )
219 100 CONTINUE
220 *
221 IF( RCMIN.EQ.ZERO ) THEN
222 *
223 * Find the first zero scale factor and return an error code.
224 *
225 DO 110 J = 1, N
226 IF( C( J ).EQ.ZERO ) THEN
227 INFO = M + J
228 RETURN
229 END IF
230 110 CONTINUE
231 ELSE
232 *
233 * Invert the scale factors.
234 *
235 DO 120 J = 1, N
236 C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM )
237 120 CONTINUE
238 *
239 * Compute COLCND = min(C(J)) / max(C(J)).
240 *
241 COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
242 END IF
243 *
244 RETURN
245 *
246 * End of DGEEQUB
247 *
248 END
2 $ 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, LDA, M, N
16 DOUBLE PRECISION AMAX, COLCND, ROWCND
17 * ..
18 * .. Array Arguments ..
19 DOUBLE PRECISION A( LDA, * ), C( * ), R( * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * DGEEQUB 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 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
53 * The M-by-N matrix whose equilibration factors are
54 * to be computed.
55 *
56 * LDA (input) INTEGER
57 * The leading dimension of the array A. LDA >= max(1,M).
58 *
59 * R (output) DOUBLE PRECISION array, dimension (M)
60 * If INFO = 0 or INFO > M, R contains the row scale factors
61 * for A.
62 *
63 * C (output) DOUBLE PRECISION array, dimension (N)
64 * If INFO = 0, C contains the column scale factors for A.
65 *
66 * ROWCND (output) DOUBLE PRECISION
67 * If INFO = 0 or INFO > M, ROWCND contains the ratio of the
68 * smallest R(i) to the largest R(i). If ROWCND >= 0.1 and
69 * AMAX is neither too large nor too small, it is not worth
70 * scaling by R.
71 *
72 * COLCND (output) DOUBLE PRECISION
73 * If INFO = 0, COLCND contains the ratio of the smallest
74 * C(i) to the largest C(i). If COLCND >= 0.1, it is not
75 * worth scaling by C.
76 *
77 * AMAX (output) DOUBLE PRECISION
78 * Absolute value of largest matrix element. If AMAX is very
79 * close to overflow or very close to underflow, the matrix
80 * should be scaled.
81 *
82 * INFO (output) INTEGER
83 * = 0: successful exit
84 * < 0: if INFO = -i, the i-th argument had an illegal value
85 * > 0: if INFO = i, and i is
86 * <= M: the i-th row of A is exactly zero
87 * > M: the (i-M)-th column of A is exactly zero
88 *
89 * =====================================================================
90 *
91 * .. Parameters ..
92 DOUBLE PRECISION ONE, ZERO
93 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
94 * ..
95 * .. Local Scalars ..
96 INTEGER I, J
97 DOUBLE PRECISION BIGNUM, RCMAX, RCMIN, SMLNUM, RADIX, LOGRDX
98 * ..
99 * .. External Functions ..
100 DOUBLE PRECISION DLAMCH
101 EXTERNAL DLAMCH
102 * ..
103 * .. External Subroutines ..
104 EXTERNAL XERBLA
105 * ..
106 * .. Intrinsic Functions ..
107 INTRINSIC ABS, MAX, MIN, LOG
108 * ..
109 * .. Executable Statements ..
110 *
111 * Test the input parameters.
112 *
113 INFO = 0
114 IF( M.LT.0 ) THEN
115 INFO = -1
116 ELSE IF( N.LT.0 ) THEN
117 INFO = -2
118 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
119 INFO = -4
120 END IF
121 IF( INFO.NE.0 ) THEN
122 CALL XERBLA( 'DGEEQUB', -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. Assume SMLNUM is a power of the radix.
136 *
137 SMLNUM = DLAMCH( 'S' )
138 BIGNUM = ONE / SMLNUM
139 RADIX = DLAMCH( 'B' )
140 LOGRDX = LOG( RADIX )
141 *
142 * Compute row scale factors.
143 *
144 DO 10 I = 1, M
145 R( I ) = ZERO
146 10 CONTINUE
147 *
148 * Find the maximum element in each row.
149 *
150 DO 30 J = 1, N
151 DO 20 I = 1, M
152 R( I ) = MAX( R( I ), ABS( A( I, J ) ) )
153 20 CONTINUE
154 30 CONTINUE
155 DO I = 1, M
156 IF( R( I ).GT.ZERO ) THEN
157 R( I ) = RADIX**INT( LOG( R( I ) ) / LOGRDX )
158 END IF
159 END DO
160 *
161 * Find the maximum and minimum scale factors.
162 *
163 RCMIN = BIGNUM
164 RCMAX = ZERO
165 DO 40 I = 1, M
166 RCMAX = MAX( RCMAX, R( I ) )
167 RCMIN = MIN( RCMIN, R( I ) )
168 40 CONTINUE
169 AMAX = RCMAX
170 *
171 IF( RCMIN.EQ.ZERO ) THEN
172 *
173 * Find the first zero scale factor and return an error code.
174 *
175 DO 50 I = 1, M
176 IF( R( I ).EQ.ZERO ) THEN
177 INFO = I
178 RETURN
179 END IF
180 50 CONTINUE
181 ELSE
182 *
183 * Invert the scale factors.
184 *
185 DO 60 I = 1, M
186 R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM )
187 60 CONTINUE
188 *
189 * Compute ROWCND = min(R(I)) / max(R(I)).
190 *
191 ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
192 END IF
193 *
194 * Compute column scale factors
195 *
196 DO 70 J = 1, N
197 C( J ) = ZERO
198 70 CONTINUE
199 *
200 * Find the maximum element in each column,
201 * assuming the row scaling computed above.
202 *
203 DO 90 J = 1, N
204 DO 80 I = 1, M
205 C( J ) = MAX( C( J ), ABS( A( I, J ) )*R( I ) )
206 80 CONTINUE
207 IF( C( J ).GT.ZERO ) THEN
208 C( J ) = RADIX**INT( LOG( C( J ) ) / LOGRDX )
209 END IF
210 90 CONTINUE
211 *
212 * Find the maximum and minimum scale factors.
213 *
214 RCMIN = BIGNUM
215 RCMAX = ZERO
216 DO 100 J = 1, N
217 RCMIN = MIN( RCMIN, C( J ) )
218 RCMAX = MAX( RCMAX, C( J ) )
219 100 CONTINUE
220 *
221 IF( RCMIN.EQ.ZERO ) THEN
222 *
223 * Find the first zero scale factor and return an error code.
224 *
225 DO 110 J = 1, N
226 IF( C( J ).EQ.ZERO ) THEN
227 INFO = M + J
228 RETURN
229 END IF
230 110 CONTINUE
231 ELSE
232 *
233 * Invert the scale factors.
234 *
235 DO 120 J = 1, N
236 C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM )
237 120 CONTINUE
238 *
239 * Compute COLCND = min(C(J)) / max(C(J)).
240 *
241 COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
242 END IF
243 *
244 RETURN
245 *
246 * End of DGEEQUB
247 *
248 END