dlaqz4()

subroutine dlaqz4	(	logical, intent(in)	ilschur,
		logical, intent(in)	ilq,
		logical, intent(in)	ilz,
		integer, intent(in)	n,
		integer, intent(in)	ilo,
		integer, intent(in)	ihi,
		integer, intent(in)	nshifts,
		integer, intent(in)	nblock_desired,
		double precision, dimension( * ), intent(inout)	sr,
		double precision, dimension( * ), intent(inout)	si,
		double precision, dimension( * ), intent(inout)	ss,
		double precision, dimension( lda, * ), intent(inout)	a,
		integer, intent(in)	lda,
		double precision, dimension( ldb, * ), intent(inout)	b,
		integer, intent(in)	ldb,
		double precision, dimension( ldq, * ), intent(inout)	q,
		integer, intent(in)	ldq,
		double precision, dimension( ldz, * ), intent(inout)	z,
		integer, intent(in)	ldz,
		double precision, dimension( ldqc, * ), intent(inout)	qc,
		integer, intent(in)	ldqc,
		double precision, dimension( ldzc, * ), intent(inout)	zc,
		integer, intent(in)	ldzc,
		double precision, dimension( * ), intent(inout)	work,
		integer, intent(in)	lwork,
		integer, intent(out)	info
	)

DLAQZ4

Download DLAQZ4 + dependencies [TGZ] [ZIP] [TXT]

Purpose:

 DLAQZ4 Executes a single multishift QZ sweep

Parameters

[in]	ILSCHUR	ILSCHUR is LOGICAL Determines whether or not to update the full Schur form
[in]	ILQ	ILQ is LOGICAL Determines whether or not to update the matrix Q
[in]	ILZ	ILZ is LOGICAL Determines whether or not to update the matrix Z
[in]	N	N is INTEGER The order of the matrices A, B, Q, and Z. N >= 0.
[in]	ILO	ILO is INTEGER
[in]	IHI	IHI is INTEGER
[in]	NSHIFTS	NSHIFTS is INTEGER The desired number of shifts to use
[in]	NBLOCK_DESIRED	NBLOCK_DESIRED is INTEGER The desired size of the computational windows
[in]	SR	SR is DOUBLE PRECISION array. SR contains the real parts of the shifts to use.
[in]	SI	SI is DOUBLE PRECISION array. SI contains the imaginary parts of the shifts to use.
[in]	SS	SS is DOUBLE PRECISION array. SS contains the scale of the shifts to use.
[in,out]	A	A is DOUBLE PRECISION array, dimension (LDA, N)
[in]	LDA	LDA is INTEGER The leading dimension of the array A. LDA >= max( 1, N ).
[in,out]	B	B is DOUBLE PRECISION array, dimension (LDB, N)
[in]	LDB	LDB is INTEGER The leading dimension of the array B. LDB >= max( 1, N ).
[in,out]	Q	Q is DOUBLE PRECISION array, dimension (LDQ, N)
[in]	LDQ	LDQ is INTEGER
[in,out]	Z	Z is DOUBLE PRECISION array, dimension (LDZ, N)
[in]	LDZ	LDZ is INTEGER
[in,out]	QC	QC is DOUBLE PRECISION array, dimension (LDQC, NBLOCK_DESIRED)
[in]	LDQC	LDQC is INTEGER
[in,out]	ZC	ZC is DOUBLE PRECISION array, dimension (LDZC, NBLOCK_DESIRED)
[in]	LDZC	LDZ is INTEGER
[out]	WORK	WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
[in]	LWORK	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.
[out]	INFO	INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value

Author

Thijs Steel, KU Leuven

Date

May 2020

Definition at line 209 of file dlaqz4.f.

      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
 
      DOUBLE PRECISION, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ),
     $                  Q( LDQ, * ), Z( LDZ, * ), QC( LDQC, * ),
     $                  ZC( LDZC, * ), WORK( * ), SR( * ), SI( * ),
     $                  SS( * )
 
      INTEGER, INTENT( OUT ) :: INFO
 
*     Parameters
      DOUBLE PRECISION :: ZERO, ONE, HALF
      parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
 
*     Local scalars
      INTEGER :: I, J, NS, ISTARTM, ISTOPM, SHEIGHT, SWIDTH, K, NP,
     $           ISTARTB, ISTOPB, ISHIFT, NBLOCK, NPOS
      DOUBLE PRECISION :: TEMP, V( 3 ), C1, S1, C2, S2, SWAP
*
*     External functions
      EXTERNAL :: xerbla, dgemm, dlaqz1, dlaqz2, dlaset, dlartg, drot,
     $            dlacpy
 
      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( 'DLAQZ4', -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 dlaset( 'FULL', ns+1, ns+1, zero, one, qc, ldqc )
      CALL dlaset( 'FULL', ns, ns, zero, one, zc, ldzc )
 
      DO i = 1, ns, 2
*        Introduce the shift
         CALL dlaqz1( 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 dlartg( temp, v( 3 ), c1, s1, v( 2 ) )
         CALL dlartg( v( 1 ), v( 2 ), c2, s2, temp )
 
         CALL drot( ns, a( ilo+1, ilo ), lda, a( ilo+2, ilo ), lda, c1,
     $              s1 )
         CALL drot( ns, a( ilo, ilo ), lda, a( ilo+1, ilo ), lda, c2,
     $              s2 )
         CALL drot( ns, b( ilo+1, ilo ), ldb, b( ilo+2, ilo ), ldb, c1,
     $              s1 )
         CALL drot( ns, b( ilo, ilo ), ldb, b( ilo+1, ilo ), ldb, c2,
     $              s2 )
         CALL drot( ns+1, qc( 1, 2 ), 1, qc( 1, 3 ), 1, c1, s1 )
         CALL drot( ns+1, qc( 1, 1 ), 1, qc( 1, 2 ), 1, c2, s2 )
 
*        Chase the shift down
         DO j = 1, ns-1-i
 
            CALL dlaqz2( .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 dgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
     $               a( ilo, ilo+ns ), lda, zero, work, sheight )
         CALL dlacpy( 'ALL', sheight, swidth, work, sheight, a( ilo,
     $                ilo+ns ), lda )
         CALL dgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
     $               b( ilo, ilo+ns ), ldb, zero, work, sheight )
         CALL dlacpy( 'ALL', sheight, swidth, work, sheight, b( ilo,
     $                ilo+ns ), ldb )
      END IF
      IF ( ilq ) THEN
         CALL dgemm( 'N', 'N', n, sheight, sheight, one, q( 1, ilo ),
     $               ldq, qc, ldqc, zero, work, n )
         CALL dlacpy( '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 dgemm( 'N', 'N', sheight, swidth, swidth, one, a( istartm,
     $               ilo ), lda, zc, ldzc, zero, work, sheight )
         CALL dlacpy( 'ALL', sheight, swidth, work, sheight, a( istartm,
     $                ilo ), lda )
         CALL dgemm( 'N', 'N', sheight, swidth, swidth, one, b( istartm,
     $               ilo ), ldb, zc, ldzc, zero, work, sheight )
         CALL dlacpy( 'ALL', sheight, swidth, work, sheight, b( istartm,
     $                ilo ), ldb )
      END IF
      IF ( ilz ) THEN
         CALL dgemm( 'N', 'N', n, swidth, swidth, one, z( 1, ilo ), ldz,
     $               zc, ldzc, zero, work, n )
         CALL dlacpy( '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 dlaset( 'FULL', ns+np, ns+np, zero, one, qc, ldqc )
         CALL dlaset( '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 dlaqz2( .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 dgemm( 'T', 'N', sheight, swidth, sheight, one, qc,
     $                  ldqc, a( k+1, k+ns+np ), lda, zero, work,
     $                  sheight )
            CALL dlacpy( 'ALL', sheight, swidth, work, sheight, a( k+1,
     $                   k+ns+np ), lda )
            CALL dgemm( 'T', 'N', sheight, swidth, sheight, one, qc,
     $                  ldqc, b( k+1, k+ns+np ), ldb, zero, work,
     $                  sheight )
            CALL dlacpy( 'ALL', sheight, swidth, work, sheight, b( k+1,
     $                   k+ns+np ), ldb )
         END IF
         IF ( ilq ) THEN
            CALL dgemm( 'N', 'N', n, nblock, nblock, one, q( 1, k+1 ),
     $                  ldq, qc, ldqc, zero, work, n )
            CALL dlacpy( '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 dgemm( 'N', 'N', sheight, swidth, swidth, one,
     $                  a( istartm, k ), lda, zc, ldzc, zero, work,
     $                  sheight )
            CALL dlacpy( 'ALL', sheight, swidth, work, sheight,
     $                   a( istartm, k ), lda )
            CALL dgemm( 'N', 'N', sheight, swidth, swidth, one,
     $                  b( istartm, k ), ldb, zc, ldzc, zero, work,
     $                  sheight )
            CALL dlacpy( 'ALL', sheight, swidth, work, sheight,
     $                   b( istartm, k ), ldb )
         END IF
         IF ( ilz ) THEN
            CALL dgemm( 'N', 'N', n, nblock, nblock, one, z( 1, k ),
     $                  ldz, zc, ldzc, zero, work, n )
            CALL dlacpy( '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 dlaset( 'FULL', ns, ns, zero, one, qc, ldqc )
      CALL dlaset( '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 dlaqz2( .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 dgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
     $               a( ihi-ns+1, ihi+1 ), lda, zero, work, sheight )
         CALL dlacpy( 'ALL', sheight, swidth, work, sheight,
     $                a( ihi-ns+1, ihi+1 ), lda )
         CALL dgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
     $               b( ihi-ns+1, ihi+1 ), ldb, zero, work, sheight )
         CALL dlacpy( 'ALL', sheight, swidth, work, sheight,
     $                b( ihi-ns+1, ihi+1 ), ldb )
      END IF
      IF ( ilq ) THEN
         CALL dgemm( 'N', 'N', n, ns, ns, one, q( 1, ihi-ns+1 ), ldq,
     $               qc, ldqc, zero, work, n )
         CALL dlacpy( '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 dgemm( 'N', 'N', sheight, swidth, swidth, one, a( istartm,
     $               ihi-ns ), lda, zc, ldzc, zero, work, sheight )
         CALL dlacpy( 'ALL', sheight, swidth, work, sheight, a( istartm,
     $                ihi-ns ), lda )
         CALL dgemm( 'N', 'N', sheight, swidth, swidth, one, b( istartm,
     $               ihi-ns ), ldb, zc, ldzc, zero, work, sheight )
         CALL dlacpy( 'ALL', sheight, swidth, work, sheight, b( istartm,
     $                ihi-ns ), ldb )
      END IF
      IF ( ilz ) THEN
         CALL dgemm( 'N', 'N', n, ns+1, ns+1, one, z( 1, ihi-ns ), ldz,
     $               zc, ldzc, zero, work, n )
         CALL dlacpy( 'ALL', n, ns+1, work, n, z( 1, ihi-ns ), ldz )
      END IF
 

Here is the call graph for this function:

Here is the caller graph for this function: