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