LAPACK 3.3.1
Linear Algebra PACKage

zgeequb.f

Go to the documentation of this file.
00001       SUBROUTINE ZGEEQUB( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
00002      $                    INFO )
00003 *
00004 *     -- LAPACK routine (version 3.2)                                 --
00005 *     -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
00006 *     -- Jason Riedy of Univ. of California Berkeley.                 --
00007 *     -- November 2008                                                --
00008 *
00009 *     -- LAPACK is a software package provided by Univ. of Tennessee, --
00010 *     -- Univ. of California Berkeley and NAG Ltd.                    --
00011 *
00012       IMPLICIT NONE
00013 *     ..
00014 *     .. Scalar Arguments ..
00015       INTEGER            INFO, LDA, M, N
00016       DOUBLE PRECISION   AMAX, COLCND, ROWCND
00017 *     ..
00018 *     .. Array Arguments ..
00019       DOUBLE PRECISION   C( * ), R( * )
00020       COMPLEX*16         A( LDA, * )
00021 *     ..
00022 *
00023 *  Purpose
00024 *  =======
00025 *
00026 *  ZGEEQUB computes row and column scalings intended to equilibrate an
00027 *  M-by-N matrix A and reduce its condition number.  R returns the row
00028 *  scale factors and C the column scale factors, chosen to try to make
00029 *  the largest element in each row and column of the matrix B with
00030 *  elements B(i,j)=R(i)*A(i,j)*C(j) have an absolute value of at most
00031 *  the radix.
00032 *
00033 *  R(i) and C(j) are restricted to be a power of the radix between
00034 *  SMLNUM = smallest safe number and BIGNUM = largest safe number.  Use
00035 *  of these scaling factors is not guaranteed to reduce the condition
00036 *  number of A but works well in practice.
00037 *
00038 *  This routine differs from ZGEEQU by restricting the scaling factors
00039 *  to a power of the radix.  Baring over- and underflow, scaling by
00040 *  these factors introduces no additional rounding errors.  However, the
00041 *  scaled entries' magnitured are no longer approximately 1 but lie
00042 *  between sqrt(radix) and 1/sqrt(radix).
00043 *
00044 *  Arguments
00045 *  =========
00046 *
00047 *  M       (input) INTEGER
00048 *          The number of rows of the matrix A.  M >= 0.
00049 *
00050 *  N       (input) INTEGER
00051 *          The number of columns of the matrix A.  N >= 0.
00052 *
00053 *  A       (input) COMPLEX*16 array, dimension (LDA,N)
00054 *          The M-by-N matrix whose equilibration factors are
00055 *          to be computed.
00056 *
00057 *  LDA     (input) INTEGER
00058 *          The leading dimension of the array A.  LDA >= max(1,M).
00059 *
00060 *  R       (output) DOUBLE PRECISION array, dimension (M)
00061 *          If INFO = 0 or INFO > M, R contains the row scale factors
00062 *          for A.
00063 *
00064 *  C       (output) DOUBLE PRECISION array, dimension (N)
00065 *          If INFO = 0,  C contains the column scale factors for A.
00066 *
00067 *  ROWCND  (output) DOUBLE PRECISION
00068 *          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
00069 *          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
00070 *          AMAX is neither too large nor too small, it is not worth
00071 *          scaling by R.
00072 *
00073 *  COLCND  (output) DOUBLE PRECISION
00074 *          If INFO = 0, COLCND contains the ratio of the smallest
00075 *          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
00076 *          worth scaling by C.
00077 *
00078 *  AMAX    (output) DOUBLE PRECISION
00079 *          Absolute value of largest matrix element.  If AMAX is very
00080 *          close to overflow or very close to underflow, the matrix
00081 *          should be scaled.
00082 *
00083 *  INFO    (output) INTEGER
00084 *          = 0:  successful exit
00085 *          < 0:  if INFO = -i, the i-th argument had an illegal value
00086 *          > 0:  if INFO = i,  and i is
00087 *                <= M:  the i-th row of A is exactly zero
00088 *                >  M:  the (i-M)-th column of A is exactly zero
00089 *
00090 *  =====================================================================
00091 *
00092 *     .. Parameters ..
00093       DOUBLE PRECISION   ONE, ZERO
00094       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
00095 *     ..
00096 *     .. Local Scalars ..
00097       INTEGER            I, J
00098       DOUBLE PRECISION   BIGNUM, RCMAX, RCMIN, SMLNUM, RADIX, LOGRDX
00099       COMPLEX*16         ZDUM
00100 *     ..
00101 *     .. External Functions ..
00102       DOUBLE PRECISION   DLAMCH
00103       EXTERNAL           DLAMCH
00104 *     ..
00105 *     .. External Subroutines ..
00106       EXTERNAL           XERBLA
00107 *     ..
00108 *     .. Intrinsic Functions ..
00109       INTRINSIC          ABS, MAX, MIN, LOG, REAL, DIMAG
00110 *     ..
00111 *     .. Statement Functions ..
00112       DOUBLE PRECISION   CABS1
00113 *     ..
00114 *     .. Statement Function definitions ..
00115       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
00116 *     ..
00117 *     .. Executable Statements ..
00118 *
00119 *     Test the input parameters.
00120 *
00121       INFO = 0
00122       IF( M.LT.0 ) THEN
00123          INFO = -1
00124       ELSE IF( N.LT.0 ) THEN
00125          INFO = -2
00126       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
00127          INFO = -4
00128       END IF
00129       IF( INFO.NE.0 ) THEN
00130          CALL XERBLA( 'ZGEEQUB', -INFO )
00131          RETURN
00132       END IF
00133 *
00134 *     Quick return if possible.
00135 *
00136       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
00137          ROWCND = ONE
00138          COLCND = ONE
00139          AMAX = ZERO
00140          RETURN
00141       END IF
00142 *
00143 *     Get machine constants.  Assume SMLNUM is a power of the radix.
00144 *
00145       SMLNUM = DLAMCH( 'S' )
00146       BIGNUM = ONE / SMLNUM
00147       RADIX = DLAMCH( 'B' )
00148       LOGRDX = LOG( RADIX )
00149 *
00150 *     Compute row scale factors.
00151 *
00152       DO 10 I = 1, M
00153          R( I ) = ZERO
00154    10 CONTINUE
00155 *
00156 *     Find the maximum element in each row.
00157 *
00158       DO 30 J = 1, N
00159          DO 20 I = 1, M
00160             R( I ) = MAX( R( I ), CABS1( A( I, J ) ) )
00161    20    CONTINUE
00162    30 CONTINUE
00163       DO I = 1, M
00164          IF( R( I ).GT.ZERO ) THEN
00165             R( I ) = RADIX**INT( LOG(R( I ) ) / LOGRDX )
00166          END IF
00167       END DO
00168 *
00169 *     Find the maximum and minimum scale factors.
00170 *
00171       RCMIN = BIGNUM
00172       RCMAX = ZERO
00173       DO 40 I = 1, M
00174          RCMAX = MAX( RCMAX, R( I ) )
00175          RCMIN = MIN( RCMIN, R( I ) )
00176    40 CONTINUE
00177       AMAX = RCMAX
00178 *
00179       IF( RCMIN.EQ.ZERO ) THEN
00180 *
00181 *        Find the first zero scale factor and return an error code.
00182 *
00183          DO 50 I = 1, M
00184             IF( R( I ).EQ.ZERO ) THEN
00185                INFO = I
00186                RETURN
00187             END IF
00188    50    CONTINUE
00189       ELSE
00190 *
00191 *        Invert the scale factors.
00192 *
00193          DO 60 I = 1, M
00194             R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM )
00195    60    CONTINUE
00196 *
00197 *        Compute ROWCND = min(R(I)) / max(R(I)).
00198 *
00199          ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
00200       END IF
00201 *
00202 *     Compute column scale factors.
00203 *
00204       DO 70 J = 1, N
00205          C( J ) = ZERO
00206    70 CONTINUE
00207 *
00208 *     Find the maximum element in each column,
00209 *     assuming the row scaling computed above.
00210 *
00211       DO 90 J = 1, N
00212          DO 80 I = 1, M
00213             C( J ) = MAX( C( J ), CABS1( A( I, J ) )*R( I ) )
00214    80    CONTINUE
00215          IF( C( J ).GT.ZERO ) THEN
00216             C( J ) = RADIX**INT( LOG( C( J ) ) / LOGRDX )
00217          END IF
00218    90 CONTINUE
00219 *
00220 *     Find the maximum and minimum scale factors.
00221 *
00222       RCMIN = BIGNUM
00223       RCMAX = ZERO
00224       DO 100 J = 1, N
00225          RCMIN = MIN( RCMIN, C( J ) )
00226          RCMAX = MAX( RCMAX, C( J ) )
00227   100 CONTINUE
00228 *
00229       IF( RCMIN.EQ.ZERO ) THEN
00230 *
00231 *        Find the first zero scale factor and return an error code.
00232 *
00233          DO 110 J = 1, N
00234             IF( C( J ).EQ.ZERO ) THEN
00235                INFO = M + J
00236                RETURN
00237             END IF
00238   110    CONTINUE
00239       ELSE
00240 *
00241 *        Invert the scale factors.
00242 *
00243          DO 120 J = 1, N
00244             C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM )
00245   120    CONTINUE
00246 *
00247 *        Compute COLCND = min(C(J)) / max(C(J)).
00248 *
00249          COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
00250       END IF
00251 *
00252       RETURN
00253 *
00254 *     End of ZGEEQUB
00255 *
00256       END
 All Files Functions