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