LAPACK: SRC/zgebal.f Source File

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