slaqz4.f

Go to the documentation of this file.

*> \brief \b SLAQZ4
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download SLAQZ4 + dependencies
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaqz4.f">
*> [TGZ]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaqz4.f">
*> [ZIP]</a>
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaqz4.f">
*> [TXT]</a>
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*      SUBROUTINE SLAQZ4( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS,
*     $    NBLOCK_DESIRED, SR, SI, SS, A, LDA, B, LDB, Q, LDQ, Z, LDZ,
*     $    QC, LDQC, ZC, LDZC, WORK, LWORK, INFO )
*      IMPLICIT NONE
*
*      Function arguments
*      LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
*      INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
*     $    NSHIFTS, NBLOCK_DESIRED, LDQC, LDZC
*
*      REAL, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
*     $    Z( LDZ, * ), QC( LDQC, * ), ZC( LDZC, * ), WORK( * ), SR( * ),
*     $    SI( * ), SS( * )
*
*      INTEGER, INTENT( OUT ) :: INFO
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SLAQZ4 Executes a single multishift QZ sweep
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] ILSCHUR
*> \verbatim
*>          ILSCHUR is LOGICAL
*>              Determines whether or not to update the full Schur form
*> \endverbatim
*>
*> \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] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrices A, B, Q, and Z.  N >= 0.
*> \endverbatim
*>
*> \param[in] ILO
*> \verbatim
*>          ILO is INTEGER
*> \endverbatim
*>
*> \param[in] IHI
*> \verbatim
*>          IHI is INTEGER
*> \endverbatim
*>
*> \param[in] NSHIFTS
*> \verbatim
*>          NSHIFTS is INTEGER
*>          The desired number of shifts to use
*> \endverbatim
*>
*> \param[in] NBLOCK_DESIRED
*> \verbatim
*>          NBLOCK_DESIRED is INTEGER
*>          The desired size of the computational windows
*> \endverbatim
*>
*> \param[in] SR
*> \verbatim
*>          SR is REAL array. SR contains
*>          the real parts of the shifts to use.
*> \endverbatim
*>
*> \param[in] SI
*> \verbatim
*>          SI is REAL array. SI contains
*>          the imaginary parts of the shifts to use.
*> \endverbatim
*>
*> \param[in] SS
*> \verbatim
*>          SS is REAL array. SS contains
*>          the scale of the shifts to use.
*> \endverbatim
*>
*> \param[in,out] A
*> \verbatim
*>          A is REAL array, dimension (LDA, N)
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A.  LDA >= max( 1, N ).
*> \endverbatim
*>
*> \param[in,out] B
*> \verbatim
*>          B is REAL array, dimension (LDB, N)
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>          The leading dimension of the array B.  LDB >= max( 1, N ).
*> \endverbatim
*>
*> \param[in,out] Q
*> \verbatim
*>          Q is REAL array, dimension (LDQ, N)
*> \endverbatim
*>
*> \param[in] LDQ
*> \verbatim
*>          LDQ is INTEGER
*> \endverbatim
*>
*> \param[in,out] Z
*> \verbatim
*>          Z is REAL array, dimension (LDZ, N)
*> \endverbatim
*>
*> \param[in] LDZ
*> \verbatim
*>          LDZ is INTEGER
*> \endverbatim
*>
*> \param[in,out] QC
*> \verbatim
*>          QC is REAL array, dimension (LDQC, NBLOCK_DESIRED)
*> \endverbatim
*>
*> \param[in] LDQC
*> \verbatim
*>          LDQC is INTEGER
*> \endverbatim
*>
*> \param[in,out] ZC
*> \verbatim
*>          ZC is REAL array, dimension (LDZC, NBLOCK_DESIRED)
*> \endverbatim
*>
*> \param[in] LDZC
*> \verbatim
*>          LDZ is INTEGER
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is REAL 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
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Thijs Steel, KU Leuven
*
*> \date May 2020
*
*> \ingroup doubleGEcomputational
*>
*  =====================================================================
      SUBROUTINE slaqz4( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS,
     $                   NBLOCK_DESIRED, SR, SI, SS, A, LDA, B, LDB, Q,
     $                   LDQ, Z, LDZ, QC, LDQC, ZC, LDZC, WORK, LWORK,
     $                   INFO )
      IMPLICIT NONE
 
*     Function arguments
      LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
      INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
     $         NSHIFTS, NBLOCK_DESIRED, LDQC, LDZC
 
      REAL, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
     $   Z( LDZ, * ), QC( LDQC, * ), ZC( LDZC, * ), WORK( * ), SR( * ),
     $   SI( * ), SS( * )
 
      INTEGER, INTENT( OUT ) :: INFO
 
