LAPACK 3.3.1
Linear Algebra PACKage

dgeequb.f

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