zlaqz3()

subroutine zlaqz3	(	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,
		complex*16, dimension( * ), intent(inout)	ALPHA,
		complex*16, dimension( * ), intent(inout)	BETA,
		complex*16, dimension( lda, * ), intent(inout)	A,
		integer, intent(in)	LDA,
		complex*16, dimension( ldb, * ), intent(inout)	B,
		integer, intent(in)	LDB,
		complex*16, dimension( ldq, * ), intent(inout)	Q,
		integer, intent(in)	LDQ,
		complex*16, dimension( ldz, * ), intent(inout)	Z,
		integer, intent(in)	LDZ,
		complex*16, dimension( ldqc, * ), intent(inout)	QC,
		integer, intent(in)	LDQC,
		complex*16, dimension( ldzc, * ), intent(inout)	ZC,
		integer, intent(in)	LDZC,
		complex*16, dimension( * ), intent(inout)	WORK,
		integer, intent(in)	LWORK,
		integer, intent(out)	INFO
	)

ZLAQZ3

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

Purpose:

 ZLAQZ3 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]	ALPHA	ALPHA is COMPLEX*16 array. SR contains the alpha parts of the shifts to use.
[in]	BETA	BETA is COMPLEX*16 array. SS contains the scale of the shifts to use.
[in,out]	A	A is COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 array, dimension (LDQ, N)
[in]	LDQ	LDQ is INTEGER
[in,out]	Z	Z is COMPLEX*16 array, dimension (LDZ, N)
[in]	LDZ	LDZ is INTEGER
[in,out]	QC	QC is COMPLEX*16 array, dimension (LDQC, NBLOCK_DESIRED)
[in]	LDQC	LDQC is INTEGER
[in,out]	ZC	ZC is COMPLEX*16 array, dimension (LDZC, NBLOCK_DESIRED)
[in]	LDZC	LDZ is INTEGER
[out]	WORK	WORK is COMPLEX*16 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 204 of file zlaqz3.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
 
      COMPLEX*16, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ,
     $   * ), Z( LDZ, * ), QC( LDQC, * ), ZC( LDZC, * ), WORK( * ),
     $   ALPHA( * ), BETA( * )
 
      INTEGER, INTENT( OUT ) :: INFO
 
*     Parameters
      COMPLEX*16         CZERO, CONE
      parameter( czero = ( 0.0d+0, 0.0d+0 ), cone = ( 1.0d+0,
     $                     0.0d+0 ) )
      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 :: SAFMIN, SAFMAX, C, SCALE
      COMPLEX*16 :: TEMP, TEMP2, TEMP3, S
 
*     External Functions
      EXTERNAL :: xerbla, dlabad, zlaset, zlartg, zrot, zlaqz1, zgemm,
     $            zlacpy
      DOUBLE PRECISION, EXTERNAL :: DLAMCH
 
      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( 'ZLAQZ3', -info )
         RETURN
      END IF
 
*
*     Executable statements
*
 
*     Get machine constants
      safmin = dlamch( 'SAFE MINIMUM' )
      safmax = one/safmin
      CALL dlabad( safmin, safmax )
 
      IF ( ilo .GE. ihi ) THEN
         RETURN
      END IF
 
      IF ( ilschur ) THEN
         istartm = 1
         istopm = n
      ELSE
         istartm = ilo
         istopm = ihi
      END IF
 
      ns = nshifts
      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 zlaset( 'FULL', ns+1, ns+1, czero, cone, qc, ldqc )
      CALL zlaset( 'FULL', ns, ns, czero, cone, zc, ldzc )
 
      DO i = 1, ns
*        Introduce the shift
         scale = sqrt( abs( alpha( i ) ) ) * sqrt( abs( beta( i ) ) )
         IF( scale .GE. safmin .AND. scale .LE. safmax ) THEN
            alpha( i ) = alpha( i )/scale
            beta( i ) = beta( i )/scale
         END IF
 
         temp2 = beta( i )*a( ilo, ilo )-alpha( i )*b( ilo, ilo )
         temp3 = beta( i )*a( ilo+1, ilo )
 
         IF ( abs( temp2 ) .GT. safmax .OR.
     $      abs( temp3 ) .GT. safmax ) THEN
            temp2 = cone
            temp3 = czero
         END IF
 
         CALL zlartg( temp2, temp3, c, s, temp )
         CALL zrot( ns, a( ilo, ilo ), lda, a( ilo+1, ilo ), lda, c,
     $              s )
         CALL zrot( ns, b( ilo, ilo ), ldb, b( ilo+1, ilo ), ldb, c,
     $              s )
         CALL zrot( ns+1, qc( 1, 1 ), 1, qc( 1, 2 ), 1, c,
     $              dconjg( s ) )
        
*        Chase the shift down
         DO j = 1, ns-i
 
            CALL zlaqz1( .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 zgemm( 'C', 'N', sheight, swidth, sheight, cone, qc, ldqc,
     $               a( ilo, ilo+ns ), lda, czero, work, sheight )
         CALL zlacpy( 'ALL', sheight, swidth, work, sheight, a( ilo,
     $                ilo+ns ), lda )
         CALL zgemm( 'C', 'N', sheight, swidth, sheight, cone, qc, ldqc,
     $               b( ilo, ilo+ns ), ldb, czero, work, sheight )
         CALL zlacpy( 'ALL', sheight, swidth, work, sheight, b( ilo,
     $                ilo+ns ), ldb )
      END IF
      IF ( ilq ) THEN
         CALL zgemm( 'N', 'N', n, sheight, sheight, cone, q( 1, ilo ),
     $               ldq, qc, ldqc, czero, work, n )
         CALL zlacpy( '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 zgemm( 'N', 'N', sheight, swidth, swidth, cone,
     $               a( istartm, ilo ), lda, zc, ldzc, czero, work,
     $               sheight )
         CALL zlacpy( 'ALL', sheight, swidth, work, sheight, a( istartm,
     $                ilo ), lda )
         CALL zgemm( 'N', 'N', sheight, swidth, swidth, cone,
     $               b( istartm, ilo ), ldb, zc, ldzc, czero, work,
     $               sheight )
         CALL zlacpy( 'ALL', sheight, swidth, work, sheight, b( istartm,
     $                ilo ), ldb )
      END IF
      IF ( ilz ) THEN
         CALL zgemm( 'N', 'N', n, swidth, swidth, cone, z( 1, ilo ),
     $               ldz, zc, ldzc, czero, work, n )
         CALL zlacpy( '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 zlaset( 'FULL', ns+np, ns+np, czero, cone, qc, ldqc )
         CALL zlaset( 'FULL', ns+np, ns+np, czero, cone, zc, ldzc )
 
*        Near the diagonal shift chase
         DO i = ns-1, 0, -1
            DO j = 0, np-1
*              Move down the block with index k+i+j, updating
*              the (ns+np x ns+np) block:
*              (k:k+ns+np,k:k+ns+np-1)
               CALL zlaqz1( .true., .true., k+i+j, 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 zgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
     $                  ldqc, a( k+1, k+ns+np ), lda, czero, work,
     $                  sheight )
            CALL zlacpy( 'ALL', sheight, swidth, work, sheight, a( k+1,
     $                   k+ns+np ), lda )
            CALL zgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
     $                  ldqc, b( k+1, k+ns+np ), ldb, czero, work,
     $                  sheight )
            CALL zlacpy( 'ALL', sheight, swidth, work, sheight, b( k+1,
     $                   k+ns+np ), ldb )
         END IF
         IF ( ilq ) THEN
            CALL zgemm( 'N', 'N', n, nblock, nblock, cone, q( 1, k+1 ),
     $                  ldq, qc, ldqc, czero, work, n )
            CALL zlacpy( '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 zgemm( 'N', 'N', sheight, swidth, swidth, cone,
     $                  a( istartm, k ), lda, zc, ldzc, czero, work,
     $                  sheight )
            CALL zlacpy( 'ALL', sheight, swidth, work, sheight,
     $                   a( istartm, k ), lda )
            CALL zgemm( 'N', 'N', sheight, swidth, swidth, cone,
     $                  b( istartm, k ), ldb, zc, ldzc, czero, work,
     $                  sheight )
            CALL zlacpy( 'ALL', sheight, swidth, work, sheight,
     $                   b( istartm, k ), ldb )
         END IF
         IF ( ilz ) THEN
            CALL zgemm( 'N', 'N', n, nblock, nblock, cone, z( 1, k ),
     $                  ldz, zc, ldzc, czero, work, n )
            CALL zlacpy( '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 zlaset( 'FULL', ns, ns, czero, cone, qc, ldqc )
      CALL zlaset( 'FULL', ns+1, ns+1, czero, cone, 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
*        Chase the shift down to the bottom right corner
         DO ishift = ihi-i, ihi-1
            CALL zlaqz1( .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 zgemm( 'C', 'N', sheight, swidth, sheight, cone, qc, ldqc,
     $               a( ihi-ns+1, ihi+1 ), lda, czero, work, sheight )
         CALL zlacpy( 'ALL', sheight, swidth, work, sheight,
     $                a( ihi-ns+1, ihi+1 ), lda )
         CALL zgemm( 'C', 'N', sheight, swidth, sheight, cone, qc, ldqc,
     $               b( ihi-ns+1, ihi+1 ), ldb, czero, work, sheight )
         CALL zlacpy( 'ALL', sheight, swidth, work, sheight,
     $                b( ihi-ns+1, ihi+1 ), ldb )
      END IF
      IF ( ilq ) THEN
         CALL zgemm( 'N', 'N', n, ns, ns, cone, q( 1, ihi-ns+1 ), ldq,
     $               qc, ldqc, czero, work, n )
         CALL zlacpy( '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 zgemm( 'N', 'N', sheight, swidth, swidth, cone,
     $               a( istartm, ihi-ns ), lda, zc, ldzc, czero, work,
     $               sheight )
         CALL zlacpy( 'ALL', sheight, swidth, work, sheight, a( istartm,
     $                ihi-ns ), lda )
         CALL zgemm( 'N', 'N', sheight, swidth, swidth, cone,
     $               b( istartm, ihi-ns ), ldb, zc, ldzc, czero, work,
     $               sheight )
         CALL zlacpy( 'ALL', sheight, swidth, work, sheight, b( istartm,
     $                ihi-ns ), ldb )
      END IF
      IF ( ilz ) THEN
         CALL zgemm( 'N', 'N', n, ns+1, ns+1, cone, z( 1, ihi-ns ), ldz,
     $               zc, ldzc, czero, work, n )
         CALL zlacpy( '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: