001:       SUBROUTINE DGBEQU( M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND,
002:      $                   AMAX, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       INTEGER            INFO, KL, KU, LDAB, M, N
010:       DOUBLE PRECISION   AMAX, COLCND, ROWCND
011: *     ..
012: *     .. Array Arguments ..
013:       DOUBLE PRECISION   AB( LDAB, * ), C( * ), R( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  DGBEQU computes row and column scalings intended to equilibrate an
020: *  M-by-N band matrix A and reduce its condition number.  R returns the
021: *  row scale factors and C the column scale factors, chosen to try to
022: *  make the largest element in each row and column of the matrix B with
023: *  elements B(i,j)=R(i)*A(i,j)*C(j) have absolute value 1.
024: *
025: *  R(i) and C(j) are restricted to be between SMLNUM = smallest safe
026: *  number and BIGNUM = largest safe number.  Use of these scaling
027: *  factors is not guaranteed to reduce the condition number of A but
028: *  works well in practice.
029: *
030: *  Arguments
031: *  =========
032: *
033: *  M       (input) INTEGER
034: *          The number of rows of the matrix A.  M >= 0.
035: *
036: *  N       (input) INTEGER
037: *          The number of columns of the matrix A.  N >= 0.
038: *
039: *  KL      (input) INTEGER
040: *          The number of subdiagonals within the band of A.  KL >= 0.
041: *
042: *  KU      (input) INTEGER
043: *          The number of superdiagonals within the band of A.  KU >= 0.
044: *
045: *  AB      (input) DOUBLE PRECISION array, dimension (LDAB,N)
046: *          The band matrix A, stored in rows 1 to KL+KU+1.  The j-th
047: *          column of A is stored in the j-th column of the array AB as
048: *          follows:
049: *          AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
050: *
051: *  LDAB    (input) INTEGER
052: *          The leading dimension of the array AB.  LDAB >= KL+KU+1.
053: *
054: *  R       (output) DOUBLE PRECISION array, dimension (M)
055: *          If INFO = 0, or INFO > M, R contains the row scale factors
056: *          for A.
057: *
058: *  C       (output) DOUBLE PRECISION array, dimension (N)
059: *          If INFO = 0, C contains the column scale factors for A.
060: *
061: *  ROWCND  (output) DOUBLE PRECISION
062: *          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
063: *          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
064: *          AMAX is neither too large nor too small, it is not worth
065: *          scaling by R.
066: *
067: *  COLCND  (output) DOUBLE PRECISION
068: *          If INFO = 0, COLCND contains the ratio of the smallest
069: *          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
070: *          worth scaling by C.
071: *
072: *  AMAX    (output) DOUBLE PRECISION
073: *          Absolute value of largest matrix element.  If AMAX is very
074: *          close to overflow or very close to underflow, the matrix
075: *          should be scaled.
076: *
077: *  INFO    (output) INTEGER
078: *          = 0:  successful exit
079: *          < 0:  if INFO = -i, the i-th argument had an illegal value
080: *          > 0:  if INFO = i, and i is
081: *                <= M:  the i-th row of A is exactly zero
082: *                >  M:  the (i-M)-th column of A is exactly zero
083: *
084: *  =====================================================================
085: *
086: *     .. Parameters ..
087:       DOUBLE PRECISION   ONE, ZERO
088:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
089: *     ..
090: *     .. Local Scalars ..
091:       INTEGER            I, J, KD
092:       DOUBLE PRECISION   BIGNUM, RCMAX, RCMIN, SMLNUM
093: *     ..
094: *     .. External Functions ..
095:       DOUBLE PRECISION   DLAMCH
096:       EXTERNAL           DLAMCH
097: *     ..
098: *     .. External Subroutines ..
099:       EXTERNAL           XERBLA
100: *     ..
101: *     .. Intrinsic Functions ..
102:       INTRINSIC          ABS, MAX, MIN
103: *     ..
104: *     .. Executable Statements ..
105: *
106: *     Test the input parameters
107: *
108:       INFO = 0
109:       IF( M.LT.0 ) THEN
110:          INFO = -1
111:       ELSE IF( N.LT.0 ) THEN
112:          INFO = -2
113:       ELSE IF( KL.LT.0 ) THEN
114:          INFO = -3
115:       ELSE IF( KU.LT.0 ) THEN
116:          INFO = -4
117:       ELSE IF( LDAB.LT.KL+KU+1 ) THEN
118:          INFO = -6
119:       END IF
120:       IF( INFO.NE.0 ) THEN
121:          CALL XERBLA( 'DGBEQU', -INFO )
122:          RETURN
123:       END IF
124: *
125: *     Quick return if possible
126: *
127:       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
128:          ROWCND = ONE
129:          COLCND = ONE
130:          AMAX = ZERO
131:          RETURN
132:       END IF
133: *
134: *     Get machine constants.
135: *
136:       SMLNUM = DLAMCH( 'S' )
137:       BIGNUM = ONE / SMLNUM
138: *
139: *     Compute row scale factors.
140: *
141:       DO 10 I = 1, M
142:          R( I ) = ZERO
143:    10 CONTINUE
144: *
145: *     Find the maximum element in each row.
146: *
147:       KD = KU + 1
148:       DO 30 J = 1, N
149:          DO 20 I = MAX( J-KU, 1 ), MIN( J+KL, M )
150:             R( I ) = MAX( R( I ), ABS( AB( KD+I-J, J ) ) )
151:    20    CONTINUE
152:    30 CONTINUE
153: *
154: *     Find the maximum and minimum scale factors.
155: *
156:       RCMIN = BIGNUM
157:       RCMAX = ZERO
158:       DO 40 I = 1, M
159:          RCMAX = MAX( RCMAX, R( I ) )
160:          RCMIN = MIN( RCMIN, R( I ) )
161:    40 CONTINUE
162:       AMAX = RCMAX
163: *
164:       IF( RCMIN.EQ.ZERO ) THEN
165: *
166: *        Find the first zero scale factor and return an error code.
167: *
168:          DO 50 I = 1, M
169:             IF( R( I ).EQ.ZERO ) THEN
170:                INFO = I
171:                RETURN
172:             END IF
173:    50    CONTINUE
174:       ELSE
175: *
176: *        Invert the scale factors.
177: *
178:          DO 60 I = 1, M
179:             R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM )
180:    60    CONTINUE
181: *
182: *        Compute ROWCND = min(R(I)) / max(R(I))
183: *
184:          ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
185:       END IF
186: *
187: *     Compute column scale factors
188: *
189:       DO 70 J = 1, N
190:          C( J ) = ZERO
191:    70 CONTINUE
192: *
193: *     Find the maximum element in each column,
194: *     assuming the row scaling computed above.
195: *
196:       KD = KU + 1
197:       DO 90 J = 1, N
198:          DO 80 I = MAX( J-KU, 1 ), MIN( J+KL, M )
199:             C( J ) = MAX( C( J ), ABS( AB( KD+I-J, J ) )*R( I ) )
200:    80    CONTINUE
201:    90 CONTINUE
202: *
203: *     Find the maximum and minimum scale factors.
204: *
205:       RCMIN = BIGNUM
206:       RCMAX = ZERO
207:       DO 100 J = 1, N
208:          RCMIN = MIN( RCMIN, C( J ) )
209:          RCMAX = MAX( RCMAX, C( J ) )
210:   100 CONTINUE
211: *
212:       IF( RCMIN.EQ.ZERO ) THEN
213: *
214: *        Find the first zero scale factor and return an error code.
215: *
216:          DO 110 J = 1, N
217:             IF( C( J ).EQ.ZERO ) THEN
218:                INFO = M + J
219:                RETURN
220:             END IF
221:   110    CONTINUE
222:       ELSE
223: *
224: *        Invert the scale factors.
225: *
226:          DO 120 J = 1, N
227:             C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM )
228:   120    CONTINUE
229: *
230: *        Compute COLCND = min(C(J)) / max(C(J))
231: *
232:          COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
233:       END IF
234: *
235:       RETURN
236: *
237: *     End of DGBEQU
238: *
239:       END
240: