1*> \brief \b DGBEQUB
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DGBEQUB + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgbequb.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgbequb.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgbequb.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DGBEQUB( M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND,
22*                           AMAX, INFO )
23*
24*       .. Scalar Arguments ..
25*       INTEGER            INFO, KL, KU, LDAB, M, N
26*       DOUBLE PRECISION   AMAX, COLCND, ROWCND
27*       ..
28*       .. Array Arguments ..
29*       DOUBLE PRECISION   AB( LDAB, * ), C( * ), R( * )
30*       ..
31*
32*
33*> \par Purpose:
34*  =============
35*>
36*> \verbatim
37*>
38*> DGBEQUB computes row and column scalings intended to equilibrate an
39*> M-by-N matrix A and reduce its condition number.  R returns the row
40*> scale factors and C the column scale factors, chosen to try to make
41*> the largest element in each row and column of the matrix B with
42*> elements B(i,j)=R(i)*A(i,j)*C(j) have an absolute value of at most
43*> the radix.
44*>
45*> R(i) and C(j) are restricted to be a power of the radix between
46*> SMLNUM = smallest safe number and BIGNUM = largest safe number.  Use
47*> of these scaling factors is not guaranteed to reduce the condition
48*> number of A but works well in practice.
49*>
50*> This routine differs from DGEEQU by restricting the scaling factors
51*> to a power of the radix.  Baring over- and underflow, scaling by
52*> these factors introduces no additional rounding errors.  However, the
53*> scaled entries' magnitured are no longer approximately 1 but lie
54*> between sqrt(radix) and 1/sqrt(radix).
55*> \endverbatim
56*
57*  Arguments:
58*  ==========
59*
60*> \param[in] M
61*> \verbatim
62*>          M is INTEGER
63*>          The number of rows of the matrix A.  M >= 0.
64*> \endverbatim
65*>
66*> \param[in] N
67*> \verbatim
68*>          N is INTEGER
69*>          The number of columns of the matrix A.  N >= 0.
70*> \endverbatim
71*>
72*> \param[in] KL
73*> \verbatim
74*>          KL is INTEGER
75*>          The number of subdiagonals within the band of A.  KL >= 0.
76*> \endverbatim
77*>
78*> \param[in] KU
79*> \verbatim
80*>          KU is INTEGER
81*>          The number of superdiagonals within the band of A.  KU >= 0.
82*> \endverbatim
83*>
84*> \param[in] AB
85*> \verbatim
86*>          AB is DOUBLE PRECISION array, dimension (LDAB,N)
87*>          On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
88*>          The j-th column of A is stored in the j-th column of the
89*>          array AB as follows:
90*>          AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
91*> \endverbatim
92*>
93*> \param[in] LDAB
94*> \verbatim
95*>          LDAB is INTEGER
96*>          The leading dimension of the array A.  LDAB >= max(1,M).
97*> \endverbatim
98*>
99*> \param[out] R
100*> \verbatim
101*>          R is DOUBLE PRECISION array, dimension (M)
102*>          If INFO = 0 or INFO > M, R contains the row scale factors
103*>          for A.
104*> \endverbatim
105*>
106*> \param[out] C
107*> \verbatim
108*>          C is DOUBLE PRECISION array, dimension (N)
109*>          If INFO = 0,  C contains the column scale factors for A.
110*> \endverbatim
111*>
112*> \param[out] ROWCND
113*> \verbatim
114*>          ROWCND is DOUBLE PRECISION
115*>          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
116*>          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
117*>          AMAX is neither too large nor too small, it is not worth
118*>          scaling by R.
119*> \endverbatim
120*>
121*> \param[out] COLCND
122*> \verbatim
123*>          COLCND is DOUBLE PRECISION
124*>          If INFO = 0, COLCND contains the ratio of the smallest
125*>          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
126*>          worth scaling by C.
127*> \endverbatim
128*>
129*> \param[out] AMAX
130*> \verbatim
131*>          AMAX is DOUBLE PRECISION
132*>          Absolute value of largest matrix element.  If AMAX is very
133*>          close to overflow or very close to underflow, the matrix
134*>          should be scaled.
135*> \endverbatim
136*>
137*> \param[out] INFO
138*> \verbatim
139*>          INFO is INTEGER
140*>          = 0:  successful exit
141*>          < 0:  if INFO = -i, the i-th argument had an illegal value
142*>          > 0:  if INFO = i,  and i is
143*>                <= M:  the i-th row of A is exactly zero
144*>                >  M:  the (i-M)-th column of A is exactly zero
145*> \endverbatim
146*
147*  Authors:
148*  ========
149*
150*> \author Univ. of Tennessee
151*> \author Univ. of California Berkeley
152*> \author Univ. of Colorado Denver
153*> \author NAG Ltd.
154*
155*> \date November 2011
156*
157*> \ingroup doubleGBcomputational
158*
159*  =====================================================================
160      SUBROUTINE DGBEQUB( M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND,
161     $                    AMAX, INFO )
162*
163*  -- LAPACK computational routine (version 3.4.0) --
164*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
165*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
166*     November 2011
167*
168*     .. Scalar Arguments ..
169      INTEGER            INFO, KL, KU, LDAB, M, N
170      DOUBLE PRECISION   AMAX, COLCND, ROWCND
171*     ..
172*     .. Array Arguments ..
173      DOUBLE PRECISION   AB( LDAB, * ), C( * ), R( * )
174*     ..
175*
176*  =====================================================================
177*
178*     .. Parameters ..
179      DOUBLE PRECISION   ONE, ZERO
180      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
181*     ..
182*     .. Local Scalars ..
183      INTEGER            I, J, KD
184      DOUBLE PRECISION   BIGNUM, RCMAX, RCMIN, SMLNUM, RADIX, LOGRDX
185*     ..
186*     .. External Functions ..
187      DOUBLE PRECISION   DLAMCH
188      EXTERNAL           DLAMCH
189*     ..
190*     .. External Subroutines ..
191      EXTERNAL           XERBLA
192*     ..
193*     .. Intrinsic Functions ..
194      INTRINSIC          ABS, MAX, MIN, LOG
195*     ..
196*     .. Executable Statements ..
197*
198*     Test the input parameters.
199*
200      INFO = 0
201      IF( M.LT.0 ) THEN
202         INFO = -1
203      ELSE IF( N.LT.0 ) THEN
204         INFO = -2
205      ELSE IF( KL.LT.0 ) THEN
206         INFO = -3
207      ELSE IF( KU.LT.0 ) THEN
208         INFO = -4
209      ELSE IF( LDAB.LT.KL+KU+1 ) THEN
210         INFO = -6
211      END IF
212      IF( INFO.NE.0 ) THEN
213         CALL XERBLA( 'DGBEQUB', -INFO )
214         RETURN
215      END IF
216*
217*     Quick return if possible.
218*
219      IF( M.EQ.0 .OR. N.EQ.0 ) THEN
220         ROWCND = ONE
221         COLCND = ONE
222         AMAX = ZERO
223         RETURN
224      END IF
225*
226*     Get machine constants.  Assume SMLNUM is a power of the radix.
227*
228      SMLNUM = DLAMCH( 'S' )
229      BIGNUM = ONE / SMLNUM
230      RADIX = DLAMCH( 'B' )
231      LOGRDX = LOG(RADIX)
232*
233*     Compute row scale factors.
234*
235      DO 10 I = 1, M
236         R( I ) = ZERO
237   10 CONTINUE
238*
239*     Find the maximum element in each row.
240*
241      KD = KU + 1
242      DO 30 J = 1, N
243         DO 20 I = MAX( J-KU, 1 ), MIN( J+KL, M )
244            R( I ) = MAX( R( I ), ABS( AB( KD+I-J, J ) ) )
245   20    CONTINUE
246   30 CONTINUE
247      DO I = 1, M
248         IF( R( I ).GT.ZERO ) THEN
249            R( I ) = RADIX**INT( LOG( R( I ) ) / LOGRDX )
250         END IF
251      END DO
252*
253*     Find the maximum and minimum scale factors.
254*
255      RCMIN = BIGNUM
256      RCMAX = ZERO
257      DO 40 I = 1, M
258         RCMAX = MAX( RCMAX, R( I ) )
259         RCMIN = MIN( RCMIN, R( I ) )
260   40 CONTINUE
261      AMAX = RCMAX
262*
263      IF( RCMIN.EQ.ZERO ) THEN
264*
265*        Find the first zero scale factor and return an error code.
266*
267         DO 50 I = 1, M
268            IF( R( I ).EQ.ZERO ) THEN
269               INFO = I
270               RETURN
271            END IF
272   50    CONTINUE
273      ELSE
274*
275*        Invert the scale factors.
276*
277         DO 60 I = 1, M
278            R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM )
279   60    CONTINUE
280*
281*        Compute ROWCND = min(R(I)) / max(R(I)).
282*
283         ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
284      END IF
285*
286*     Compute column scale factors.
287*
288      DO 70 J = 1, N
289         C( J ) = ZERO
290   70 CONTINUE
291*
292*     Find the maximum element in each column,
293*     assuming the row scaling computed above.
294*
295      DO 90 J = 1, N
296         DO 80 I = MAX( J-KU, 1 ), MIN( J+KL, M )
297            C( J ) = MAX( C( J ), ABS( AB( KD+I-J, J ) )*R( I ) )
298   80    CONTINUE
299         IF( C( J ).GT.ZERO ) THEN
300            C( J ) = RADIX**INT( LOG( C( J ) ) / LOGRDX )
301         END IF
302   90 CONTINUE
303*
304*     Find the maximum and minimum scale factors.
305*
306      RCMIN = BIGNUM
307      RCMAX = ZERO
308      DO 100 J = 1, N
309         RCMIN = MIN( RCMIN, C( J ) )
310         RCMAX = MAX( RCMAX, C( J ) )
311  100 CONTINUE
312*
313      IF( RCMIN.EQ.ZERO ) THEN
314*
315*        Find the first zero scale factor and return an error code.
316*
317         DO 110 J = 1, N
318            IF( C( J ).EQ.ZERO ) THEN
319               INFO = M + J
320               RETURN
321            END IF
322  110    CONTINUE
323      ELSE
324*
325*        Invert the scale factors.
326*
327         DO 120 J = 1, N
328            C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM )
329  120    CONTINUE
330*
331*        Compute COLCND = min(C(J)) / max(C(J)).
332*
333         COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
334      END IF
335*
336      RETURN
337*
338*     End of DGBEQUB
339*
340      END
341