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