SRC/zgeequb.f

   1 *> \brief \b ZGEEQUB
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZGEEQUB + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgeequb.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgeequb.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgeequb.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZGEEQUB( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
  22 *                           INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            INFO, LDA, M, N
  26 *       DOUBLE PRECISION   AMAX, COLCND, ROWCND
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       DOUBLE PRECISION   C( * ), R( * )
  30 *       COMPLEX*16         A( LDA, * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> ZGEEQUB computes row and column scalings intended to equilibrate an
  40 *> M-by-N matrix A and reduce its condition number.  R returns the row
  41 *> scale factors and C the column scale factors, chosen to try to make
  42 *> the largest element in each row and column of the matrix B with
  43 *> elements B(i,j)=R(i)*A(i,j)*C(j) have an absolute value of at most
  44 *> the radix.
  45 *>
  46 *> R(i) and C(j) are restricted to be a power of the radix between
  47 *> SMLNUM = smallest safe number and BIGNUM = largest safe number.  Use
  48 *> of these scaling factors is not guaranteed to reduce the condition
  49 *> number of A but works well in practice.
  50 *>
  51 *> This routine differs from ZGEEQU by restricting the scaling factors
  52 *> to a power of the radix.  Baring over- and underflow, scaling by
  53 *> these factors introduces no additional rounding errors.  However, the
  54 *> scaled entries' magnitured are no longer approximately 1 but lie
  55 *> between sqrt(radix) and 1/sqrt(radix).
  56 *> \endverbatim
  57 *
  58 *  Arguments:
  59 *  ==========
  60 *
  61 *> \param[in] M
  62 *> \verbatim
  63 *>          M is INTEGER
  64 *>          The number of rows of the matrix A.  M >= 0.
  65 *> \endverbatim
  66 *>
  67 *> \param[in] N
  68 *> \verbatim
  69 *>          N is INTEGER
  70 *>          The number of columns of the matrix A.  N >= 0.
  71 *> \endverbatim
  72 *>
  73 *> \param[in] A
  74 *> \verbatim
  75 *>          A is COMPLEX*16 array, dimension (LDA,N)
  76 *>          The M-by-N matrix whose equilibration factors are
  77 *>          to be computed.
  78 *> \endverbatim
  79 *>
  80 *> \param[in] LDA
  81 *> \verbatim
  82 *>          LDA is INTEGER
  83 *>          The leading dimension of the array A.  LDA >= max(1,M).
  84 *> \endverbatim
  85 *>
  86 *> \param[out] R
  87 *> \verbatim
  88 *>          R is DOUBLE PRECISION array, dimension (M)
  89 *>          If INFO = 0 or INFO > M, R contains the row scale factors
  90 *>          for A.
  91 *> \endverbatim
  92 *>
  93 *> \param[out] C
  94 *> \verbatim
  95 *>          C is DOUBLE PRECISION array, dimension (N)
  96 *>          If INFO = 0,  C contains the column scale factors for A.
  97 *> \endverbatim
  98 *>
  99 *> \param[out] ROWCND
 100 *> \verbatim
 101 *>          ROWCND is DOUBLE PRECISION
 102 *>          If INFO = 0 or INFO > M, ROWCND contains the ratio of the
 103 *>          smallest R(i) to the largest R(i).  If ROWCND >= 0.1 and
 104 *>          AMAX is neither too large nor too small, it is not worth
 105 *>          scaling by R.
 106 *> \endverbatim
 107 *>
 108 *> \param[out] COLCND
 109 *> \verbatim
 110 *>          COLCND is DOUBLE PRECISION
 111 *>          If INFO = 0, COLCND contains the ratio of the smallest
 112 *>          C(i) to the largest C(i).  If COLCND >= 0.1, it is not
 113 *>          worth scaling by C.
 114 *> \endverbatim
 115 *>
 116 *> \param[out] AMAX
 117 *> \verbatim
 118 *>          AMAX is DOUBLE PRECISION
 119 *>          Absolute value of largest matrix element.  If AMAX is very
 120 *>          close to overflow or very close to underflow, the matrix
 121 *>          should be scaled.
 122 *> \endverbatim
 123 *>
 124 *> \param[out] INFO
 125 *> \verbatim
 126 *>          INFO is INTEGER
 127 *>          = 0:  successful exit
 128 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
 129 *>          > 0:  if INFO = i,  and i is
 130 *>                <= M:  the i-th row of A is exactly zero
 131 *>                >  M:  the (i-M)-th column of A is exactly zero
 132 *> \endverbatim
 133 *
 134 *  Authors:
 135 *  ========
 136 *
 137 *> \author Univ. of Tennessee
 138 *> \author Univ. of California Berkeley
 139 *> \author Univ. of Colorado Denver
 140 *> \author NAG Ltd.
 141 *
 142 *> \date November 2011
 143 *
 144 *> \ingroup complex16GEcomputational
 145 *
 146 *  =====================================================================
 147       SUBROUTINE ZGEEQUB( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
 148      $                    INFO )
 149 *
 150 *  -- LAPACK computational routine (version 3.4.0) --
 151 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 152 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 153 *     November 2011
 154 *
 155 *     .. Scalar Arguments ..
 156       INTEGER            INFO, LDA, M, N
 157       DOUBLE PRECISION   AMAX, COLCND, ROWCND
 158 *     ..
 159 *     .. Array Arguments ..
 160       DOUBLE PRECISION   C( * ), R( * )
 161       COMPLEX*16         A( LDA, * )
 162 *     ..
 163 *
 164 *  =====================================================================
 165 *
 166 *     .. Parameters ..
 167       DOUBLE PRECISION   ONE, ZERO
 168       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 169 *     ..
 170 *     .. Local Scalars ..
 171       INTEGER            I, J
 172       DOUBLE PRECISION   BIGNUM, RCMAX, RCMIN, SMLNUM, RADIX, LOGRDX
 173       COMPLEX*16         ZDUM
 174 *     ..
 175 *     .. External Functions ..
 176       DOUBLE PRECISION   DLAMCH
 177       EXTERNAL           DLAMCH
 178 *     ..
 179 *     .. External Subroutines ..
 180       EXTERNAL           XERBLA
 181 *     ..
 182 *     .. Intrinsic Functions ..
 183       INTRINSIC          ABS, MAX, MIN, LOG, DBLE, DIMAG
 184 *     ..
 185 *     .. Statement Functions ..
 186       DOUBLE PRECISION   CABS1
 187 *     ..
 188 *     .. Statement Function definitions ..
 189       CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
 190 *     ..
 191 *     .. Executable Statements ..
 192 *
 193 *     Test the input parameters.
 194 *
 195       INFO = 0
 196       IF( M.LT.0 ) THEN
 197          INFO = -1
 198       ELSE IF( N.LT.0 ) THEN
 199          INFO = -2
 200       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
 201          INFO = -4
 202       END IF
 203       IF( INFO.NE.0 ) THEN
 204          CALL XERBLA( 'ZGEEQUB', -INFO )
 205          RETURN
 206       END IF
 207 *
 208 *     Quick return if possible.
 209 *
 210       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
 211          ROWCND = ONE
 212          COLCND = ONE
 213          AMAX = ZERO
 214          RETURN
 215       END IF
 216 *
 217 *     Get machine constants.  Assume SMLNUM is a power of the radix.
 218 *
 219       SMLNUM = DLAMCH( 'S' )
 220       BIGNUM = ONE / SMLNUM
 221       RADIX = DLAMCH( 'B' )
 222       LOGRDX = LOG( RADIX )
 223 *
 224 *     Compute row scale factors.
 225 *
 226       DO 10 I = 1, M
 227          R( I ) = ZERO
 228    10 CONTINUE
 229 *
 230 *     Find the maximum element in each row.
 231 *
 232       DO 30 J = 1, N
 233          DO 20 I = 1, M
 234             R( I ) = MAX( R( I ), CABS1( A( I, J ) ) )
 235    20    CONTINUE
 236    30 CONTINUE
 237       DO I = 1, M
 238          IF( R( I ).GT.ZERO ) THEN
 239             R( I ) = RADIX**INT( LOG(R( I ) ) / LOGRDX )
 240          END IF
 241       END DO
 242 *
 243 *     Find the maximum and minimum scale factors.
 244 *
 245       RCMIN = BIGNUM
 246       RCMAX = ZERO
 247       DO 40 I = 1, M
 248          RCMAX = MAX( RCMAX, R( I ) )
 249          RCMIN = MIN( RCMIN, R( I ) )
 250    40 CONTINUE
 251       AMAX = RCMAX
 252 *
 253       IF( RCMIN.EQ.ZERO ) THEN
 254 *
 255 *        Find the first zero scale factor and return an error code.
 256 *
 257          DO 50 I = 1, M
 258             IF( R( I ).EQ.ZERO ) THEN
 259                INFO = I
 260                RETURN
 261             END IF
 262    50    CONTINUE
 263       ELSE
 264 *
 265 *        Invert the scale factors.
 266 *
 267          DO 60 I = 1, M
 268             R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM )
 269    60    CONTINUE
 270 *
 271 *        Compute ROWCND = min(R(I)) / max(R(I)).
 272 *
 273          ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
 274       END IF
 275 *
 276 *     Compute column scale factors.
 277 *
 278       DO 70 J = 1, N
 279          C( J ) = ZERO
 280    70 CONTINUE
 281 *
 282 *     Find the maximum element in each column,
 283 *     assuming the row scaling computed above.
 284 *
 285       DO 90 J = 1, N
 286          DO 80 I = 1, M
 287             C( J ) = MAX( C( J ), CABS1( A( I, J ) )*R( I ) )
 288    80    CONTINUE
 289          IF( C( J ).GT.ZERO ) THEN
 290             C( J ) = RADIX**INT( LOG( C( J ) ) / LOGRDX )
 291          END IF
 292    90 CONTINUE
 293 *
 294 *     Find the maximum and minimum scale factors.
 295 *
 296       RCMIN = BIGNUM
 297       RCMAX = ZERO
 298       DO 100 J = 1, N
 299          RCMIN = MIN( RCMIN, C( J ) )
 300          RCMAX = MAX( RCMAX, C( J ) )
 301   100 CONTINUE
 302 *
 303       IF( RCMIN.EQ.ZERO ) THEN
 304 *
 305 *        Find the first zero scale factor and return an error code.
 306 *
 307          DO 110 J = 1, N
 308             IF( C( J ).EQ.ZERO ) THEN
 309                INFO = M + J
 310                RETURN
 311             END IF
 312   110    CONTINUE
 313       ELSE
 314 *
 315 *        Invert the scale factors.
 316 *
 317          DO 120 J = 1, N
 318             C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM )
 319   120    CONTINUE
 320 *
 321 *        Compute COLCND = min(C(J)) / max(C(J)).
 322 *
 323          COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
 324       END IF
 325 *
 326       RETURN
 327 *
 328 *     End of ZGEEQUB
 329 *
 330       END