*     Parameters
      REAL :: ZERO, ONE, HALF
      PARAMETER( ZERO = 0.0, one = 1.0, half = 0.5 )
 
*     Local scalars
      INTEGER :: I, J, NS, ISTARTM, ISTOPM, SHEIGHT, SWIDTH, K, NP,
     $           ISTARTB, ISTOPB, ISHIFT, NBLOCK, NPOS
      REAL :: TEMP, V( 3 ), C1, S1, C2, S2, SWAP
*
*     External functions
      EXTERNAL :: xerbla, sgemm, slaqz1, slaqz2, slaset, slartg, srot,
     $            slacpy
 
      info = 0
      IF ( nblock_desired .LT. nshifts+1 ) THEN
         info = -8
      END IF
      IF ( lwork .EQ.-1 ) THEN
*        workspace query, quick return
         work( 1 ) = n*nblock_desired
         RETURN
      ELSE IF ( lwork .LT. n*nblock_desired ) THEN
         info = -25
      END IF
 
      IF( info.NE.0 ) THEN
         CALL xerbla( 'SLAQZ4', -info )
         RETURN
      END IF
 
*     Executable statements
 
      IF ( nshifts .LT. 2 ) THEN
         RETURN
      END IF
 
      IF ( ilo .GE. ihi ) THEN
         RETURN
      END IF
 
      IF ( ilschur ) THEN
         istartm = 1
         istopm = n
      ELSE
         istartm = ilo
         istopm = ihi
      END IF
 
*     Shuffle shifts into pairs of real shifts and pairs
*     of complex conjugate shifts assuming complex
*     conjugate shifts are already adjacent to one
*     another
 
      DO i = 1, nshifts-2, 2
         IF( si( i ).NE.-si( i+1 ) ) THEN
*
            swap = sr( i )
            sr( i ) = sr( i+1 )
            sr( i+1 ) = sr( i+2 )
            sr( i+2 ) = swap
 
            swap = si( i )
            si( i ) = si( i+1 )
            si( i+1 ) = si( i+2 )
            si( i+2 ) = swap
            
            swap = ss( i )
            ss( i ) = ss( i+1 )
            ss( i+1 ) = ss( i+2 )
            ss( i+2 ) = swap
         END IF
      END DO
 
*     NSHFTS is supposed to be even, but if it is odd,
*     then simply reduce it by one.  The shuffle above
*     ensures that the dropped shift is real and that
*     the remaining shifts are paired.
 
      ns = nshifts-mod( nshifts, 2 )
      npos = max( nblock_desired-ns, 1 )
 
*     The following block introduces the shifts and chases
*     them down one by one just enough to make space for
*     the other shifts. The near-the-diagonal block is
*     of size (ns+1) x ns.
 
      CALL slaset( 'FULL', ns+1, ns+1, zero, one, qc, ldqc )
      CALL slaset( 'FULL', ns, ns, zero, one, zc, ldzc )
 
      DO i = 1, ns, 2
*        Introduce the shift
         CALL slaqz1( a( ilo, ilo ), lda, b( ilo, ilo ), ldb, sr( i ),
     $                sr( i+1 ), si( i ), ss( i ), ss( i+1 ), v )
 
         temp = v( 2 )
         CALL slartg( temp, v( 3 ), c1, s1, v( 2 ) )
         CALL slartg( v( 1 ), v( 2 ), c2, s2, temp )
 
         CALL srot( ns, a( ilo+1, ilo ), lda, a( ilo+2, ilo ), lda, c1,
     $              s1 )
         CALL srot( ns, a( ilo, ilo ), lda, a( ilo+1, ilo ), lda, c2,
     $              s2 )
         CALL srot( ns, b( ilo+1, ilo ), ldb, b( ilo+2, ilo ), ldb, c1,
     $              s1 )
         CALL srot( ns, b( ilo, ilo ), ldb, b( ilo+1, ilo ), ldb, c2,
     $              s2 )
         CALL srot( ns+1, qc( 1, 2 ), 1, qc( 1, 3 ), 1, c1, s1 )
         CALL srot( ns+1, qc( 1, 1 ), 1, qc( 1, 2 ), 1, c2, s2 )
 
*        Chase the shift down
         DO j = 1, ns-1-i
 
            CALL slaqz2( .true., .true., j, 1, ns, ihi-ilo+1, a( ilo,
     $                   ilo ), lda, b( ilo, ilo ), ldb, ns+1, 1, qc,
     $                   ldqc, ns, 1, zc, ldzc )
 
         END DO
 
      END DO
 
*     Update the rest of the pencil
 
