1       SUBROUTINE DGEEQU( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
  2      $                   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, LDA, M, N
 11       DOUBLE PRECISION   AMAX, COLCND, ROWCND
 12 *     ..
 13 *     .. Array Arguments ..
 14       DOUBLE PRECISION   A( LDA, * ), C( * ), R( * )
 15 *     ..
 16 *
 17 *  Purpose
 18 *  =======
 19 *
 20 *  DGEEQU computes row and column scalings intended to equilibrate an
 21 *  M-by-N matrix A and reduce its condition number.  R returns the row
 22 *  scale factors and C the column scale factors, chosen to try to make
 23 *  the largest element in each row and column of the matrix B with
 24 *  elements B(i,j)=R(i)*A(i,j)*C(j) have absolute value 1.
 25 *
 26 *  R(i) and C(j) are restricted to be between SMLNUM = smallest safe
 27 *  number and BIGNUM = largest safe number.  Use of these scaling
 28 *  factors is not guaranteed to reduce the condition number of A but
 29 *  works well in practice.
 30 *
 31 *  Arguments
 32 *  =========
 33 *
 34 *  M       (input) INTEGER
 35 *          The number of rows of the matrix A.  M >= 0.
 36 *
 37 *  N       (input) INTEGER
 38 *          The number of columns of the matrix A.  N >= 0.
 39 *
 40 *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
 41 *          The M-by-N matrix whose equilibration factors are
 42 *          to be computed.
 43 *
 44 *  LDA     (input) INTEGER
 45 *          The leading dimension of the array A.  LDA >= max(1,M).
 46 *
 47 *  R       (output) DOUBLE PRECISION array, dimension (M)
 48 *          If INFO = 0 or INFO > M, R contains the row scale factors
 49 *          for A.
 50 *
 51 *  C       (output) DOUBLE PRECISION array, dimension (N)
 52 *          If INFO = 0,  C contains the column scale factors for A.
 53 *
 54 *  ROWCND  (output) DOUBLE PRECISION
 55 *          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
 56 *          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
 57 *          AMAX is neither too large nor too small, it is not worth
 58 *          scaling by R.
 59 *
 60 *  COLCND  (output) DOUBLE PRECISION
 61 *          If INFO = 0, COLCND contains the ratio of the smallest
 62 *          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
 63 *          worth scaling by C.
 64 *
 65 *  AMAX    (output) DOUBLE PRECISION
 66 *          Absolute value of largest matrix element.  If AMAX is very
 67 *          close to overflow or very close to underflow, the matrix
 68 *          should be scaled.
 69 *
 70 *  INFO    (output) INTEGER
 71 *          = 0:  successful exit
 72 *          < 0:  if INFO = -i, the i-th argument had an illegal value
 73 *          > 0:  if INFO = i,  and i is
 74 *                <= M:  the i-th row of A is exactly zero
 75 *                >  M:  the (i-M)-th column of A is exactly zero
 76 *
 77 *  =====================================================================
 78 *
 79 *     .. Parameters ..
 80       DOUBLE PRECISION   ONE, ZERO
 81       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 82 *     ..
 83 *     .. Local Scalars ..
 84       INTEGER            I, J
 85       DOUBLE PRECISION   BIGNUM, RCMAX, RCMIN, SMLNUM
 86 *     ..
 87 *     .. External Functions ..
 88       DOUBLE PRECISION   DLAMCH
 89       EXTERNAL           DLAMCH
 90 *     ..
 91 *     .. External Subroutines ..
 92       EXTERNAL           XERBLA
 93 *     ..
 94 *     .. Intrinsic Functions ..
 95       INTRINSIC          ABSMAXMIN
 96 *     ..
 97 *     .. Executable Statements ..
 98 *
 99 *     Test the input parameters.
100 *
101       INFO = 0
102       IF( M.LT.0 ) THEN
103          INFO = -1
104       ELSE IF( N.LT.0 ) THEN
105          INFO = -2
106       ELSE IF( LDA.LT.MAX1, M ) ) THEN
107          INFO = -4
108       END IF
109       IF( INFO.NE.0 ) THEN
110          CALL XERBLA( 'DGEEQU'-INFO )
111          RETURN
112       END IF
113 *
114 *     Quick return if possible
115 *
116       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
117          ROWCND = ONE
118          COLCND = ONE
119          AMAX = ZERO
120          RETURN
121       END IF
122 *
123 *     Get machine constants.
124 *
125       SMLNUM = DLAMCH( 'S' )
126       BIGNUM = ONE / SMLNUM
127 *
128 *     Compute row scale factors.
129 *
130       DO 10 I = 1, M
131          R( I ) = ZERO
132    10 CONTINUE
133 *
134 *     Find the maximum element in each row.
135 *
136       DO 30 J = 1, N
137          DO 20 I = 1, M
138             R( I ) = MAX( R( I ), ABS( A( I, J ) ) )
139    20    CONTINUE
140    30 CONTINUE
141 *
142 *     Find the maximum and minimum scale factors.
143 *
144       RCMIN = BIGNUM
145       RCMAX = ZERO
146       DO 40 I = 1, M
147          RCMAX = MAX( RCMAX, R( I ) )
148          RCMIN = MIN( RCMIN, R( I ) )
149    40 CONTINUE
150       AMAX = RCMAX
151 *
152       IF( RCMIN.EQ.ZERO ) THEN
153 *
154 *        Find the first zero scale factor and return an error code.
155 *
156          DO 50 I = 1, M
157             IF( R( I ).EQ.ZERO ) THEN
158                INFO = I
159                RETURN
160             END IF
161    50    CONTINUE
162       ELSE
163 *
164 *        Invert the scale factors.
165 *
166          DO 60 I = 1, M
167             R( I ) = ONE / MINMAX( R( I ), SMLNUM ), BIGNUM )
168    60    CONTINUE
169 *
170 *        Compute ROWCND = min(R(I)) / max(R(I))
171 *
172          ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
173       END IF
174 *
175 *     Compute column scale factors
176 *
177       DO 70 J = 1, N
178          C( J ) = ZERO
179    70 CONTINUE
180 *
181 *     Find the maximum element in each column,
182 *     assuming the row scaling computed above.
183 *
184       DO 90 J = 1, N
185          DO 80 I = 1, M
186             C( J ) = MAX( C( J ), ABS( A( I, J ) )*R( I ) )
187    80    CONTINUE
188    90 CONTINUE
189 *
190 *     Find the maximum and minimum scale factors.
191 *
192       RCMIN = BIGNUM
193       RCMAX = ZERO
194       DO 100 J = 1, N
195          RCMIN = MIN( RCMIN, C( J ) )
196          RCMAX = MAX( RCMAX, C( J ) )
197   100 CONTINUE
198 *
199       IF( RCMIN.EQ.ZERO ) THEN
200 *
201 *        Find the first zero scale factor and return an error code.
202 *
203          DO 110 J = 1, N
204             IF( C( J ).EQ.ZERO ) THEN
205                INFO = M + J
206                RETURN
207             END IF
208   110    CONTINUE
209       ELSE
210 *
211 *        Invert the scale factors.
212 *
213          DO 120 J = 1, N
214             C( J ) = ONE / MINMAX( C( J ), SMLNUM ), BIGNUM )
215   120    CONTINUE
216 *
217 *        Compute COLCND = min(C(J)) / max(C(J))
218 *
219          COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
220       END IF
221 *
222       RETURN
223 *
224 *     End of DGEEQU
225 *
226       END