001:       SUBROUTINE DGBEQUB( M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND,
002:      $                    AMAX, INFO )
003: *
004: *     -- LAPACK routine (version 3.2)                                 --
005: *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
006: *     -- Jason Riedy of Univ. of California Berkeley.                 --
007: *     -- November 2008                                                --
008: *
009: *     -- LAPACK is a software package provided by Univ. of Tennessee, --
010: *     -- Univ. of California Berkeley and NAG Ltd.                    --
011: *
012:       IMPLICIT NONE
013: *     ..
014: *     .. Scalar Arguments ..
015:       INTEGER            INFO, KL, KU, LDAB, M, N
016:       DOUBLE PRECISION   AMAX, COLCND, ROWCND
017: *     ..
018: *     .. Array Arguments ..
019:       DOUBLE PRECISION   AB( LDAB, * ), C( * ), R( * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  DGBEQUB computes row and column scalings intended to equilibrate an
026: *  M-by-N matrix A and reduce its condition number.  R returns the row
027: *  scale factors and C the column scale factors, chosen to try to make
028: *  the largest element in each row and column of the matrix B with
029: *  elements B(i,j)=R(i)*A(i,j)*C(j) have an absolute value of at most
030: *  the radix.
031: *
032: *  R(i) and C(j) are restricted to be a power of the radix between
033: *  SMLNUM = smallest safe number and BIGNUM = largest safe number.  Use
034: *  of these scaling factors is not guaranteed to reduce the condition
035: *  number of A but works well in practice.
036: *
037: *  This routine differs from DGEEQU by restricting the scaling factors
038: *  to a power of the radix.  Baring over- and underflow, scaling by
039: *  these factors introduces no additional rounding errors.  However, the
040: *  scaled entries' magnitured are no longer approximately 1 but lie
041: *  between sqrt(radix) and 1/sqrt(radix).
042: *
043: *  Arguments
044: *  =========
045: *
046: *  M       (input) INTEGER
047: *          The number of rows of the matrix A.  M >= 0.
048: *
049: *  N       (input) INTEGER
050: *          The number of columns of the matrix A.  N >= 0.
051: *
052: *  KL      (input) INTEGER
053: *          The number of subdiagonals within the band of A.  KL >= 0.
054: *
055: *  KU      (input) INTEGER
056: *          The number of superdiagonals within the band of A.  KU >= 0.
057: *
058: *  AB      (input) DOUBLE PRECISION array, dimension (LDAB,N)
059: *          On entry, the matrix A in band storage, in rows 1 to KL+KU+1.
060: *          The j-th column of A is stored in the j-th column of the
061: *          array AB as follows:
062: *          AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl)
063: *
064: *  LDAB    (input) INTEGER
065: *          The leading dimension of the array A.  LDAB >= max(1,M).
066: *
067: *  R       (output) DOUBLE PRECISION array, dimension (M)
068: *          If INFO = 0 or INFO > M, R contains the row scale factors
069: *          for A.
070: *
071: *  C       (output) DOUBLE PRECISION array, dimension (N)
072: *          If INFO = 0,  C contains the column scale factors for A.
073: *
074: *  ROWCND  (output) DOUBLE PRECISION
075: *          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
076: *          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
077: *          AMAX is neither too large nor too small, it is not worth
078: *          scaling by R.
079: *
080: *  COLCND  (output) DOUBLE PRECISION
081: *          If INFO = 0, COLCND contains the ratio of the smallest
082: *          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
083: *          worth scaling by C.
084: *
085: *  AMAX    (output) DOUBLE PRECISION
086: *          Absolute value of largest matrix element.  If AMAX is very
087: *          close to overflow or very close to underflow, the matrix
088: *          should be scaled.
089: *
090: *  INFO    (output) INTEGER
091: *          = 0:  successful exit
092: *          < 0:  if INFO = -i, the i-th argument had an illegal value
093: *          > 0:  if INFO = i,  and i is
094: *                <= M:  the i-th row of A is exactly zero
095: *                >  M:  the (i-M)-th column of A is exactly zero
096: *
097: *  =====================================================================
098: *
099: *     .. 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
262: