LAPACK: SRC/dgebal.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
00002 *
00003 *  -- LAPACK routine (version 3.2.2) --
00004 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00005 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00006 *     June 2010
00007 *
00008 *     .. Scalar Arguments ..
00009       CHARACTER          JOB
00010       INTEGER            IHI, ILO, INFO, LDA, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   A( LDA, * ), SCALE( * )
00014 *     ..
00015 *
00016 *  Purpose
00017 *  =======
00018 *
00019 *  DGEBAL balances a general real matrix A.  This involves, first,
00020 *  permuting A by a similarity transformation to isolate eigenvalues
00021 *  in the first 1 to ILO-1 and last IHI+1 to N elements on the
00022 *  diagonal; and second, applying a diagonal similarity transformation
00023 *  to rows and columns ILO to IHI to make the rows and columns as
00024 *  close in norm as possible.  Both steps are optional.
00025 *
00026 *  Balancing may reduce the 1-norm of the matrix, and improve the
00027 *  accuracy of the computed eigenvalues and/or eigenvectors.
00028 *
00029 *  Arguments
00030 *  =========
00031 *
00032 *  JOB     (input) CHARACTER*1
00033 *          Specifies the operations to be performed on A:
00034 *          = 'N':  none:  simply set ILO = 1, IHI = N, SCALE(I) = 1.0
00035 *                  for i = 1,...,N;
00036 *          = 'P':  permute only;
00037 *          = 'S':  scale only;
00038 *          = 'B':  both permute and scale.
00039 *
00040 *  N       (input) INTEGER
00041 *          The order of the matrix A.  N >= 0.
00042 *
00043 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
00044 *          On entry, the input matrix A.
00045 *          On exit,  A is overwritten by the balanced matrix.
00046 *          If JOB = 'N', A is not referenced.
00047 *          See Further Details.
00048 *
00049 *  LDA     (input) INTEGER
00050 *          The leading dimension of the array A.  LDA >= max(1,N).
00051 *
00052 *  ILO     (output) INTEGER
00053 *  IHI     (output) INTEGER
00054 *          ILO and IHI are set to integers such that on exit
00055 *          A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
00056 *          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
00057 *
00058 *  SCALE   (output) DOUBLE PRECISION array, dimension (N)
00059 *          Details of the permutations and scaling factors applied to
00060 *          A.  If P(j) is the index of the row and column interchanged
00061 *          with row and column j and D(j) is the scaling factor
00062 *          applied to row and column j, then
00063 *          SCALE(j) = P(j)    for j = 1,...,ILO-1
00064 *                   = D(j)    for j = ILO,...,IHI
00065 *                   = P(j)    for j = IHI+1,...,N.
00066 *          The order in which the interchanges are made is N to IHI+1,
00067 *          then 1 to ILO-1.
00068 *
00069 *  INFO    (output) INTEGER
00070 *          = 0:  successful exit.
00071 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00072 *
00073 *  Further Details
00074 *  ===============
00075 *
00076 *  The permutations consist of row and column interchanges which put
00077 *  the matrix in the form
00078 *
00079 *             ( T1   X   Y  )
00080 *     P A P = (  0   B   Z  )
00081 *             (  0   0   T2 )
00082 *
00083 *  where T1 and T2 are upper triangular matrices whose eigenvalues lie
00084 *  along the diagonal.  The column indices ILO and IHI mark the starting
00085 *  and ending columns of the submatrix B. Balancing consists of applying
00086 *  a diagonal similarity transformation inv(D) * B * D to make the
00087 *  1-norms of each row of B and its corresponding column nearly equal.
00088 *  The output matrix is
00089 *
00090 *     ( T1     X*D          Y    )
00091 *     (  0  inv(D)*B*D  inv(D)*Z ).
00092 *     (  0      0           T2   )
00093 *
00094 *  Information about the permutations P and the diagonal matrix D is
00095 *  returned in the vector SCALE.
00096 *
00097 *  This subroutine is based on the EISPACK routine BALANC.
00098 *
00099 *  Modified by Tzu-Yi Chen, Computer Science Division, University of
00100 *    California at Berkeley, USA
00101 *
00102 *  =====================================================================
00103 *
00104 *     .. Parameters ..
00105       DOUBLE PRECISION   ZERO, ONE
00106       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
00107       DOUBLE PRECISION   SCLFAC
00108       PARAMETER          ( SCLFAC = 2.0D+0 )
00109       DOUBLE PRECISION   FACTOR
00110       PARAMETER          ( FACTOR = 0.95D+0 )
00111 *     ..
00112 *     .. Local Scalars ..
00113       LOGICAL            NOCONV
00114       INTEGER            I, ICA, IEXC, IRA, J, K, L, M
00115       DOUBLE PRECISION   C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
00116      $                   SFMIN2
00117 *     ..
00118 *     .. External Functions ..
00119       LOGICAL            DISNAN, LSAME
00120       INTEGER            IDAMAX
00121       DOUBLE PRECISION   DLAMCH
00122       EXTERNAL           DISNAN, LSAME, IDAMAX, DLAMCH
00123 *     ..
00124 *     .. External Subroutines ..
00125       EXTERNAL           DSCAL, DSWAP, XERBLA
00126 *     ..
00127 *     .. Intrinsic Functions ..
00128       INTRINSIC          ABS, MAX, MIN
00129 *     ..
00130 *     .. Executable Statements ..
00131 *
00132 *     Test the input parameters
00133 *
00134       INFO = 0
00135       IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
00136      $    .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
00137          INFO = -1
00138       ELSE IF( N.LT.0 ) THEN
00139          INFO = -2
00140       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00141          INFO = -4
00142       END IF
00143       IF( INFO.NE.0 ) THEN
00144          CALL XERBLA( 'DGEBAL', -INFO )
00145          RETURN
00146       END IF
00147 *
00148       K = 1
00149       L = N
00150 *
00151       IF( N.EQ.0 )
00152      $   GO TO 210
00153 *
00154       IF( LSAME( JOB, 'N' ) ) THEN
00155          DO 10 I = 1, N
00156             SCALE( I ) = ONE
00157    10    CONTINUE
00158          GO TO 210
00159       END IF
00160 *
00161       IF( LSAME( JOB, 'S' ) )
00162      $   GO TO 120
00163 *
00164 *     Permutation to isolate eigenvalues if possible
00165 *
00166       GO TO 50
00167 *
00168 *     Row and column exchange.
00169 *
00170    20 CONTINUE
00171       SCALE( M ) = J
00172       IF( J.EQ.M )
00173      $   GO TO 30
00174 *
00175       CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
00176       CALL DSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
00177 *
00178    30 CONTINUE
00179       GO TO ( 40, 80 )IEXC
00180 *
00181 *     Search for rows isolating an eigenvalue and push them down.
00182 *
00183    40 CONTINUE
00184       IF( L.EQ.1 )
00185      $   GO TO 210
00186       L = L - 1
00187 *
00188    50 CONTINUE
00189       DO 70 J = L, 1, -1
00190 *
00191          DO 60 I = 1, L
00192             IF( I.EQ.J )
00193      $         GO TO 60
00194             IF( A( J, I ).NE.ZERO )
00195      $         GO TO 70
00196    60    CONTINUE
00197 *
00198          M = L
00199          IEXC = 1
00200          GO TO 20
00201    70 CONTINUE
00202 *
00203       GO TO 90
00204 *
00205 *     Search for columns isolating an eigenvalue and push them left.
00206 *
00207    80 CONTINUE
00208       K = K + 1
00209 *
00210    90 CONTINUE
00211       DO 110 J = K, L
00212 *
00213          DO 100 I = K, L
00214             IF( I.EQ.J )
00215      $         GO TO 100
00216             IF( A( I, J ).NE.ZERO )
00217      $         GO TO 110
00218   100    CONTINUE
00219 *
00220          M = K
00221          IEXC = 2
00222          GO TO 20
00223   110 CONTINUE
00224 *
00225   120 CONTINUE
00226       DO 130 I = K, L
00227          SCALE( I ) = ONE
00228   130 CONTINUE
00229 *
00230       IF( LSAME( JOB, 'P' ) )
00231      $   GO TO 210
00232 *
00233 *     Balance the submatrix in rows K to L.
00234 *
00235 *     Iterative loop for norm reduction
00236 *
00237       SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' )
00238       SFMAX1 = ONE / SFMIN1
00239       SFMIN2 = SFMIN1*SCLFAC
00240       SFMAX2 = ONE / SFMIN2
00241   140 CONTINUE
00242       NOCONV = .FALSE.
00243 *
00244       DO 200 I = K, L
00245          C = ZERO
00246          R = ZERO
00247 *
00248          DO 150 J = K, L
00249             IF( J.EQ.I )
00250      $         GO TO 150
00251             C = C + ABS( A( J, I ) )
00252             R = R + ABS( A( I, J ) )
00253   150    CONTINUE
00254          ICA = IDAMAX( L, A( 1, I ), 1 )
00255          CA = ABS( A( ICA, I ) )
00256          IRA = IDAMAX( N-K+1, A( I, K ), LDA )
00257          RA = ABS( A( I, IRA+K-1 ) )
00258 *
00259 *        Guard against zero C or R due to underflow.
00260 *
00261          IF( C.EQ.ZERO .OR. R.EQ.ZERO )
00262      $      GO TO 200
00263          G = R / SCLFAC
00264          F = ONE
00265          S = C + R
00266   160    CONTINUE
00267          IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
00268      $       MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
00269             IF( DISNAN( C+F+CA+R+G+RA ) ) THEN
00270 *
00271 *           Exit if NaN to avoid infinite loop
00272 *
00273             INFO = -3
00274             CALL XERBLA( 'DGEBAL', -INFO )
00275             RETURN
00276          END IF
00277          F = F*SCLFAC
00278          C = C*SCLFAC
00279          CA = CA*SCLFAC
00280          R = R / SCLFAC
00281          G = G / SCLFAC
00282          RA = RA / SCLFAC
00283          GO TO 160
00284 *
00285   170    CONTINUE
00286          G = C / SCLFAC
00287   180    CONTINUE
00288          IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR.
00289      $       MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190
00290          F = F / SCLFAC
00291          C = C / SCLFAC
00292          G = G / SCLFAC
00293          CA = CA / SCLFAC
00294          R = R*SCLFAC
00295          RA = RA*SCLFAC
00296          GO TO 180
00297 *
00298 *        Now balance.
00299 *
00300   190    CONTINUE
00301          IF( ( C+R ).GE.FACTOR*S )
00302      $      GO TO 200
00303          IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
00304             IF( F*SCALE( I ).LE.SFMIN1 )
00305      $         GO TO 200
00306          END IF
00307          IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
00308             IF( SCALE( I ).GE.SFMAX1 / F )
00309      $         GO TO 200
00310          END IF
00311          G = ONE / F
00312          SCALE( I ) = SCALE( I )*F
00313          NOCONV = .TRUE.
00314 *
00315          CALL DSCAL( N-K+1, G, A( I, K ), LDA )
00316          CALL DSCAL( L, F, A( 1, I ), 1 )
00317 *
00318   200 CONTINUE
00319 *
00320       IF( NOCONV )
00321      $   GO TO 140
00322 *
00323   210 CONTINUE
00324       ILO = K
00325       IHI = L
00326 *
00327       RETURN
00328 *
00329 *     End of DGEBAL
00330 *
00331       END