SRC/dgbequb.f

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