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

Diff for /rpl/lapack/lapack/zhgeqz.f between versions 1.8 and 1.20

version 1.8, 2011/11/21 20:43:12 version 1.20, 2023/08/07 08:39:25
Line 2 Line 2
 *  *
 *  =========== DOCUMENTATION ===========  *  =========== DOCUMENTATION ===========
 *  *
 * Online html documentation available at   * Online html documentation available at
 *            http://www.netlib.org/lapack/explore-html/   *            http://www.netlib.org/lapack/explore-html/
 *  *
 *> \htmlonly  *> \htmlonly
 *> Download ZHGEQZ + dependencies   *> Download ZHGEQZ + dependencies
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhgeqz.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhgeqz.f">
 *> [TGZ]</a>   *> [TGZ]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhgeqz.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhgeqz.f">
 *> [ZIP]</a>   *> [ZIP]</a>
 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhgeqz.f">   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhgeqz.f">
 *> [TXT]</a>  *> [TXT]</a>
 *> \endhtmlonly   *> \endhtmlonly
 *  *
 *  Definition:  *  Definition:
 *  ===========  *  ===========
Line 21 Line 21
 *       SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,  *       SUBROUTINE ZHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
 *                          ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK,  *                          ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK,
 *                          RWORK, INFO )  *                          RWORK, INFO )
 *   *
 *       .. Scalar Arguments ..  *       .. Scalar Arguments ..
 *       CHARACTER          COMPQ, COMPZ, JOB  *       CHARACTER          COMPQ, COMPZ, JOB
 *       INTEGER            IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N  *       INTEGER            IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
Line 32 Line 32
 *      $                   Q( LDQ, * ), T( LDT, * ), WORK( * ),  *      $                   Q( LDQ, * ), T( LDT, * ), WORK( * ),
 *      $                   Z( LDZ, * )  *      $                   Z( LDZ, * )
 *       ..  *       ..
 *    *
 *  *
 *> \par Purpose:  *> \par Purpose:
 *  =============  *  =============
Line 44 Line 44
 *> using the single-shift QZ method.  *> using the single-shift QZ method.
 *> Matrix pairs of this type are produced by the reduction to  *> Matrix pairs of this type are produced by the reduction to
 *> generalized upper Hessenberg form of a complex matrix pair (A,B):  *> generalized upper Hessenberg form of a complex matrix pair (A,B):
 *>   *>
 *>    A = Q1*H*Z1**H,  B = Q1*T*Z1**H,  *>    A = Q1*H*Z1**H,  B = Q1*T*Z1**H,
 *>   *>
 *> as computed by ZGGHRD.  *> as computed by ZGGHRD.
 *>   *>
 *> If JOB='S', then the Hessenberg-triangular pair (H,T) is  *> If JOB='S', then the Hessenberg-triangular pair (H,T) is
 *> also reduced to generalized Schur form,  *> also reduced to generalized Schur form,
 *>   *>
 *>    H = Q*S*Z**H,  T = Q*P*Z**H,  *>    H = Q*S*Z**H,  T = Q*P*Z**H,
 *>   *>
 *> where Q and Z are unitary matrices and S and P are upper triangular.  *> where Q and Z are unitary matrices and S and P are upper triangular.
 *>   *>
 *> Optionally, the unitary matrix Q from the generalized Schur  *> Optionally, the unitary matrix Q from the generalized Schur
 *> factorization may be postmultiplied into an input matrix Q1, and the  *> factorization may be postmultiplied into an input matrix Q1, and the
 *> unitary matrix Z may be postmultiplied into an input matrix Z1.  *> unitary matrix Z may be postmultiplied into an input matrix Z1.
Line 63 Line 63
 *> the matrix pair (A,B) to generalized Hessenberg form, then the output  *> the matrix pair (A,B) to generalized Hessenberg form, then the output
 *> matrices Q1*Q and Z1*Z are the unitary factors from the generalized  *> matrices Q1*Q and Z1*Z are the unitary factors from the generalized
 *> Schur factorization of (A,B):  *> Schur factorization of (A,B):
 *>   *>
 *>    A = (Q1*Q)*S*(Z1*Z)**H,  B = (Q1*Q)*P*(Z1*Z)**H.  *>    A = (Q1*Q)*S*(Z1*Z)**H,  B = (Q1*Q)*P*(Z1*Z)**H.
 *>   *>
 *> To avoid overflow, eigenvalues of the matrix pair (H,T)  *> To avoid overflow, eigenvalues of the matrix pair (H,T)
 *> (equivalently, of (A,B)) are computed as a pair of complex values  *> (equivalently, of (A,B)) are computed as a pair of complex values
 *> (alpha,beta).  If beta is nonzero, lambda = alpha / beta is an  *> (alpha,beta).  If beta is nonzero, lambda = alpha / beta is an
Line 190 Line 190
 *> \param[in,out] Q  *> \param[in,out] Q
 *> \verbatim  *> \verbatim
 *>          Q is COMPLEX*16 array, dimension (LDQ, N)  *>          Q is COMPLEX*16 array, dimension (LDQ, N)
 *>          On entry, if COMPZ = 'V', the unitary matrix Q1 used in the  *>          On entry, if COMPQ = 'V', the unitary matrix Q1 used in the
 *>          reduction of (A,B) to generalized Hessenberg form.  *>          reduction of (A,B) to generalized Hessenberg form.
 *>          On exit, if COMPZ = 'I', the unitary matrix of left Schur  *>          On exit, if COMPQ = 'I', the unitary matrix of left Schur
 *>          vectors of (H,T), and if COMPZ = 'V', the unitary matrix of  *>          vectors of (H,T), and if COMPQ = 'V', the unitary matrix of
 *>          left Schur vectors of (A,B).  *>          left Schur vectors of (A,B).
 *>          Not referenced if COMPZ = 'N'.  *>          Not referenced if COMPQ = 'N'.
 *> \endverbatim  *> \endverbatim
 *>  *>
 *> \param[in] LDQ  *> \param[in] LDQ
Line 261 Line 261
 *  Authors:  *  Authors:
 *  ========  *  ========
 *  *
 *> \author Univ. of Tennessee   *> \author Univ. of Tennessee
 *> \author Univ. of California Berkeley   *> \author Univ. of California Berkeley
 *> \author Univ. of Colorado Denver   *> \author Univ. of Colorado Denver
 *> \author NAG Ltd.   *> \author NAG Ltd.
 *  
 *> \date November 2011  
 *  *
 *> \ingroup complex16GEcomputational  *> \ingroup complex16GEcomputational
 *  *
Line 284 Line 282
      $                   ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK,       $                   ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK,
      $                   RWORK, INFO )       $                   RWORK, INFO )
 *  *
 *  -- LAPACK computational routine (version 3.4.0) --  *  -- 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..--
 *     November 2011  
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          COMPQ, COMPZ, JOB        CHARACTER          COMPQ, COMPZ, JOB
Line 319 Line 316
       DOUBLE PRECISION   ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,        DOUBLE PRECISION   ABSB, ANORM, ASCALE, ATOL, BNORM, BSCALE, BTOL,
      $                   C, SAFMIN, TEMP, TEMP2, TEMPR, ULP       $                   C, SAFMIN, TEMP, TEMP2, TEMPR, ULP
       COMPLEX*16         ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,        COMPLEX*16         ABI22, AD11, AD12, AD21, AD22, CTEMP, CTEMP2,
      $                   CTEMP3, ESHIFT, RTDISC, S, SHIFT, SIGNBC, T1,       $                   CTEMP3, ESHIFT, S, SHIFT, SIGNBC,
      $                   U12, X       $                   U12, X, ABI12, Y
 *     ..  *     ..
 *     .. External Functions ..  *     .. External Functions ..
         COMPLEX*16         ZLADIV
       LOGICAL            LSAME        LOGICAL            LSAME
       DOUBLE PRECISION   DLAMCH, ZLANHS        DOUBLE PRECISION   DLAMCH, ZLANHS
       EXTERNAL           LSAME, DLAMCH, ZLANHS        EXTERNAL           ZLADIV, LSAME, DLAMCH, ZLANHS
 *     ..  *     ..
 *     .. External Subroutines ..  *     .. External Subroutines ..
       EXTERNAL           XERBLA, ZLARTG, ZLASET, ZROT, ZSCAL        EXTERNAL           XERBLA, ZLARTG, ZLASET, ZROT, ZSCAL
