LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cgebal.f
Go to the documentation of this file.
00001 *> \brief \b CGEBAL
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download CGEBAL + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgebal.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgebal.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgebal.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       CHARACTER          JOB
00025 *       INTEGER            IHI, ILO, INFO, LDA, N
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       REAL               SCALE( * )
00029 *       COMPLEX            A( LDA, * )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> CGEBAL balances a general complex matrix A.  This involves, first,
00039 *> permuting A by a similarity transformation to isolate eigenvalues
00040 *> in the first 1 to ILO-1 and last IHI+1 to N elements on the
00041 *> diagonal; and second, applying a diagonal similarity transformation
00042 *> to rows and columns ILO to IHI to make the rows and columns as
00043 *> close in norm as possible.  Both steps are optional.
00044 *>
00045 *> Balancing may reduce the 1-norm of the matrix, and improve the
00046 *> accuracy of the computed eigenvalues and/or eigenvectors.
00047 *> \endverbatim
00048 *
00049 *  Arguments:
00050 *  ==========
00051 *
00052 *> \param[in] JOB
00053 *> \verbatim
00054 *>          JOB is CHARACTER*1
00055 *>          Specifies the operations to be performed on A:
00056 *>          = 'N':  none:  simply set ILO = 1, IHI = N, SCALE(I) = 1.0
00057 *>                  for i = 1,...,N;
00058 *>          = 'P':  permute only;
00059 *>          = 'S':  scale only;
00060 *>          = 'B':  both permute and scale.
00061 *> \endverbatim
00062 *>
00063 *> \param[in] N
00064 *> \verbatim
00065 *>          N is INTEGER
00066 *>          The order of the matrix A.  N >= 0.
00067 *> \endverbatim
00068 *>
00069 *> \param[in,out] A
00070 *> \verbatim
00071 *>          A is COMPLEX array, dimension (LDA,N)
00072 *>          On entry, the input matrix A.
00073 *>          On exit,  A is overwritten by the balanced matrix.
00074 *>          If JOB = 'N', A is not referenced.
00075 *>          See Further Details.
00076 *> \endverbatim
00077 *>
00078 *> \param[in] LDA
00079 *> \verbatim
00080 *>          LDA is INTEGER
00081 *>          The leading dimension of the array A.  LDA >= max(1,N).
00082 *> \endverbatim
00083 *>
00084 *> \param[out] ILO
00085 *> \verbatim
00086 *>          ILO is INTEGER
00087 *> \endverbatim
00088 *> \param[out] IHI
00089 *> \verbatim
00090 *>          IHI is INTEGER
00091 *>          ILO and IHI are set to integers such that on exit
00092 *>          A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
00093 *>          If JOB = 'N' or 'S', ILO = 1 and IHI = N.
00094 *> \endverbatim
00095 *>
00096 *> \param[out] SCALE
00097 *> \verbatim
00098 *>          SCALE is REAL array, dimension (N)
00099 *>          Details of the permutations and scaling factors applied to
00100 *>          A.  If P(j) is the index of the row and column interchanged
00101 *>          with row and column j and D(j) is the scaling factor
00102 *>          applied to row and column j, then
00103 *>          SCALE(j) = P(j)    for j = 1,...,ILO-1
00104 *>                   = D(j)    for j = ILO,...,IHI
00105 *>                   = P(j)    for j = IHI+1,...,N.
00106 *>          The order in which the interchanges are made is N to IHI+1,
00107 *>          then 1 to ILO-1.
00108 *> \endverbatim
00109 *>
00110 *> \param[out] INFO
00111 *> \verbatim
00112 *>          INFO is INTEGER
00113 *>          = 0:  successful exit.
00114 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00115 *> \endverbatim
00116 *
00117 *  Authors:
00118 *  ========
00119 *
00120 *> \author Univ. of Tennessee 
00121 *> \author Univ. of California Berkeley 
00122 *> \author Univ. of Colorado Denver 
00123 *> \author NAG Ltd. 
00124 *
00125 *> \date November 2011
00126 *
00127 *> \ingroup complexGEcomputational
00128 *
00129 *> \par Further Details:
00130 *  =====================
00131 *>
00132 *> \verbatim
00133 *>
00134 *>  The permutations consist of row and column interchanges which put
00135 *>  the matrix in the form
00136 *>
00137 *>             ( T1   X   Y  )
00138 *>     P A P = (  0   B   Z  )
00139 *>             (  0   0   T2 )
00140 *>
00141 *>  where T1 and T2 are upper triangular matrices whose eigenvalues lie
00142 *>  along the diagonal.  The column indices ILO and IHI mark the starting
00143 *>  and ending columns of the submatrix B. Balancing consists of applying
00144 *>  a diagonal similarity transformation inv(D) * B * D to make the
00145 *>  1-norms of each row of B and its corresponding column nearly equal.
00146 *>  The output matrix is
00147 *>
00148 *>     ( T1     X*D          Y    )
00149 *>     (  0  inv(D)*B*D  inv(D)*Z ).
00150 *>     (  0      0           T2   )
00151 *>
00152 *>  Information about the permutations P and the diagonal matrix D is
00153 *>  returned in the vector SCALE.
00154 *>
00155 *>  This subroutine is based on the EISPACK routine CBAL.
00156 *>
00157 *>  Modified by Tzu-Yi Chen, Computer Science Division, University of
00158 *>    California at Berkeley, USA
00159 *> \endverbatim
00160 *>
00161 *  =====================================================================
00162       SUBROUTINE CGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
00163 *
00164 *  -- LAPACK computational routine (version 3.4.0) --
00165 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00166 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00167 *     November 2011
00168 *
00169 *     .. Scalar Arguments ..
00170       CHARACTER          JOB
00171       INTEGER            IHI, ILO, INFO, LDA, N
00172 *     ..
00173 *     .. Array Arguments ..
00174       REAL               SCALE( * )
00175       COMPLEX            A( LDA, * )
00176 *     ..
00177 *
00178 *  =====================================================================
00179 *
00180 *     .. Parameters ..
00181       REAL               ZERO, ONE
00182       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00183       REAL               SCLFAC
00184       PARAMETER          ( SCLFAC = 2.0E+0 )
00185       REAL               FACTOR
00186       PARAMETER          ( FACTOR = 0.95E+0 )
00187 *     ..
00188 *     .. Local Scalars ..
00189       LOGICAL            NOCONV
00190       INTEGER            I, ICA, IEXC, IRA, J, K, L, M
00191       REAL               C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
00192      $                   SFMIN2
00193       COMPLEX            CDUM
00194 *     ..
00195 *     .. External Functions ..
00196       LOGICAL            SISNAN, LSAME
00197       INTEGER            ICAMAX
00198       REAL               SLAMCH
00199       EXTERNAL           SISNAN, LSAME, ICAMAX, SLAMCH
00200 *     ..
00201 *     .. External Subroutines ..
00202       EXTERNAL           CSSCAL, CSWAP, XERBLA
00203 *     ..
00204 *     .. Intrinsic Functions ..
00205       INTRINSIC          ABS, AIMAG, MAX, MIN, REAL
00206 *     ..
00207 *     .. Statement Functions ..
00208       REAL               CABS1
00209 *     ..
00210 *     .. Statement Function definitions ..
00211       CABS1( CDUM ) = ABS( REAL( CDUM ) ) + ABS( AIMAG( CDUM ) )
00212 *     ..
00213 *     .. Executable Statements ..
00214 *
00215 *     Test the input parameters
00216 *
00217       INFO = 0
00218       IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
00219      $    .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
00220          INFO = -1
00221       ELSE IF( N.LT.0 ) THEN
00222          INFO = -2
00223       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00224          INFO = -4
00225       END IF
00226       IF( INFO.NE.0 ) THEN
00227          CALL XERBLA( 'CGEBAL', -INFO )
00228          RETURN
00229       END IF
00230 *
00231       K = 1
00232       L = N
00233 *
00234       IF( N.EQ.0 )
00235      $   GO TO 210
00236 *
00237       IF( LSAME( JOB, 'N' ) ) THEN
00238          DO 10 I = 1, N
00239             SCALE( I ) = ONE
00240    10    CONTINUE
00241          GO TO 210
00242       END IF
00243 *
00244       IF( LSAME( JOB, 'S' ) )
00245      $   GO TO 120
00246 *
00247 *     Permutation to isolate eigenvalues if possible
00248 *
00249       GO TO 50
00250 *
00251 *     Row and column exchange.
00252 *
00253    20 CONTINUE
00254       SCALE( M ) = J
00255       IF( J.EQ.M )
00256      $   GO TO 30
00257 *
00258       CALL CSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
00259       CALL CSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
00260 *
00261    30 CONTINUE
00262       GO TO ( 40, 80 )IEXC
00263 *
00264 *     Search for rows isolating an eigenvalue and push them down.
00265 *
00266    40 CONTINUE
00267       IF( L.EQ.1 )
00268      $   GO TO 210
00269       L = L - 1
00270 *
00271    50 CONTINUE
00272       DO 70 J = L, 1, -1
00273 *
00274          DO 60 I = 1, L
00275             IF( I.EQ.J )
00276      $         GO TO 60
00277             IF( REAL( A( J, I ) ).NE.ZERO .OR. AIMAG( A( J, I ) ).NE.
00278      $          ZERO )GO TO 70
00279    60    CONTINUE
00280 *
00281          M = L
00282          IEXC = 1
00283          GO TO 20
00284    70 CONTINUE
00285 *
00286       GO TO 90
00287 *
00288 *     Search for columns isolating an eigenvalue and push them left.
00289 *
00290    80 CONTINUE
00291       K = K + 1
00292 *
00293    90 CONTINUE
00294       DO 110 J = K, L
00295 *
00296          DO 100 I = K, L
00297             IF( I.EQ.J )
00298      $         GO TO 100
00299             IF( REAL( A( I, J ) ).NE.ZERO .OR. AIMAG( A( I, J ) ).NE.
00300      $          ZERO )GO TO 110
00301   100    CONTINUE
00302 *
00303          M = K
00304          IEXC = 2
00305          GO TO 20
00306   110 CONTINUE
00307 *
00308   120 CONTINUE
00309       DO 130 I = K, L
00310          SCALE( I ) = ONE
00311   130 CONTINUE
00312 *
00313       IF( LSAME( JOB, 'P' ) )
00314      $   GO TO 210
00315 *
00316 *     Balance the submatrix in rows K to L.
00317 *
00318 *     Iterative loop for norm reduction
00319 *
00320       SFMIN1 = SLAMCH( 'S' ) / SLAMCH( 'P' )
00321       SFMAX1 = ONE / SFMIN1
00322       SFMIN2 = SFMIN1*SCLFAC
00323       SFMAX2 = ONE / SFMIN2
00324   140 CONTINUE
00325       NOCONV = .FALSE.
00326 *
00327       DO 200 I = K, L
00328          C = ZERO
00329          R = ZERO
00330 *
00331          DO 150 J = K, L
00332             IF( J.EQ.I )
00333      $         GO TO 150
00334             C = C + CABS1( A( J, I ) )
00335             R = R + CABS1( A( I, J ) )
00336   150    CONTINUE
00337          ICA = ICAMAX( L, A( 1, I ), 1 )
00338          CA = ABS( A( ICA, I ) )
00339          IRA = ICAMAX( N-K+1, A( I, K ), LDA )
00340          RA = ABS( A( I, IRA+K-1 ) )
00341 *
00342 *        Guard against zero C or R due to underflow.
00343 *
00344          IF( C.EQ.ZERO .OR. R.EQ.ZERO )
00345      $      GO TO 200
00346          G = R / SCLFAC
00347          F = ONE
00348          S = C + R
00349   160    CONTINUE
00350          IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
00351      $       MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
00352             IF( SISNAN( C+F+CA+R+G+RA ) ) THEN
00353 *
00354 *           Exit if NaN to avoid infinite loop
00355 *
00356             INFO = -3
00357             CALL XERBLA( 'CGEBAL', -INFO )
00358             RETURN
00359          END IF
00360          F = F*SCLFAC
00361          C = C*SCLFAC
00362          CA = CA*SCLFAC
00363          R = R / SCLFAC
00364          G = G / SCLFAC
00365          RA = RA / SCLFAC
00366          GO TO 160
00367 *
00368   170    CONTINUE
00369          G = C / SCLFAC
00370   180    CONTINUE
00371          IF( G.LT.R .OR. MAX( R, RA ).GE.SFMAX2 .OR.
00372      $       MIN( F, C, G, CA ).LE.SFMIN2 )GO TO 190
00373          F = F / SCLFAC
00374          C = C / SCLFAC
00375          G = G / SCLFAC
00376          CA = CA / SCLFAC
00377          R = R*SCLFAC
00378          RA = RA*SCLFAC
00379          GO TO 180
00380 *
00381 *        Now balance.
00382 *
00383   190    CONTINUE
00384          IF( ( C+R ).GE.FACTOR*S )
00385      $      GO TO 200
00386          IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
00387             IF( F*SCALE( I ).LE.SFMIN1 )
00388      $         GO TO 200
00389          END IF
00390          IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
00391             IF( SCALE( I ).GE.SFMAX1 / F )
00392      $         GO TO 200
00393          END IF
00394          G = ONE / F
00395          SCALE( I ) = SCALE( I )*F
00396          NOCONV = .TRUE.
00397 *
00398          CALL CSSCAL( N-K+1, G, A( I, K ), LDA )
00399          CALL CSSCAL( L, F, A( 1, I ), 1 )
00400 *
00401   200 CONTINUE
00402 *
00403       IF( NOCONV )
00404      $   GO TO 140
00405 *
00406   210 CONTINUE
00407       ILO = K
00408       IHI = L
00409 *
00410       RETURN
00411 *
00412 *     End of CGEBAL
00413 *
00414       END
 All Files Functions