[BACK]Return to dhgeqz.f CVS log [TXT] [DIR] Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/dhgeqz.f between versions 1.7 and 1.21

version 1.7, 2010/12/21 13:53:27 version 1.21, 2023/08/07 08:38:52
Line 1 Line 1
   *> \brief \b DHGEQZ
   *
   *  =========== DOCUMENTATION ===========
   *
   * Online html documentation available at
   *            http://www.netlib.org/lapack/explore-html/
   *
   *> \htmlonly
   *> Download DHGEQZ + dependencies
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dhgeqz.f">
   *> [TGZ]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dhgeqz.f">
   *> [ZIP]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dhgeqz.f">
   *> [TXT]</a>
   *> \endhtmlonly
   *
   *  Definition:
   *  ===========
   *
   *       SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
   *                          ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
   *                          LWORK, INFO )
   *
   *       .. Scalar Arguments ..
   *       CHARACTER          COMPQ, COMPZ, JOB
   *       INTEGER            IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
   *       ..
   *       .. Array Arguments ..
   *       DOUBLE PRECISION   ALPHAI( * ), ALPHAR( * ), BETA( * ),
   *      $                   H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
   *      $                   WORK( * ), Z( LDZ, * )
   *       ..
   *
   *
   *> \par Purpose:
   *  =============
   *>
   *> \verbatim
   *>
   *> DHGEQZ computes the eigenvalues of a real matrix pair (H,T),
   *> where H is an upper Hessenberg matrix and T is upper triangular,
   *> using the double-shift QZ method.
   *> Matrix pairs of this type are produced by the reduction to
   *> generalized upper Hessenberg form of a real matrix pair (A,B):
   *>
   *>    A = Q1*H*Z1**T,  B = Q1*T*Z1**T,
   *>
   *> as computed by DGGHRD.
   *>
   *> If JOB='S', then the Hessenberg-triangular pair (H,T) is
   *> also reduced to generalized Schur form,
   *>
   *>    H = Q*S*Z**T,  T = Q*P*Z**T,
   *>
   *> where Q and Z are orthogonal matrices, P is an upper triangular
   *> matrix, and S is a quasi-triangular matrix with 1-by-1 and 2-by-2
   *> diagonal blocks.
   *>
   *> The 1-by-1 blocks correspond to real eigenvalues of the matrix pair
   *> (H,T) and the 2-by-2 blocks correspond to complex conjugate pairs of
   *> eigenvalues.
   *>
   *> Additionally, the 2-by-2 upper triangular diagonal blocks of P
   *> corresponding to 2-by-2 blocks of S are reduced to positive diagonal
   *> form, i.e., if S(j+1,j) is non-zero, then P(j+1,j) = P(j,j+1) = 0,
   *> P(j,j) > 0, and P(j+1,j+1) > 0.
   *>
   *> Optionally, the orthogonal matrix Q from the generalized Schur
   *> factorization may be postmultiplied into an input matrix Q1, and the
   *> orthogonal matrix Z may be postmultiplied into an input matrix Z1.
   *> If Q1 and Z1 are the orthogonal matrices from DGGHRD that reduced
   *> the matrix pair (A,B) to generalized upper Hessenberg form, then the
   *> output matrices Q1*Q and Z1*Z are the orthogonal factors from the
   *> generalized Schur factorization of (A,B):
   *>
   *>    A = (Q1*Q)*S*(Z1*Z)**T,  B = (Q1*Q)*P*(Z1*Z)**T.
   *>
   *> To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently,
   *> of (A,B)) are computed as a pair of values (alpha,beta), where alpha is
   *> complex and beta real.
   *> If beta is nonzero, lambda = alpha / beta is an eigenvalue of the
   *> generalized nonsymmetric eigenvalue problem (GNEP)
   *>    A*x = lambda*B*x
   *> and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the
   *> alternate form of the GNEP
   *>    mu*A*y = B*y.
   *> Real eigenvalues can be read directly from the generalized Schur
   *> form:
   *>   alpha = S(i,i), beta = P(i,i).
   *>
   *> Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
   *>      Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
   *>      pp. 241--256.
   *> \endverbatim
   *
   *  Arguments:
   *  ==========
   *
   *> \param[in] JOB
   *> \verbatim
   *>          JOB is CHARACTER*1
   *>          = 'E': Compute eigenvalues only;
   *>          = 'S': Compute eigenvalues and the Schur form.
   *> \endverbatim
   *>
   *> \param[in] COMPQ
   *> \verbatim
   *>          COMPQ is CHARACTER*1
   *>          = 'N': Left Schur vectors (Q) are not computed;
   *>          = 'I': Q is initialized to the unit matrix and the matrix Q
   *>                 of left Schur vectors of (H,T) is returned;
   *>          = 'V': Q must contain an orthogonal matrix Q1 on entry and
   *>                 the product Q1*Q is returned.
   *> \endverbatim
   *>
   *> \param[in] COMPZ
   *> \verbatim
   *>          COMPZ is CHARACTER*1
   *>          = 'N': Right Schur vectors (Z) are not computed;
   *>          = 'I': Z is initialized to the unit matrix and the matrix Z
   *>                 of right Schur vectors of (H,T) is returned;
   *>          = 'V': Z must contain an orthogonal matrix Z1 on entry and
   *>                 the product Z1*Z is returned.
   *> \endverbatim
   *>
   *> \param[in] N
   *> \verbatim
   *>          N is INTEGER
   *>          The order of the matrices H, T, Q, and Z.  N >= 0.
   *> \endverbatim
   *>
   *> \param[in] ILO
   *> \verbatim
   *>          ILO is INTEGER
   *> \endverbatim
   *>
   *> \param[in] IHI
   *> \verbatim
   *>          IHI is INTEGER
   *>          ILO and IHI mark the rows and columns of H which are in
   *>          Hessenberg form.  It is assumed that A is already upper
   *>          triangular in rows and columns 1:ILO-1 and IHI+1:N.
   *>          If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0.
   *> \endverbatim
   *>
   *> \param[in,out] H
   *> \verbatim
   *>          H is DOUBLE PRECISION array, dimension (LDH, N)
   *>          On entry, the N-by-N upper Hessenberg matrix H.
   *>          On exit, if JOB = 'S', H contains the upper quasi-triangular
   *>          matrix S from the generalized Schur factorization.
   *>          If JOB = 'E', the diagonal blocks of H match those of S, but
   *>          the rest of H is unspecified.
   *> \endverbatim
   *>
   *> \param[in] LDH
   *> \verbatim
   *>          LDH is INTEGER
   *>          The leading dimension of the array H.  LDH >= max( 1, N ).
   *> \endverbatim
   *>
   *> \param[in,out] T
   *> \verbatim
   *>          T is DOUBLE PRECISION array, dimension (LDT, N)
   *>          On entry, the N-by-N upper triangular matrix T.
   *>          On exit, if JOB = 'S', T contains the upper triangular
   *>          matrix P from the generalized Schur factorization;
   *>          2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks of S
   *>          are reduced to positive diagonal form, i.e., if H(j+1,j) is
   *>          non-zero, then T(j+1,j) = T(j,j+1) = 0, T(j,j) > 0, and
   *>          T(j+1,j+1) > 0.
   *>          If JOB = 'E', the diagonal blocks of T match those of P, but
   *>          the rest of T is unspecified.
   *> \endverbatim
   *>
   *> \param[in] LDT
   *> \verbatim
   *>          LDT is INTEGER
   *>          The leading dimension of the array T.  LDT >= max( 1, N ).
   *> \endverbatim
   *>
   *> \param[out] ALPHAR
   *> \verbatim
   *>          ALPHAR is DOUBLE PRECISION array, dimension (N)
   *>          The real parts of each scalar alpha defining an eigenvalue
   *>          of GNEP.
   *> \endverbatim
   *>
   *> \param[out] ALPHAI
   *> \verbatim
   *>          ALPHAI is DOUBLE PRECISION array, dimension (N)
   *>          The imaginary parts of each scalar alpha defining an
   *>          eigenvalue of GNEP.
   *>          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
   *>          positive, then the j-th and (j+1)-st eigenvalues are a
   *>          complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).
   *> \endverbatim
   *>
   *> \param[out] BETA
   *> \verbatim
   *>          BETA is DOUBLE PRECISION array, dimension (N)
   *>          The scalars beta that define the eigenvalues of GNEP.
   *>          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
   *>          beta = BETA(j) represent the j-th eigenvalue of the matrix
   *>          pair (A,B), in one of the forms lambda = alpha/beta or
   *>          mu = beta/alpha.  Since either lambda or mu may overflow,
   *>          they should not, in general, be computed.
   *> \endverbatim
   *>
   *> \param[in,out] Q
   *> \verbatim
   *>          Q is DOUBLE PRECISION array, dimension (LDQ, N)
   *>          On entry, if COMPQ = 'V', the orthogonal matrix Q1 used in
   *>          the reduction of (A,B) to generalized Hessenberg form.
   *>          On exit, if COMPQ = 'I', the orthogonal matrix of left Schur
   *>          vectors of (H,T), and if COMPQ = 'V', the orthogonal matrix
   *>          of left Schur vectors of (A,B).
   *>          Not referenced if COMPQ = 'N'.
   *> \endverbatim
   *>
   *> \param[in] LDQ
   *> \verbatim
   *>          LDQ is INTEGER
   *>          The leading dimension of the array Q.  LDQ >= 1.
   *>          If COMPQ='V' or 'I', then LDQ >= N.
   *> \endverbatim
   *>
   *> \param[in,out] Z
   *> \verbatim
   *>          Z is DOUBLE PRECISION array, dimension (LDZ, N)
   *>          On entry, if COMPZ = 'V', the orthogonal matrix Z1 used in
   *>          the reduction of (A,B) to generalized Hessenberg form.
   *>          On exit, if COMPZ = 'I', the orthogonal matrix of
   *>          right Schur vectors of (H,T), and if COMPZ = 'V', the
   *>          orthogonal matrix of right Schur vectors of (A,B).
   *>          Not referenced if COMPZ = 'N'.
   *> \endverbatim
   *>
   *> \param[in] LDZ
   *> \verbatim
   *>          LDZ is INTEGER
   *>          The leading dimension of the array Z.  LDZ >= 1.
   *>          If COMPZ='V' or 'I', then LDZ >= N.
   *> \endverbatim
   *>
   *> \param[out] WORK
   *> \verbatim
   *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
   *>          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
   *> \endverbatim
   *>
   *> \param[in] LWORK
   *> \verbatim
   *>          LWORK is INTEGER
   *>          The dimension of the array WORK.  LWORK >= max(1,N).
   *>
   *>          If LWORK = -1, then a workspace query is assumed; the routine
   *>          only calculates the optimal size of the WORK array, returns
   *>          this value as the first entry of the WORK array, and no error
   *>          message related to LWORK is issued by XERBLA.
   *> \endverbatim
   *>
   *> \param[out] INFO
   *> \verbatim
   *>          INFO is INTEGER
   *>          = 0: successful exit
   *>          < 0: if INFO = -i, the i-th argument had an illegal value
   *>          = 1,...,N: the QZ iteration did not converge.  (H,T) is not
   *>                     in Schur form, but ALPHAR(i), ALPHAI(i), and
   *>                     BETA(i), i=INFO+1,...,N should be correct.
   *>          = N+1,...,2*N: the shift calculation failed.  (H,T) is not
   *>                     in Schur form, but ALPHAR(i), ALPHAI(i), and
   *>                     BETA(i), i=INFO-N+1,...,N should be correct.
   *> \endverbatim
   *
   *  Authors:
   *  ========
   *
   *> \author Univ. of Tennessee
   *> \author Univ. of California Berkeley
   *> \author Univ. of Colorado Denver
   *> \author NAG Ltd.
   *
   *> \ingroup doubleGEcomputational
   *
   *> \par Further Details:
   *  =====================
   *>
   *> \verbatim
   *>
   *>  Iteration counters:
   *>
   *>  JITER  -- counts iterations.
   *>  IITER  -- counts iterations run since ILAST was last
   *>            changed.  This is therefore reset only when a 1-by-1 or
   *>            2-by-2 block deflates off the bottom.
   *> \endverbatim
   *>
   *  =====================================================================
       SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,        SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
      $                   ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,       $                   ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
      $                   LWORK, INFO )       $                   LWORK, INFO )
 *  *
 *  -- LAPACK routine (version 3.2.1)                                  --  *  -- LAPACK computational routine --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 *  -- April 2009                                                      --  
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          COMPQ, COMPZ, JOB        CHARACTER          COMPQ, COMPZ, JOB
Line 17 Line 316
      $                   WORK( * ), Z( LDZ, * )       $                   WORK( * ), Z( LDZ, * )
 *     ..  *     ..
 *  *
 *  Purpose  
 *  =======  
 *  
 *  DHGEQZ computes the eigenvalues of a real matrix pair (H,T),  
 *  where H is an upper Hessenberg matrix and T is upper triangular,  
 *  using the double-shift QZ method.  
 *  Matrix pairs of this type are produced by the reduction to  
 *  generalized upper Hessenberg form of a real matrix pair (A,B):  
 *  
 *     A = Q1*H*Z1**T,  B = Q1*T*Z1**T,  
 *  
 *  as computed by DGGHRD.  
 *  
 *  If JOB='S', then the Hessenberg-triangular pair (H,T) is  
 *  also reduced to generalized Schur form,  
 *    
 *     H = Q*S*Z**T,  T = Q*P*Z**T,  
 *    
 *  where Q and Z are orthogonal matrices, P is an upper triangular  
 *  matrix, and S is a quasi-triangular matrix with 1-by-1 and 2-by-2  
 *  diagonal blocks.  
 *  
 *  The 1-by-1 blocks correspond to real eigenvalues of the matrix pair  
 *  (H,T) and the 2-by-2 blocks correspond to complex conjugate pairs of  
 *  eigenvalues.  
 *  
 *  Additionally, the 2-by-2 upper triangular diagonal blocks of P  
 *  corresponding to 2-by-2 blocks of S are reduced to positive diagonal  
 *  form, i.e., if S(j+1,j) is non-zero, then P(j+1,j) = P(j,j+1) = 0,  
 *  P(j,j) > 0, and P(j+1,j+1) > 0.  
 *  
 *  Optionally, the orthogonal matrix Q from the generalized Schur  
 *  factorization may be postmultiplied into an input matrix Q1, and the  
 *  orthogonal matrix Z may be postmultiplied into an input matrix Z1.  
 *  If Q1 and Z1 are the orthogonal matrices from DGGHRD that reduced  
 *  the matrix pair (A,B) to generalized upper Hessenberg form, then the  
 *  output matrices Q1*Q and Z1*Z are the orthogonal factors from the  
 *  generalized Schur factorization of (A,B):  
 *  
 *     A = (Q1*Q)*S*(Z1*Z)**T,  B = (Q1*Q)*P*(Z1*Z)**T.  
 *    
 *  To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently,  
 *  of (A,B)) are computed as a pair of values (alpha,beta), where alpha is  
 *  complex and beta real.  
 *  If beta is nonzero, lambda = alpha / beta is an eigenvalue of the  
 *  generalized nonsymmetric eigenvalue problem (GNEP)  
 *     A*x = lambda*B*x  
 *  and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the  
 *  alternate form of the GNEP  
 *     mu*A*y = B*y.  
 *  Real eigenvalues can be read directly from the generalized Schur  
 *  form:   
 *    alpha = S(i,i), beta = P(i,i).  
 *  
 *  Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix  
 *       Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),  
 *       pp. 241--256.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  JOB     (input) CHARACTER*1  
 *          = 'E': Compute eigenvalues only;  
 *          = 'S': Compute eigenvalues and the Schur form.   
 *  
 *  COMPQ   (input) CHARACTER*1  
 *          = 'N': Left Schur vectors (Q) are not computed;  
 *          = 'I': Q is initialized to the unit matrix and the matrix Q  
 *                 of left Schur vectors of (H,T) is returned;  
 *          = 'V': Q must contain an orthogonal matrix Q1 on entry and  
 *                 the product Q1*Q is returned.  
 *  
 *  COMPZ   (input) CHARACTER*1  
 *          = 'N': Right Schur vectors (Z) are not computed;  
 *          = 'I': Z is initialized to the unit matrix and the matrix Z  
 *                 of right Schur vectors of (H,T) is returned;  
 *          = 'V': Z must contain an orthogonal matrix Z1 on entry and  
 *                 the product Z1*Z is returned.  
 *  
 *  N       (input) INTEGER  
 *          The order of the matrices H, T, Q, and Z.  N >= 0.  
 *  
 *  ILO     (input) INTEGER  
 *  IHI     (input) INTEGER  
 *          ILO and IHI mark the rows and columns of H which are in  
 *          Hessenberg form.  It is assumed that A is already upper  
 *          triangular in rows and columns 1:ILO-1 and IHI+1:N.  
 *          If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0.  
 *  
 *  H       (input/output) DOUBLE PRECISION array, dimension (LDH, N)  
 *          On entry, the N-by-N upper Hessenberg matrix H.  
 *          On exit, if JOB = 'S', H contains the upper quasi-triangular  
 *          matrix S from the generalized Schur factorization;  
 *          2-by-2 diagonal blocks (corresponding to complex conjugate  
 *          pairs of eigenvalues) are returned in standard form, with  
 *          H(i,i) = H(i+1,i+1) and H(i+1,i)*H(i,i+1) < 0.  
 *          If JOB = 'E', the diagonal blocks of H match those of S, but  
 *          the rest of H is unspecified.  
 *  
 *  LDH     (input) INTEGER  
 *          The leading dimension of the array H.  LDH >= max( 1, N ).  
 *  
 *  T       (input/output) DOUBLE PRECISION array, dimension (LDT, N)  
 *          On entry, the N-by-N upper triangular matrix T.  
 *          On exit, if JOB = 'S', T contains the upper triangular  
 *          matrix P from the generalized Schur factorization;  
 *          2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks of S  
 *          are reduced to positive diagonal form, i.e., if H(j+1,j) is  
 *          non-zero, then T(j+1,j) = T(j,j+1) = 0, T(j,j) > 0, and  
 *          T(j+1,j+1) > 0.  
 *          If JOB = 'E', the diagonal blocks of T match those of P, but  
 *          the rest of T is unspecified.  
 *  
 *  LDT     (input) INTEGER  
 *          The leading dimension of the array T.  LDT >= max( 1, N ).  
 *  
 *  ALPHAR  (output) DOUBLE PRECISION array, dimension (N)  
 *          The real parts of each scalar alpha defining an eigenvalue  
 *          of GNEP.  
 *  
 *  ALPHAI  (output) DOUBLE PRECISION array, dimension (N)  
 *          The imaginary parts of each scalar alpha defining an  
 *          eigenvalue of GNEP.  
 *          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if  
 *          positive, then the j-th and (j+1)-st eigenvalues are a  
 *          complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).  
 *  
 *  BETA    (output) DOUBLE PRECISION array, dimension (N)  
 *          The scalars beta that define the eigenvalues of GNEP.  
 *          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and  
 *          beta = BETA(j) represent the j-th eigenvalue of the matrix  
 *          pair (A,B), in one of the forms lambda = alpha/beta or  
 *          mu = beta/alpha.  Since either lambda or mu may overflow,  
 *          they should not, in general, be computed.  
 *  
 *  Q       (input/output) DOUBLE PRECISION array, dimension (LDQ, N)  
 *          On entry, if COMPZ = 'V', the orthogonal matrix Q1 used in  
 *          the reduction of (A,B) to generalized Hessenberg form.  
 *          On exit, if COMPZ = 'I', the orthogonal matrix of left Schur  
 *          vectors of (H,T), and if COMPZ = 'V', the orthogonal matrix  
 *          of left Schur vectors of (A,B).  
 *          Not referenced if COMPZ = 'N'.  
 *  
 *  LDQ     (input) INTEGER  
 *          The leading dimension of the array Q.  LDQ >= 1.  
 *          If COMPQ='V' or 'I', then LDQ >= N.  
 *  
 *  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ, N)  
 *          On entry, if COMPZ = 'V', the orthogonal matrix Z1 used in  
 *          the reduction of (A,B) to generalized Hessenberg form.  
 *          On exit, if COMPZ = 'I', the orthogonal matrix of  
 *          right Schur vectors of (H,T), and if COMPZ = 'V', the  
 *          orthogonal matrix of right Schur vectors of (A,B).  
 *          Not referenced if COMPZ = 'N'.  
 *  
 *  LDZ     (input) INTEGER  
 *          The leading dimension of the array Z.  LDZ >= 1.  
 *          If COMPZ='V' or 'I', then LDZ >= N.  
 *  
 *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))  
 *          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.  
 *  
 *  LWORK   (input) INTEGER  
 *          The dimension of the array WORK.  LWORK >= max(1,N).  
 *  
 *          If LWORK = -1, then a workspace query is assumed; the routine  
 *          only calculates the optimal size of the WORK array, returns  
 *          this value as the first entry of the WORK array, and no error  
 *          message related to LWORK is issued by XERBLA.  
 *  
 *  INFO    (output) INTEGER  
 *          = 0: successful exit  
 *          < 0: if INFO = -i, the i-th argument had an illegal value  
 *          = 1,...,N: the QZ iteration did not converge.  (H,T) is not  
 *                     in Schur form, but ALPHAR(i), ALPHAI(i), and  
 *                     BETA(i), i=INFO+1,...,N should be correct.  
 *          = N+1,...,2*N: the shift calculation failed.  (H,T) is not  
 *                     in Schur form, but ALPHAR(i), ALPHAI(i), and  
 *                     BETA(i), i=INFO-N+1,...,N should be correct.  
 *  
 *  Further Details  
 *  ===============  
 *  
 *  Iteration counters:  
 *  
 *  JITER  -- counts iterations.  
 *  IITER  -- counts iterations run since ILAST was last  
 *            changed.  This is therefore reset only when a 1-by-1 or  
 *            2-by-2 block deflates off the bottom.  
 *  
 *  =====================================================================  *  =====================================================================
 *  *
 *     .. Parameters ..  *     .. Parameters ..