Line 351 Line 349
          ILSCHR = .TRUE.           ILSCHR = .TRUE.
          ISCHUR = 2           ISCHUR = 2
       ELSE        ELSE
            ILSCHR = .TRUE.
          ISCHUR = 0           ISCHUR = 0
       END IF        END IF
 *  *
Line 364 Line 363
          ILQ = .TRUE.           ILQ = .TRUE.
          ICOMPQ = 3           ICOMPQ = 3
       ELSE        ELSE
            ILQ = .TRUE.
          ICOMPQ = 0           ICOMPQ = 0
       END IF        END IF
 *  *
Line 377 Line 377
          ILZ = .TRUE.           ILZ = .TRUE.
          ICOMPZ = 3           ICOMPZ = 3
       ELSE        ELSE
            ILZ = .TRUE.
          ICOMPZ = 0           ICOMPZ = 0
       END IF        END IF
 *  *
Line 454 Line 455
                CALL ZSCAL( J-1, SIGNBC, T( 1, J ), 1 )                 CALL ZSCAL( J-1, SIGNBC, T( 1, J ), 1 )
                CALL ZSCAL( J, SIGNBC, H( 1, J ), 1 )                 CALL ZSCAL( J, SIGNBC, H( 1, J ), 1 )
             ELSE              ELSE
                H( J, J ) = H( J, J )*SIGNBC                 CALL ZSCAL( 1, SIGNBC, H( J, J ), 1 )
             END IF              END IF
             IF( ILZ )              IF( ILZ )
      $         CALL ZSCAL( N, SIGNBC, Z( 1, J ), 1 )       $         CALL ZSCAL( N, SIGNBC, Z( 1, J ), 1 )
Line 515 Line 516
          IF( ILAST.EQ.ILO ) THEN           IF( ILAST.EQ.ILO ) THEN
             GO TO 60              GO TO 60
          ELSE           ELSE
             IF( ABS1( H( ILAST, ILAST-1 ) ).LE.ATOL ) THEN              IF( ABS1( H( ILAST, ILAST-1 ) ).LE.MAX( SAFMIN, ULP*( 
        $         ABS1( H( ILAST, ILAST ) ) + ABS1( H( ILAST-1, ILAST-1 ) 
        $         ) ) ) ) THEN
                H( ILAST, ILAST-1 ) = CZERO                 H( ILAST, ILAST-1 ) = CZERO
                GO TO 60                 GO TO 60
             END IF              END IF
Line 535 Line 538
             IF( J.EQ.ILO ) THEN              IF( J.EQ.ILO ) THEN
                ILAZRO = .TRUE.                 ILAZRO = .TRUE.
             ELSE              ELSE
                IF( ABS1( H( J, J-1 ) ).LE.ATOL ) THEN                 IF( ABS1( H( J, J-1 ) ).LE.MAX( SAFMIN, ULP*( 
        $            ABS1( H( J, J ) ) + ABS1( H( J-1, J-1 ) ) 
        $            ) ) ) THEN
                   H( J, J-1 ) = CZERO                    H( J, J-1 ) = CZERO
                   ILAZRO = .TRUE.                    ILAZRO = .TRUE.
                ELSE                 ELSE
Line 666 Line 671
                CALL ZSCAL( ILAST+1-IFRSTM, SIGNBC, H( IFRSTM, ILAST ),                 CALL ZSCAL( ILAST+1-IFRSTM, SIGNBC, H( IFRSTM, ILAST ),
      $                     1 )       $                     1 )
             ELSE              ELSE
                H( ILAST, ILAST ) = H( ILAST, ILAST )*SIGNBC                 CALL ZSCAL( 1, SIGNBC, H( ILAST, ILAST ), 1 )
             END IF              END IF
             IF( ILZ )              IF( ILZ )
      $         CALL ZSCAL( N, SIGNBC, Z( 1, ILAST ), 1 )       $         CALL ZSCAL( N, SIGNBC, Z( 1, ILAST ), 1 )
Line 730 Line 735
             AD22 = ( ASCALE*H( ILAST, ILAST ) ) /              AD22 = ( ASCALE*H( ILAST, ILAST ) ) /
      $             ( BSCALE*T( ILAST, ILAST ) )       $             ( BSCALE*T( ILAST, ILAST ) )
             ABI22 = AD22 - U12*AD21              ABI22 = AD22 - U12*AD21
               ABI12 = AD12 - U12*AD11
 *  *
             T1 = HALF*( AD11+ABI22 )              SHIFT = ABI22
             RTDISC = SQRT( T1**2+AD12*AD21-AD11*AD22 )              CTEMP = SQRT( ABI12 )*SQRT( AD21 )
             TEMP = DBLE( T1-ABI22 )*DBLE( RTDISC ) +              TEMP = ABS1( CTEMP )
      $             DIMAG( T1-ABI22 )*DIMAG( RTDISC )              IF( CTEMP.NE.ZERO ) THEN
             IF( TEMP.LE.ZERO ) THEN                 X = HALF*( AD11-SHIFT )
                SHIFT = T1 + RTDISC                 TEMP2 = ABS1( X )
             ELSE                 TEMP = MAX( TEMP, ABS1( X ) )
                SHIFT = T1 - RTDISC                 Y = TEMP*SQRT( ( X / TEMP )**2+( CTEMP / TEMP )**2 )
                  IF( TEMP2.GT.ZERO ) THEN
                     IF( DBLE( X / TEMP2 )*DBLE( Y )+
        $                DIMAG( X / TEMP2 )*DIMAG( Y ).LT.ZERO )Y = -Y
                  END IF
                  SHIFT = SHIFT - CTEMP*ZLADIV( CTEMP, ( X+Y ) )
             END IF              END IF
          ELSE           ELSE
 *  *
 *           Exceptional shift.  Chosen for no particularly good reason.  *           Exceptional shift.  Chosen for no particularly good reason.
 *  *
             ESHIFT = ESHIFT + DCONJG( ( ASCALE*H( ILAST-1, ILAST ) ) /              IF( ( IITER / 20 )*20.EQ.IITER .AND. 
      $               ( BSCALE*T( ILAST-1, ILAST-1 ) ) )       $         BSCALE*ABS1(T( ILAST, ILAST )).GT.SAFMIN ) THEN
                  ESHIFT = ESHIFT + ( ASCALE*H( ILAST,
        $            ILAST ) )/( BSCALE*T( ILAST, ILAST ) )
               ELSE
                  ESHIFT = ESHIFT + ( ASCALE*H( ILAST,
        $            ILAST-1 ) )/( BSCALE*T( ILAST-1, ILAST-1 ) )
               END IF
             SHIFT = ESHIFT              SHIFT = ESHIFT
          END IF           END IF
 *  *
Line 850 Line 867
                CALL ZSCAL( J-1, SIGNBC, T( 1, J ), 1 )                 CALL ZSCAL( J-1, SIGNBC, T( 1, J ), 1 )
                CALL ZSCAL( J, SIGNBC, H( 1, J ), 1 )                 CALL ZSCAL( J, SIGNBC, H( 1, J ), 1 )
             ELSE              ELSE
                H( J, J ) = H( J, J )*SIGNBC                 CALL ZSCAL( 1, SIGNBC, H( J, J ), 1 )
             END IF              END IF
             IF( ILZ )              IF( ILZ )
      $         CALL ZSCAL( N, SIGNBC, Z( 1, J ), 1 )       $         CALL ZSCAL( N, SIGNBC, Z( 1, J ), 1 )
Removed from v.1.8  
changed lines
  Added in v.1.20

CVSweb interface <joel.bertrand@systella.fr>