001:       SUBROUTINE DGEEQU( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
002:      $                   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, LDA, M, N
010:       DOUBLE PRECISION   AMAX, COLCND, ROWCND
011: *     ..
012: *     .. Array Arguments ..
013:       DOUBLE PRECISION   A( LDA, * ), C( * ), R( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  DGEEQU computes row and column scalings intended to equilibrate an
020: *  M-by-N matrix A and reduce its condition number.  R returns the row
021: *  scale factors and C the column scale factors, chosen to try to make
022: *  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: *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)
040: *          The M-by-N matrix whose equilibration factors are
041: *          to be computed.
042: *
043: *  LDA     (input) INTEGER
044: *          The leading dimension of the array A.  LDA >= max(1,M).
045: *
046: *  R       (output) DOUBLE PRECISION array, dimension (M)
047: *          If INFO = 0 or INFO > M, R contains the row scale factors
048: *          for A.
049: *
050: *  C       (output) DOUBLE PRECISION array, dimension (N)
051: *          If INFO = 0,  C contains the column scale factors for A.
052: *
053: *  ROWCND  (output) DOUBLE PRECISION
054: *          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
055: *          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
056: *          AMAX is neither too large nor too small, it is not worth
057: *          scaling by R.
058: *
059: *  COLCND  (output) DOUBLE PRECISION
060: *          If INFO = 0, COLCND contains the ratio of the smallest
061: *          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
062: *          worth scaling by C.
063: *
064: *  AMAX    (output) DOUBLE PRECISION
065: *          Absolute value of largest matrix element.  If AMAX is very
066: *          close to overflow or very close to underflow, the matrix
067: *          should be scaled.
068: *
069: *  INFO    (output) INTEGER
070: *          = 0:  successful exit
071: *          < 0:  if INFO = -i, the i-th argument had an illegal value
072: *          > 0:  if INFO = i,  and i is
073: *                <= M:  the i-th row of A is exactly zero
074: *                >  M:  the (i-M)-th column of A is exactly zero
075: *
076: *  =====================================================================
077: *
078: *     .. Parameters ..
079:       DOUBLE PRECISION   ONE, ZERO
080:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
081: *     ..
082: *     .. Local Scalars ..
083:       INTEGER            I, J
084:       DOUBLE PRECISION   BIGNUM, RCMAX, RCMIN, SMLNUM
085: *     ..
086: *     .. External Functions ..
087:       DOUBLE PRECISION   DLAMCH
088:       EXTERNAL           DLAMCH
089: *     ..
090: *     .. External Subroutines ..
091:       EXTERNAL           XERBLA
092: *     ..
093: *     .. Intrinsic Functions ..
094:       INTRINSIC          ABS, MAX, MIN
095: *     ..
096: *     .. Executable Statements ..
097: *
098: *     Test the input parameters.
099: *
100:       INFO = 0
101:       IF( M.LT.0 ) THEN
102:          INFO = -1
103:       ELSE IF( N.LT.0 ) THEN
104:          INFO = -2
105:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
106:          INFO = -4
107:       END IF
108:       IF( INFO.NE.0 ) THEN
109:          CALL XERBLA( 'DGEEQU', -INFO )
110:          RETURN
111:       END IF
112: *
113: *     Quick return if possible
114: *
115:       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
116:          ROWCND = ONE
117:          COLCND = ONE
118:          AMAX = ZERO
119:          RETURN
120:       END IF
121: *
122: *     Get machine constants.
123: *
124:       SMLNUM = DLAMCH( 'S' )
125:       BIGNUM = ONE / SMLNUM
126: *
127: *     Compute row scale factors.
128: *
129:       DO 10 I = 1, M
130:          R( I ) = ZERO
131:    10 CONTINUE
132: *
133: *     Find the maximum element in each row.
134: *
135:       DO 30 J = 1, N
136:          DO 20 I = 1, M
137:             R( I ) = MAX( R( I ), ABS( A( I, J ) ) )
138:    20    CONTINUE
139:    30 CONTINUE
140: *
141: *     Find the maximum and minimum scale factors.
142: *
143:       RCMIN = BIGNUM
144:       RCMAX = ZERO
145:       DO 40 I = 1, M
146:          RCMAX = MAX( RCMAX, R( I ) )
147:          RCMIN = MIN( RCMIN, R( I ) )
148:    40 CONTINUE
149:       AMAX = RCMAX
150: *
151:       IF( RCMIN.EQ.ZERO ) THEN
152: *
153: *        Find the first zero scale factor and return an error code.
154: *
155:          DO 50 I = 1, M
156:             IF( R( I ).EQ.ZERO ) THEN
157:                INFO = I
158:                RETURN
159:             END IF
160:    50    CONTINUE
161:       ELSE
162: *
163: *        Invert the scale factors.
164: *
165:          DO 60 I = 1, M
166:             R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM )
167:    60    CONTINUE
168: *
169: *        Compute ROWCND = min(R(I)) / max(R(I))
170: *
171:          ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
172:       END IF
173: *
174: *     Compute column scale factors
175: *
176:       DO 70 J = 1, N
177:          C( J ) = ZERO
178:    70 CONTINUE
179: *
180: *     Find the maximum element in each column,
181: *     assuming the row scaling computed above.
182: *
183:       DO 90 J = 1, N
184:          DO 80 I = 1, M
185:             C( J ) = MAX( C( J ), ABS( A( I, J ) )*R( I ) )
186:    80    CONTINUE
187:    90 CONTINUE
188: *
189: *     Find the maximum and minimum scale factors.
190: *
191:       RCMIN = BIGNUM
192:       RCMAX = ZERO
193:       DO 100 J = 1, N
194:          RCMIN = MIN( RCMIN, C( J ) )
195:          RCMAX = MAX( RCMAX, C( J ) )
196:   100 CONTINUE
197: *
198:       IF( RCMIN.EQ.ZERO ) THEN
199: *
200: *        Find the first zero scale factor and return an error code.
201: *
202:          DO 110 J = 1, N
203:             IF( C( J ).EQ.ZERO ) THEN
204:                INFO = M + J
205:                RETURN
206:             END IF
207:   110    CONTINUE
208:       ELSE
209: *
210: *        Invert the scale factors.
211: *
212:          DO 120 J = 1, N
213:             C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM )
214:   120    CONTINUE
215: *
216: *        Compute COLCND = min(C(J)) / max(C(J))
217: *
218:          COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
219:       END IF
220: *
221:       RETURN
222: *
223: *     End of DGEEQU
224: *
225:       END
226: