dlaqz2.f

Go to the documentation of this file.

*> \brief \b DLAQZ2
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download DLAQZ2 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqz2.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqz2.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqz2.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*      SUBROUTINE DLAQZ2( ILQ, ILZ, K, ISTARTM, ISTOPM, IHI, A, LDA, B,
*     $    LDB, NQ, QSTART, Q, LDQ, NZ, ZSTART, Z, LDZ )
*      IMPLICIT NONE
*
*      Arguments
*      LOGICAL, INTENT( IN ) :: ILQ, ILZ
*      INTEGER, INTENT( IN ) :: K, LDA, LDB, LDQ, LDZ, ISTARTM, ISTOPM,
*     $    NQ, NZ, QSTART, ZSTART, IHI
*      DOUBLE PRECISION :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ), Z( LDZ,
*     $    * )
*      ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*>      DLAQZ2 chases a 2x2 shift bulge in a matrix pencil down a single position
*> \endverbatim
*
*
*  Arguments:
*  ==========
*
*>
*> \param[in] ILQ
*> \verbatim
*>          ILQ is LOGICAL
*>              Determines whether or not to update the matrix Q
*> \endverbatim
*>
*> \param[in] ILZ
*> \verbatim
*>          ILZ is LOGICAL
*>              Determines whether or not to update the matrix Z
*> \endverbatim
*>
*> \param[in] K
*> \verbatim
*>          K is INTEGER
*>              Index indicating the position of the bulge.
*>              On entry, the bulge is located in
*>              (A(k+1:k+2,k:k+1),B(k+1:k+2,k:k+1)).
*>              On exit, the bulge is located in
*>              (A(k+2:k+3,k+1:k+2),B(k+2:k+3,k+1:k+2)).
*> \endverbatim
*>
*> \param[in] ISTARTM
*> \verbatim
*>          ISTARTM is INTEGER
*> \endverbatim
*>
*> \param[in] ISTOPM
*> \verbatim
*>          ISTOPM is INTEGER
*>              Updates to (A,B) are restricted to
*>              (istartm:k+3,k:istopm). It is assumed
*>              without checking that istartm <= k+1 and
*>              k+2 <= istopm
*> \endverbatim
*>
*> \param[in] IHI
*> \verbatim
*>          IHI is INTEGER
*> \endverbatim
*>
*> \param[inout] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (LDA,N)
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>              The leading dimension of A as declared in
*>              the calling procedure.
*> \endverbatim
*
*> \param[inout] B
*> \verbatim
*>          B is DOUBLE PRECISION array, dimension (LDB,N)
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>              The leading dimension of B as declared in
*>              the calling procedure.
*> \endverbatim
*>
*> \param[in] NQ
*> \verbatim
*>          NQ is INTEGER
*>              The order of the matrix Q
*> \endverbatim
*>
*> \param[in] QSTART
*> \verbatim
*>          QSTART is INTEGER
*>              Start index of the matrix Q. Rotations are applied
*>              To columns k+2-qStart:k+4-qStart of Q.
*> \endverbatim
*
*> \param[inout] Q
*> \verbatim
*>          Q is DOUBLE PRECISION array, dimension (LDQ,NQ)
*> \endverbatim
*>
*> \param[in] LDQ
*> \verbatim
*>          LDQ is INTEGER
*>              The leading dimension of Q as declared in
*>              the calling procedure.
*> \endverbatim
*>
*> \param[in] NZ
*> \verbatim
*>          NZ is INTEGER
*>              The order of the matrix Z
*> \endverbatim
*>
*> \param[in] ZSTART
*> \verbatim
*>          ZSTART is INTEGER
*>              Start index of the matrix Z. Rotations are applied
*>              To columns k+1-qStart:k+3-qStart of Z.
*> \endverbatim
*
*> \param[inout] Z
*> \verbatim
*>          Z is DOUBLE PRECISION array, dimension (LDZ,NZ)
*> \endverbatim
*>
*> \param[in] LDZ
*> \verbatim
*>          LDZ is INTEGER
*>              The leading dimension of Q as declared in
*>              the calling procedure.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Thijs Steel, KU Leuven
*
*> \date May 2020
*
*> \ingroup laqz2
*>
*  =====================================================================
      SUBROUTINE dlaqz2( ILQ, ILZ, K, ISTARTM, ISTOPM, IHI, A, LDA, B,
     $                   LDB, NQ, QSTART, Q, LDQ, NZ, ZSTART, Z, LDZ )
      IMPLICIT NONE
*
*     Arguments
      LOGICAL, INTENT( IN ) :: ILQ, ILZ
      INTEGER, INTENT( IN ) :: K, LDA, LDB, LDQ, LDZ, ISTARTM, ISTOPM,
     $         nq, nz, qstart, zstart, ihi
      DOUBLE PRECISION :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ), Z( LDZ,
     $                    * )
*
*     Parameters
      DOUBLE PRECISION :: ZERO, ONE, HALF
      parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
*
*     Local variables
      DOUBLE PRECISION :: H( 2, 3 ), C1, S1, C2, S2, TEMP
*
*     External functions
      EXTERNAL :: dlartg, drot
*
      IF( k+2 .EQ. ihi ) THEN
*        Shift is located on the edge of the matrix, remove it
         h = b( ihi-1:ihi, ihi-2:ihi )
*        Make H upper triangular
         CALL dlartg( h( 1, 1 ), h( 2, 1 ), c1, s1, temp )
         h( 2, 1 ) = zero
         h( 1, 1 ) = temp
         CALL drot( 2, h( 1, 2 ), 2, h( 2, 2 ), 2, c1, s1 )
*
         CALL dlartg( h( 2, 3 ), h( 2, 2 ), c1, s1, temp )
         CALL drot( 1, h( 1, 3 ), 1, h( 1, 2 ), 1, c1, s1 )
         CALL dlartg( h( 1, 2 ), h( 1, 1 ), c2, s2, temp )
*
         CALL drot( ihi-istartm+1, b( istartm, ihi ), 1, b( istartm,
     $              ihi-1 ), 1, c1, s1 )
         CALL drot( ihi-istartm+1, b( istartm, ihi-1 ), 1, b( istartm,
     $              ihi-2 ), 1, c2, s2 )
         b( ihi-1, ihi-2 ) = zero
         b( ihi, ihi-2 ) = zero
         CALL drot( ihi-istartm+1, a( istartm, ihi ), 1, a( istartm,
     $              ihi-1 ), 1, c1, s1 )
         CALL drot( ihi-istartm+1, a( istartm, ihi-1 ), 1, a( istartm,
     $              ihi-2 ), 1, c2, s2 )
         IF ( ilz ) THEN
            CALL drot( nz, z( 1, ihi-zstart+1 ), 1, z( 1, ihi-1-zstart+
     $                 1 ), 1, c1, s1 )
            CALL drot( nz, z( 1, ihi-1-zstart+1 ), 1, z( 1,
     $                 ihi-2-zstart+1 ), 1, c2, s2 )
         END IF
*
         CALL dlartg( a( ihi-1, ihi-2 ), a( ihi, ihi-2 ), c1, s1,
     $                temp )
         a( ihi-1, ihi-2 ) = temp
         a( ihi, ihi-2 ) = zero
         CALL drot( istopm-ihi+2, a( ihi-1, ihi-1 ), lda, a( ihi,
     $              ihi-1 ), lda, c1, s1 )
         CALL drot( istopm-ihi+2, b( ihi-1, ihi-1 ), ldb, b( ihi,
     $              ihi-1 ), ldb, c1, s1 )
         IF ( ilq ) THEN
            CALL drot( nq, q( 1, ihi-1-qstart+1 ), 1, q( 1, ihi-qstart+
     $                 1 ), 1, c1, s1 )
         END IF
*
         CALL dlartg( b( ihi, ihi ), b( ihi, ihi-1 ), c1, s1, temp )
         b( ihi, ihi ) = temp
         b( ihi, ihi-1 ) = zero
         CALL drot( ihi-istartm, b( istartm, ihi ), 1, b( istartm,
     $              ihi-1 ), 1, c1, s1 )
         CALL drot( ihi-istartm+1, a( istartm, ihi ), 1, a( istartm,
     $              ihi-1 ), 1, c1, s1 )
         IF ( ilz ) THEN
            CALL drot( nz, z( 1, ihi-zstart+1 ), 1, z( 1, ihi-1-zstart+
     $                 1 ), 1, c1, s1 )
         END IF
*
      ELSE
*
*        Normal operation, move bulge down
*
         h = b( k+1:k+2, k:k+2 )
*
*        Make H upper triangular
*
         CALL dlartg( h( 1, 1 ), h( 2, 1 ), c1, s1, temp )
         h( 2, 1 ) = zero
         h( 1, 1 ) = temp
         CALL drot( 2, h( 1, 2 ), 2, h( 2, 2 ), 2, c1, s1 )
*
*        Calculate Z1 and Z2
*
         CALL dlartg( h( 2, 3 ), h( 2, 2 ), c1, s1, temp )
         CALL drot( 1, h( 1, 3 ), 1, h( 1, 2 ), 1, c1, s1 )
         CALL dlartg( h( 1, 2 ), h( 1, 1 ), c2, s2, temp )
*
*        Apply transformations from the right
*
         CALL drot( k+3-istartm+1, a( istartm, k+2 ), 1, a( istartm,
     $              k+1 ), 1, c1, s1 )
         CALL drot( k+3-istartm+1, a( istartm, k+1 ), 1, a( istartm,
     $              k ), 1, c2, s2 )
         CALL drot( k+2-istartm+1, b( istartm, k+2 ), 1, b( istartm,
     $              k+1 ), 1, c1, s1 )
         CALL drot( k+2-istartm+1, b( istartm, k+1 ), 1, b( istartm,
     $              k ), 1, c2, s2 )
         IF ( ilz ) THEN
            CALL drot( nz, z( 1, k+2-zstart+1 ), 1, z( 1, k+1-zstart+
     $                 1 ), 1, c1, s1 )
            CALL drot( nz, z( 1, k+1-zstart+1 ), 1, z( 1, k-zstart+1 ),
     $                 1, c2, s2 )
         END IF
         b( k+1, k ) = zero
         b( k+2, k ) = zero
*
*        Calculate Q1 and Q2
*
         CALL dlartg( a( k+2, k ), a( k+3, k ), c1, s1, temp )
         a( k+2, k ) = temp
         a( k+3, k ) = zero
         CALL dlartg( a( k+1, k ), a( k+2, k ), c2, s2, temp )
         a( k+1, k ) = temp
         a( k+2, k ) = zero
*
*        Apply transformations from the left
*
         CALL drot( istopm-k, a( k+2, k+1 ), lda, a( k+3, k+1 ), lda,
     $              c1, s1 )
         CALL drot( istopm-k, a( k+1, k+1 ), lda, a( k+2, k+1 ), lda,
     $              c2, s2 )
*
         CALL drot( istopm-k, b( k+2, k+1 ), ldb, b( k+3, k+1 ), ldb,
     $              c1, s1 )
         CALL drot( istopm-k, b( k+1, k+1 ), ldb, b( k+2, k+1 ), ldb,
     $              c2, s2 )
         IF ( ilq ) THEN
            CALL drot( nq, q( 1, k+2-qstart+1 ), 1, q( 1, k+3-qstart+
     $                 1 ), 1, c1, s1 )
            CALL drot( nq, q( 1, k+1-qstart+1 ), 1, q( 1, k+2-qstart+
     $                 1 ), 1, c2, s2 )
         END IF
*
      END IF
*
*     End of DLAQZ2
*
      END SUBROUTINE