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