SRC/dgebal.f

   1 *> \brief \b DGEBAL
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DGEBAL + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgebal.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgebal.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgebal.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       CHARACTER          JOB
  25 *       INTEGER            IHI, ILO, INFO, LDA, N
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       DOUBLE PRECISION   A( LDA, * ), SCALE( * )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> DGEBAL balances a general real matrix A.  This involves, first,
  38 *> permuting A by a similarity transformation to isolate eigenvalues
  39 *> in the first 1 to ILO-1 and last IHI+1 to N elements on the
  40 *> diagonal; and second, applying a diagonal similarity transformation
  41 *> to rows and columns ILO to IHI to make the rows and columns as
  42 *> close in norm as possible.  Both steps are optional.
  43 *>
  44 *> Balancing may reduce the 1-norm of the matrix, and improve the
  45 *> accuracy of the computed eigenvalues and/or eigenvectors.
  46 *> \endverbatim
  47 *
  48 *  Arguments:
  49 *  ==========
  50 *
  51 *> \param[in] JOB
  52 *> \verbatim
  53 *>          JOB is CHARACTER*1
  54 *>          Specifies the operations to be performed on A:
  55 *>          = 'N':  none:  simply set ILO = 1, IHI = N, SCALE(I) = 1.0
  56 *>                  for i = 1,...,N;
  57 *>          = 'P':  permute only;
  58 *>          = 'S':  scale only;
  59 *>          = 'B':  both permute and scale.
  60 *> \endverbatim
  61 *>
  62 *> \param[in] N
  63 *> \verbatim
  64 *>          N is INTEGER
  65 *>          The order of the matrix A.  N >= 0.
  66 *> \endverbatim
  67 *>
  68 *> \param[in,out] A
  69 *> \verbatim
  70 *>          A is DOUBLE array, dimension (LDA,N)
  71 *>          On entry, the input matrix A.
  72 *>          On exit,  A is overwritten by the balanced matrix.
  73 *>          If JOB = 'N', A is not referenced.
  74 *>          See Further Details.
  75 *> \endverbatim
  76 *>
  77 *> \param[in] LDA
  78 *> \verbatim
  79 *>          LDA is INTEGER
  80 *>          The leading dimension of the array A.  LDA >= max(1,N).
  81 *> \endverbatim
  82 *>
  83 *> \param[out] ILO
  84 *> \verbatim
  85 *>          ILO is INTEGER
  86 *> \endverbatim
  87 *> \param[out] IHI
  88 *> \verbatim
  89 *>          IHI is INTEGER
  90 *>          ILO and IHI are set to integers such that on exit
  91 *>          A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
  92 *>          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
  93 *> \endverbatim
  94 *>
  95 *> \param[out] SCALE
  96 *> \verbatim
  97 *>          SCALE is DOUBLE array, dimension (N)
  98 *>          Details of the permutations and scaling factors applied to
  99 *>          A.  If P(j) is the index of the row and column interchanged
 100 *>          with row and column j and D(j) is the scaling factor
 101 *>          applied to row and column j, then
 102 *>          SCALE(j) = P(j)    for j = 1,...,ILO-1
 103 *>                   = D(j)    for j = ILO,...,IHI
 104 *>                   = P(j)    for j = IHI+1,...,N.
 105 *>          The order in which the interchanges are made is N to IHI+1,
 106 *>          then 1 to ILO-1.
 107 *> \endverbatim
 108 *>
 109 *> \param[out] INFO
 110 *> \verbatim
 111 *>          INFO is INTEGER
 112 *>          = 0:  successful exit.
 113 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 114 *> \endverbatim
 115 *
 116 *  Authors:
 117 *  ========
 118 *
 119 *> \author Univ. of Tennessee
 120 *> \author Univ. of California Berkeley
 121 *> \author Univ. of Colorado Denver
 122 *> \author NAG Ltd.
 123 *
 124 *> \date November 2015
 125 *
 126 *> \ingroup doubleGEcomputational
 127 *
 128 *> \par Further Details:
 129 *  =====================
 130 *>
 131 *> \verbatim
 132 *>
 133 *>  The permutations consist of row and column interchanges which put
 134 *>  the matrix in the form
 135 *>
 136 *>             ( T1   X   Y  )
 137 *>     P A P = (  0   B   Z  )
 138 *>             (  0   0   T2 )
 139 *>
 140 *>  where T1 and T2 are upper triangular matrices whose eigenvalues lie
 141 *>  along the diagonal.  The column indices ILO and IHI mark the starting
 142 *>  and ending columns of the submatrix B. Balancing consists of applying
 143 *>  a diagonal similarity transformation inv(D) * B * D to make the
 144 *>  1-norms of each row of B and its corresponding column nearly equal.
 145 *>  The output matrix is
 146 *>
 147 *>     ( T1     X*D          Y    )
 148 *>     (  0  inv(D)*B*D  inv(D)*Z ).
 149 *>     (  0      0           T2   )
 150 *>
 151 *>  Information about the permutations P and the diagonal matrix D is
 152 *>  returned in the vector SCALE.
 153 *>
 154 *>  This subroutine is based on the EISPACK routine BALANC.
 155 *>
 156 *>  Modified by Tzu-Yi Chen, Computer Science Division, University of
 157 *>    California at Berkeley, USA
 158 *> \endverbatim
 159 *>
 160 *  =====================================================================
 161       SUBROUTINE DGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
 162 *
 163 *  -- LAPACK computational routine (version 3.6.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 2015
 167 *
 168 *     .. Scalar Arguments ..
 169       CHARACTER          JOB
 170       INTEGER            IHI, ILO, INFO, LDA, N
 171 *     ..
 172 *     .. Array Arguments ..
 173       DOUBLE PRECISION   A( LDA, * ), SCALE( * )
 174 *     ..
 175 *
 176 *  =====================================================================
 177 *
 178 *     .. Parameters ..
 179       DOUBLE PRECISION   ZERO, ONE
 180       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 181       DOUBLE PRECISION   SCLFAC
 182       PARAMETER          ( SCLFAC = 2.0D+0 )
 183       DOUBLE PRECISION   FACTOR
 184       PARAMETER          ( FACTOR = 0.95D+0 )
 185 *     ..
 186 *     .. Local Scalars ..
 187       LOGICAL            NOCONV
 188       INTEGER            I, ICA, IEXC, IRA, J, K, L, M
 189       DOUBLE PRECISION   C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
 190      $                   SFMIN2
 191 *     ..
 192 *     .. External Functions ..
 193       LOGICAL            DISNAN, LSAME
 194       INTEGER            IDAMAX
 195       DOUBLE PRECISION   DLAMCH, DNRM2
 196       EXTERNAL           DISNAN, LSAME, IDAMAX, DLAMCH, DNRM2
 197 *     ..
 198 *     .. External Subroutines ..
 199       EXTERNAL           DSCAL, DSWAP, XERBLA
 200 *     ..
 201 *     .. Intrinsic Functions ..
 202       INTRINSIC          ABS, MAX, MIN
 203 *     ..
 204 *     Test the input parameters
 205 *
 206       INFO = 0
 207       IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
 208      $    .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
 209          INFO = -1
 210       ELSE IF( N.LT.0 ) THEN
 211          INFO = -2
 212       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 213          INFO = -4
 214       END IF
 215       IF( INFO.NE.0 ) THEN
 216          CALL XERBLA( 'DGEBAL', -INFO )
 217          RETURN
 218       END IF
 219 *
 220       K = 1
 221       L = N
 222 *
 223       IF( N.EQ.0 )
 224      $   GO TO 210
 225 *
 226       IF( LSAME( JOB, 'N' ) ) THEN
 227          DO 10 I = 1, N
 228             SCALE( I ) = ONE
 229    10    CONTINUE
 230          GO TO 210
 231       END IF
 232 *
 233       IF( LSAME( JOB, 'S' ) )
 234      $   GO TO 120
 235 *
 236 *     Permutation to isolate eigenvalues if possible
 237 *
 238       GO TO 50
 239 *
 240 *     Row and column exchange.
 241 *
 242    20 CONTINUE
 243       SCALE( M ) = J
 244       IF( J.EQ.M )
 245      $   GO TO 30
 246 *
 247       CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
 248       CALL DSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
 249 *
 250    30 CONTINUE
 251       GO TO ( 40, 80 )IEXC
 252 *
 253 *     Search for rows isolating an eigenvalue and push them down.
 254 *
 255    40 CONTINUE
 256       IF( L.EQ.1 )
 257      $   GO TO 210
 258       L = L - 1
 259 *
 260    50 CONTINUE
 261       DO 70 J = L, 1, -1
 262 *
 263          DO 60 I = 1, L
 264             IF( I.EQ.J )
 265      $         GO TO 60
 266             IF( A( J, I ).NE.ZERO )
 267      $         GO TO 70
 268    60    CONTINUE
 269 *
 270          M = L
 271          IEXC = 1
 272          GO TO 20
 273    70 CONTINUE
 274 *
 275       GO TO 90
 276 *
 277 *     Search for columns isolating an eigenvalue and push them left.
 278 *
 279    80 CONTINUE
 280       K = K + 1
 281 *
 282    90 CONTINUE
 283       DO 110 J = K, L
 284 *
 285          DO 100 I = K, L
 286             IF( I.EQ.J )
 287      $         GO TO 100
 288             IF( A( I, J ).NE.ZERO )
 289      $         GO TO 110
 290   100    CONTINUE
 291 *
 292          M = K
 293          IEXC = 2
 294          GO TO 20
 295   110 CONTINUE
 296 *
 297   120 CONTINUE
 298       DO 130 I = K, L
 299          SCALE( I ) = ONE
 300   130 CONTINUE
 301 *
 302       IF( LSAME( JOB, 'P' ) )
 303      $   GO TO 210
 304 *
 305 *     Balance the submatrix in rows K to L.
 306 *
 307 *     Iterative loop for norm reduction
 308 *
 309       SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' )
 310       SFMAX1 = ONE / SFMIN1
 311       SFMIN2 = SFMIN1*SCLFAC
 312       SFMAX2 = ONE / SFMIN2
 313 *
 314   140 CONTINUE
 315       NOCONV = .FALSE.
 316 *
 317       DO 200 I = K, L
 318 *
 319          C = DNRM2( L-K+1, A( K, I ), 1 )
 320          R = DNRM2( L-K+1, A( I, K ), LDA )
 321          ICA = IDAMAX( L, A( 1, I ), 1 )
 322          CA = ABS( A( ICA, I ) )
 323          IRA = IDAMAX( N-K+1, A( I, K ), LDA )
 324          RA = ABS( A( I, IRA+K-1 ) )
 325 *
 326 *        Guard against zero C or R due to underflow.
 327 *
 328          IF( C.EQ.ZERO .OR. R.EQ.ZERO )
 329      $      GO TO 200
 330          G = R / SCLFAC
 331          F = ONE
 332          S = C + R
 333   160    CONTINUE
 334          IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
 335      $       MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
 336             IF( DISNAN( C+F+CA+R+G+RA ) ) THEN
 337 *
 338 *           Exit if NaN to avoid infinite loop
 339 *
 340             INFO = -3
 341             CALL XERBLA( 'DGEBAL', -INFO )
 342             RETURN
 343          END IF
 344          F = F*SCLFAC
 345          C = C*SCLFAC
 346          CA = CA*SCLFAC
 347          R = R / SCLFAC
 348          G = G / SCLFAC
 349          RA = RA / SCLFAC
 350          GO TO 160
 351 *
 352   170    CONTINUE
 353          G = C / SCLFAC
 354   180    CONTINUE
 355          IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR.
 356      $       MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190
 357          F = F / SCLFAC
 358          C = C / SCLFAC
 359          G = G / SCLFAC
 360          CA = CA / SCLFAC
 361          R = R*SCLFAC
 362          RA = RA*SCLFAC
 363          GO TO 180
 364 *
 365 *        Now balance.
 366 *
 367   190    CONTINUE
 368          IF( ( C+R ).GE.FACTOR*S )
 369      $      GO TO 200
 370          IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
 371             IF( F*SCALE( I ).LE.SFMIN1 )
 372      $         GO TO 200
 373          END IF
 374          IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
 375             IF( SCALE( I ).GE.SFMAX1 / F )
 376      $         GO TO 200
 377          END IF
 378          G = ONE / F
 379          SCALE( I ) = SCALE( I )*F
 380          NOCONV = .TRUE.
 381 *
 382          CALL DSCAL( N-K+1, G, A( I, K ), LDA )
 383          CALL DSCAL( L, F, A( 1, I ), 1 )
 384 *
 385   200 CONTINUE
 386 *
 387       IF( NOCONV )
 388      $   GO TO 140
 389 *
 390   210 CONTINUE
 391       ILO = K
 392       IHI = L
 393 *
 394       RETURN
 395 *
 396 *     End of DGEBAL
 397 *
 398       END