1       SUBROUTINE ZGBEQU( M, N, KL, KU, AB, LDAB, R, C, ROWCND, COLCND,
  2      $                   AMAX, INFO )
  3 *
  4 *  -- LAPACK routine (version 3.2) --
  5 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  6 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  7 *     November 2006
  8 *
  9 *     .. Scalar Arguments ..
 10       INTEGER            INFO, KL, KU, LDAB, M, N
 11       DOUBLE PRECISION   AMAX, COLCND, ROWCND
 12 *     ..
 13 *     .. Array Arguments ..
 14       DOUBLE PRECISION   C( * ), R( * )
 15       COMPLEX*16         AB( LDAB, * )
 16 *     ..
 17 *
 18 *  Purpose
 19 *  =======
 20 *
 21 *  ZGBEQU computes row and column scalings intended to equilibrate an
 22 *  M-by-N band matrix A and reduce its condition number.  R returns the
 23 *  row scale factors and C the column scale factors, chosen to try to
 24 *  make the largest element in each row and column of the matrix B with
 25 *  elements B(i,j)=R(i)*A(i,j)*C(j) have absolute value 1.
 26 *
 27 *  R(i) and C(j) are restricted to be between SMLNUM = smallest safe
 28 *  number and BIGNUM = largest safe number.  Use of these scaling
 29 *  factors is not guaranteed to reduce the condition number of A but
 30 *  works well in practice.
 31 *
 32 *  Arguments
 33 *  =========
 34 *
 35 *  M       (input) INTEGER
 36 *          The number of rows of the matrix A.  M >= 0.
 37 *
 38 *  N       (input) INTEGER
 39 *          The number of columns of the matrix A.  N >= 0.
 40 *
 41 *  KL      (input) INTEGER
 42 *          The number of subdiagonals within the band of A.  KL >= 0.
 43 *
 44 *  KU      (input) INTEGER
 45 *          The number of superdiagonals within the band of A.  KU >= 0.
 46 *
 47 *  AB      (input) COMPLEX*16 array, dimension (LDAB,N)
 48 *          The band matrix A, stored in rows 1 to KL+KU+1.  The j-th
 49 *          column of A is stored in the j-th column of the array AB as
 50 *          follows:
 51 *          AB(ku+1+i-j,j) = A(i,j) for max(1,j-ku)<=i<=min(m,j+kl).
 52 *
 53 *  LDAB    (input) INTEGER
 54 *          The leading dimension of the array AB.  LDAB >= KL+KU+1.
 55 *
 56 *  R       (output) DOUBLE PRECISION array, dimension (M)
 57 *          If INFO = 0, or INFO > M, R contains the row scale factors
 58 *          for A.
 59 *
 60 *  C       (output) DOUBLE PRECISION array, dimension (N)
 61 *          If INFO = 0, C contains the column scale factors for A.
 62 *
 63 *  ROWCND  (output) DOUBLE PRECISION
 64 *          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
 65 *          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
 66 *          AMAX is neither too large nor too small, it is not worth
 67 *          scaling by R.
 68 *
 69 *  COLCND  (output) DOUBLE PRECISION
 70 *          If INFO = 0, COLCND contains the ratio of the smallest
 71 *          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
 72 *          worth scaling by C.
 73 *
 74 *  AMAX    (output) DOUBLE PRECISION
 75 *          Absolute value of largest matrix element.  If AMAX is very
 76 *          close to overflow or very close to underflow, the matrix
 77 *          should be scaled.
 78 *
 79 *  INFO    (output) INTEGER
 80 *          = 0:  successful exit
 81 *          < 0:  if INFO = -i, the i-th argument had an illegal value
 82 *          > 0:  if INFO = i, and i is
 83 *                <= M:  the i-th row of A is exactly zero
 84 *                >  M:  the (i-M)-th column of A is exactly zero
 85 *
 86 *  =====================================================================
 87 *
 88 *     .. Parameters ..
 89       DOUBLE PRECISION   ONE, ZERO
 90       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 91 *     ..
 92 *     .. Local Scalars ..
 93       INTEGER            I, J, KD
 94       DOUBLE PRECISION   BIGNUM, RCMAX, RCMIN, SMLNUM
 95       COMPLEX*16         ZDUM
 96 *     ..
 97 *     .. External Functions ..
 98       DOUBLE PRECISION   DLAMCH
 99       EXTERNAL           DLAMCH
100 *     ..
101 *     .. External Subroutines ..
102       EXTERNAL           XERBLA
103 *     ..
104 *     .. Intrinsic Functions ..
105       INTRINSIC          ABSDBLEDIMAGMAXMIN
106 *     ..
107 *     .. Statement Functions ..
108       DOUBLE PRECISION   CABS1
109 *     ..
110 *     .. Statement Function definitions ..
111       CABS1( ZDUM ) = ABSDBLE( ZDUM ) ) + ABSDIMAG( 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 / MINMAX( 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 / MINMAX( 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