Line 419 Line 528
 *  *
             GO TO 80              GO TO 80
          ELSE           ELSE
             IF( ABS( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN              IF( ABS( H( ILAST, ILAST-1 ) ).LE.MAX( SAFMIN, ULP*( 
        $         ABS( H( ILAST, ILAST ) ) + ABS( H( ILAST-1, ILAST-1 ) ) 
        $         ) ) ) THEN
                H( ILAST, ILAST-1 ) = ZERO                 H( ILAST, ILAST-1 ) = ZERO
                GO TO 80                 GO TO 80
             END IF              END IF
Line 439 Line 550
             IF( J.EQ.ILO ) THEN              IF( J.EQ.ILO ) THEN
                ILAZRO = .TRUE.                 ILAZRO = .TRUE.
             ELSE              ELSE
                IF( ABS( H( J, J-1 ) ).LE.ATOL ) THEN                 IF( ABS( H( J, J-1 ) ).LE.MAX( SAFMIN, ULP*( 
        $         ABS( H( J, J ) ) + ABS( H( J-1, J-1 ) ) 
        $         ) ) ) THEN
                   H( J, J-1 ) = ZERO                    H( J, J-1 ) = ZERO
                   ILAZRO = .TRUE.                    ILAZRO = .TRUE.
                ELSE                 ELSE
Line 627 Line 740
 *           Exceptional shift.  Chosen for no particularly good reason.  *           Exceptional shift.  Chosen for no particularly good reason.
 *           (Single shift only.)  *           (Single shift only.)
 *  *
             IF( ( DBLE( MAXIT )*SAFMIN )*ABS( H( ILAST-1, ILAST ) ).LT.              IF( ( DBLE( MAXIT )*SAFMIN )*ABS( H( ILAST, ILAST-1 ) ).LT.
      $          ABS( T( ILAST-1, ILAST-1 ) ) ) THEN       $          ABS( T( ILAST-1, ILAST-1 ) ) ) THEN
                ESHIFT = ESHIFT + H( ILAST-1, ILAST ) /                 ESHIFT = H( ILAST, ILAST-1 ) /
      $                  T( ILAST-1, ILAST-1 )       $                  T( ILAST-1, ILAST-1 )
             ELSE              ELSE
                ESHIFT = ESHIFT + ONE / ( SAFMIN*DBLE( MAXIT ) )                 ESHIFT = ESHIFT + ONE / ( SAFMIN*DBLE( MAXIT ) )
Line 647 Line 760
      $                  T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1,       $                  T( ILAST-1, ILAST-1 ), LDT, SAFMIN*SAFETY, S1,
      $                  S2, WR, WR2, WI )       $                  S2, WR, WR2, WI )
 *  *
               IF ( ABS( (WR/S1)*T( ILAST, ILAST ) - H( ILAST, ILAST ) )
        $         .GT. ABS( (WR2/S2)*T( ILAST, ILAST )
        $         - H( ILAST, ILAST ) ) ) THEN
                  TEMP = WR
                  WR = WR2
                  WR2 = TEMP
                  TEMP = S1
                  S1 = S2
                  S2 = TEMP
               END IF
             TEMP = MAX( S1, SAFMIN*MAX( ONE, ABS( WR ), ABS( WI ) ) )              TEMP = MAX( S1, SAFMIN*MAX( ONE, ABS( WR ), ABS( WI ) ) )
             IF( WI.NE.ZERO )              IF( WI.NE.ZERO )
      $         GO TO 200       $         GO TO 200
Line 808 Line 931
                      Z( J, ILAST ) = -Z( J, ILAST )                       Z( J, ILAST ) = -Z( J, ILAST )
   220             CONTINUE    220             CONTINUE
                END IF                 END IF
                  B22 = -B22
             END IF              END IF
 *  *
 *           Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)  *           Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
Removed from v.1.7  
changed lines
  Added in v.1.21

CVSweb interface <joel.bertrand@systella.fr>