*     Update A(ilo:ilo+ns,ilo+ns:istopm) and B(ilo:ilo+ns,ilo+ns:istopm)
*     from the left with Qc(1:ns+1,1:ns+1)'
      sheight = ns+1
      swidth = istopm-( ilo+ns )+1
      IF ( swidth > 0 ) THEN
         CALL sgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
     $               a( ilo, ilo+ns ), lda, zero, work, sheight )
         CALL slacpy( 'ALL', sheight, swidth, work, sheight, a( ilo,
     $                ilo+ns ), lda )
         CALL sgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
     $               b( ilo, ilo+ns ), ldb, zero, work, sheight )
         CALL slacpy( 'ALL', sheight, swidth, work, sheight, b( ilo,
     $                ilo+ns ), ldb )
      END IF
      IF ( ilq ) THEN
         CALL sgemm( 'N', 'N', n, sheight, sheight, one, q( 1, ilo ),
     $               ldq, qc, ldqc, zero, work, n )
         CALL slacpy( 'ALL', n, sheight, work, n, q( 1, ilo ), ldq )
      END IF
 
*     Update A(istartm:ilo-1,ilo:ilo+ns-1) and B(istartm:ilo-1,ilo:ilo+ns-1)
*     from the right with Zc(1:ns,1:ns)
      sheight = ilo-1-istartm+1
      swidth = ns
      IF ( sheight > 0 ) THEN
         CALL sgemm( 'N', 'N', sheight, swidth, swidth, one, a( istartm,
     $               ilo ), lda, zc, ldzc, zero, work, sheight )
         CALL slacpy( 'ALL', sheight, swidth, work, sheight, a( istartm,
     $                ilo ), lda )
         CALL sgemm( 'N', 'N', sheight, swidth, swidth, one, b( istartm,
     $               ilo ), ldb, zc, ldzc, zero, work, sheight )
         CALL slacpy( 'ALL', sheight, swidth, work, sheight, b( istartm,
     $                ilo ), ldb )
      END IF
      IF ( ilz ) THEN
         CALL sgemm( 'N', 'N', n, swidth, swidth, one, z( 1, ilo ), ldz,
     $               zc, ldzc, zero, work, n )
         CALL slacpy( 'ALL', n, swidth, work, n, z( 1, ilo ), ldz )
      END IF
 
*     The following block chases the shifts down to the bottom
*     right block. If possible, a shift is moved down npos
*     positions at a time
 
      k = ilo
      DO WHILE ( k < ihi-ns )
         np = min( ihi-ns-k, npos )
*        Size of the near-the-diagonal block
         nblock = ns+np
*        istartb points to the first row we will be updating
         istartb = k+1
*        istopb points to the last column we will be updating
         istopb = k+nblock-1
 
         CALL slaset( 'FULL', ns+np, ns+np, zero, one, qc, ldqc )
         CALL slaset( 'FULL', ns+np, ns+np, zero, one, zc, ldzc )
 
*        Near the diagonal shift chase
         DO i = ns-1, 0, -2
            DO j = 0, np-1
*              Move down the block with index k+i+j-1, updating
*              the (ns+np x ns+np) block:
*              (k:k+ns+np,k:k+ns+np-1)
               CALL slaqz2( .true., .true., k+i+j-1, istartb, istopb,
     $                      ihi, a, lda, b, ldb, nblock, k+1, qc, ldqc,
     $                      nblock, k, zc, ldzc )
            END DO
         END DO
 
*        Update rest of the pencil
 
*        Update A(k+1:k+ns+np, k+ns+np:istopm) and
*        B(k+1:k+ns+np, k+ns+np:istopm)
*        from the left with Qc(1:ns+np,1:ns+np)'
         sheight = ns+np
         swidth = istopm-( k+ns+np )+1
         IF ( swidth > 0 ) THEN
            CALL sgemm( 'T', 'N', sheight, swidth, sheight, one, qc,
     $                  ldqc, a( k+1, k+ns+np ), lda, zero, work,
     $                  sheight )
            CALL slacpy( 'ALL', sheight, swidth, work, sheight, a( k+1,
     $                   k+ns+np ), lda )
            CALL sgemm( 'T', 'N', sheight, swidth, sheight, one, qc,
     $                  ldqc, b( k+1, k+ns+np ), ldb, zero, work,
     $                  sheight )
            CALL slacpy( 'ALL', sheight, swidth, work, sheight, b( k+1,
     $                   k+ns+np ), ldb )
         END IF
         IF ( ilq ) THEN
            CALL sgemm( 'N', 'N', n, nblock, nblock, one, q( 1, k+1 ),
     $                  ldq, qc, ldqc, zero, work, n )
            CALL slacpy( 'ALL', n, nblock, work, n, q( 1, k+1 ), ldq )
         END IF
 
*        Update A(istartm:k,k:k+ns+npos-1) and B(istartm:k,k:k+ns+npos-1)
*        from the right with Zc(1:ns+np,1:ns+np)
         sheight = k-istartm+1
         swidth = nblock
         IF ( sheight > 0 ) THEN
            CALL sgemm( 'N', 'N', sheight, swidth, swidth, one,
     $                  a( istartm, k ), lda, zc, ldzc, zero, work,
     $                  sheight )
            CALL slacpy( 'ALL', sheight, swidth, work, sheight,
     $                   a( istartm, k ), lda )
            CALL sgemm( 'N', 'N', sheight, swidth, swidth, one,
     $                  b( istartm, k ), ldb, zc, ldzc, zero, work,
     $                  sheight )
            CALL slacpy( 'ALL', sheight, swidth, work, sheight,
     $                   b( istartm, k ), ldb )
         END IF
         IF ( ilz ) THEN
            CALL sgemm( 'N', 'N', n, nblock, nblock, one, z( 1, k ),
     $                  ldz, zc, ldzc, zero, work, n )
            CALL slacpy( 'ALL', n, nblock, work, n, z( 1, k ), ldz )
         END IF
 
         k = k+np
 
      END DO
 
*     The following block removes the shifts from the bottom right corner
*     one by one. Updates are initially applied to A(ihi-ns+1:ihi,ihi-ns:ihi).
 
      CALL slaset( 'FULL', ns, ns, zero, one, qc, ldqc )
      CALL slaset( 'FULL', ns+1, ns+1, zero, one, zc, ldzc )
 
*     istartb points to the first row we will be updating
      istartb = ihi-ns+1
*     istopb points to the last column we will be updating
      istopb = ihi
 
      DO i = 1, ns, 2
*        Chase the shift down to the bottom right corner
         DO ishift = ihi-i-1, ihi-2
            CALL slaqz2( .true., .true., ishift, istartb, istopb, ihi,
     $                   a, lda, b, ldb, ns, ihi-ns+1, qc, ldqc, ns+1,
     $                   ihi-ns, zc, ldzc )
         END DO
         
      END DO
 
*     Update rest of the pencil
 
*     Update A(ihi-ns+1:ihi, ihi+1:istopm)
*     from the left with Qc(1:ns,1:ns)'
      sheight = ns
      swidth = istopm-( ihi+1 )+1
      IF ( swidth > 0 ) THEN
         CALL sgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
     $               a( ihi-ns+1, ihi+1 ), lda, zero, work, sheight )
         CALL slacpy( 'ALL', sheight, swidth, work, sheight,
     $                a( ihi-ns+1, ihi+1 ), lda )
         CALL sgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
     $               b( ihi-ns+1, ihi+1 ), ldb, zero, work, sheight )
         CALL slacpy( 'ALL', sheight, swidth, work, sheight,
     $                b( ihi-ns+1, ihi+1 ), ldb )
      END IF
      IF ( ilq ) THEN
         CALL sgemm( 'N', 'N', n, ns, ns, one, q( 1, ihi-ns+1 ), ldq,
     $               qc, ldqc, zero, work, n )
         CALL slacpy( 'ALL', n, ns, work, n, q( 1, ihi-ns+1 ), ldq )
      END IF
 
*     Update A(istartm:ihi-ns,ihi-ns:ihi)
*     from the right with Zc(1:ns+1,1:ns+1)
      sheight = ihi-ns-istartm+1
      swidth = ns+1
      IF ( sheight > 0 ) THEN
         CALL sgemm( 'N', 'N', sheight, swidth, swidth, one, a( istartm,
     $               ihi-ns ), lda, zc, ldzc, zero, work, sheight )
         CALL slacpy( 'ALL', sheight, swidth, work, sheight, a( istartm,
     $                ihi-ns ), lda )
         CALL sgemm( 'N', 'N', sheight, swidth, swidth, one, b( istartm,
     $               ihi-ns ), ldb, zc, ldzc, zero, work, sheight )
         CALL slacpy( 'ALL', sheight, swidth, work, sheight, b( istartm,
     $                ihi-ns ), ldb )
      END IF
      IF ( ilz ) THEN
      CALL sgemm( 'N', 'N', n, ns+1, ns+1, one, z( 1, ihi-ns ), ldz, zc,
     $            ldzc, zero, work, n )
         CALL slacpy( 'ALL', n, ns+1, work, n, z( 1, ihi-ns ), ldz )
      END IF
 
      END SUBROUTINE