SRC/zgbequb.f

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