claqz3()

subroutine claqz3	(	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, dimension( * ), intent(inout)	alpha,
		complex, dimension( * ), intent(inout)	beta,
		complex, dimension( lda, * ), intent(inout)	a,
		integer, intent(in)	lda,
		complex, dimension( ldb, * ), intent(inout)	b,
		integer, intent(in)	ldb,
		complex, dimension( ldq, * ), intent(inout)	q,
		integer, intent(in)	ldq,
		complex, dimension( ldz, * ), intent(inout)	z,
		integer, intent(in)	ldz,
		complex, dimension( ldqc, * ), intent(inout)	qc,
		integer, intent(in)	ldqc,
		complex, dimension( ldzc, * ), intent(inout)	zc,
		integer, intent(in)	ldzc,
		complex, dimension( * ), intent(inout)	work,
		integer, intent(in)	lwork,
		integer, intent(out)	info )

CLAQZ3

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

Purpose:

!>
!> CLAQZ3 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 array. SR contains !> the alpha parts of the shifts to use. !>
[in]	BETA	!> BETA is COMPLEX array. SS contains !> the scale of the shifts to use. !>
[in,out]	A	!> A is COMPLEX 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 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 array, dimension (LDQ, N) !>
[in]	LDQ	!> LDQ is INTEGER !>
[in,out]	Z	!> Z is COMPLEX array, dimension (LDZ, N) !>
[in]	LDZ	!> LDZ is INTEGER !>
[in,out]	QC	!> QC is COMPLEX array, dimension (LDQC, NBLOCK_DESIRED) !>
[in]	LDQC	!> LDQC is INTEGER !>
[in,out]	ZC	!> ZC is COMPLEX array, dimension (LDZC, NBLOCK_DESIRED) !>
[in]	LDZC	!> LDZ is INTEGER !>
[out]	WORK	!> WORK is COMPLEX 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 201 of file claqz3.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, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
     $   Z( LDZ, * ), QC( LDQC, * ), ZC( LDZC, * ), WORK( * ),
     $   ALPHA( * ), BETA( * )
 
      INTEGER, INTENT( OUT ) :: INFO
 
*     Parameters
      COMPLEX         CZERO, CONE
      parameter( czero = ( 0.0, 0.0 ), cone = ( 1.0, 0.0 ) )
      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 :: SAFMIN, SAFMAX, C, SCALE
      COMPLEX :: TEMP, TEMP2, TEMP3, S
 
*     External Functions
      EXTERNAL :: xerbla, claset, clartg, crot,
     $            claqz1, cgemm, clacpy
      REAL, EXTERNAL :: SLAMCH
 
      info = 0
      IF ( nblock_desired .LT. nshifts+1 ) THEN
         info = -8
      END IF
      IF ( lwork .EQ.-1 ) THEN
*        workspace query, quick return
         work( 1 ) = cmplx( n*nblock_desired )
         RETURN
      ELSE IF ( lwork .LT. n*nblock_desired ) THEN
         info = -25
      END IF
 
      IF( info.NE.0 ) THEN
         CALL xerbla( 'CLAQZ3', -info )
         RETURN
      END IF
 
*
*     Executable statements
*
 
*     Get machine constants
      safmin = slamch( 'SAFE MINIMUM' )
      safmax = one/safmin
 
      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 claset( 'FULL', ns+1, ns+1, czero, cone, qc, ldqc )
      CALL claset( '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 clartg( temp2, temp3, c, s, temp )
         CALL crot( ns, a( ilo, ilo ), lda, a( ilo+1, ilo ), lda,
     $              c, s )
         CALL crot( ns, b( ilo, ilo ), ldb, b( ilo+1, ilo ), ldb,
     $              c, s )
         CALL crot( ns+1, qc( 1, 1 ), 1, qc( 1, 2 ), 1, c,
     $              conjg( s ) )
        
*        Chase the shift down
         DO j = 1, ns-i
 
            CALL claqz1( .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 cgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
     $               ldqc,
     $               a( ilo, ilo+ns ), lda, czero, work, sheight )
         CALL clacpy( 'ALL', sheight, swidth, work, sheight, a( ilo,
     $                ilo+ns ), lda )
         CALL cgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
     $               ldqc,
     $               b( ilo, ilo+ns ), ldb, czero, work, sheight )
         CALL clacpy( 'ALL', sheight, swidth, work, sheight, b( ilo,
     $                ilo+ns ), ldb )
      END IF
      IF ( ilq ) THEN
        CALL cgemm( 'N', 'N', n, sheight, sheight, cone, q( 1, ilo ),
     $              ldq, qc, ldqc, czero, work, n )
         CALL clacpy( '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 cgemm( 'N', 'N', sheight, swidth, swidth, cone,
     $               a( istartm, ilo ), lda, zc, ldzc, czero, work,
     $               sheight )
         CALL clacpy( 'ALL', sheight, swidth, work, sheight,
     $                a( istartm,
     $                ilo ), lda )
         CALL cgemm( 'N', 'N', sheight, swidth, swidth, cone,
     $               b( istartm, ilo ), ldb, zc, ldzc, czero, work,
     $               sheight )
         CALL clacpy( 'ALL', sheight, swidth, work, sheight,
     $                b( istartm,
     $                ilo ), ldb )
      END IF
      IF ( ilz ) THEN
         CALL cgemm( 'N', 'N', n, swidth, swidth, cone, z( 1, ilo ),
     $               ldz, zc, ldzc, czero, work, n )
         CALL clacpy( '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 claset( 'FULL', ns+np, ns+np, czero, cone, qc, ldqc )
         CALL claset( '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 claqz1( .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 cgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
     $                  ldqc, a( k+1, k+ns+np ), lda, czero, work,
     $                  sheight )
            CALL clacpy( 'ALL', sheight, swidth, work, sheight,
     $                   a( k+1,
     $                   k+ns+np ), lda )
            CALL cgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
     $                  ldqc, b( k+1, k+ns+np ), ldb, czero, work,
     $                  sheight )
            CALL clacpy( 'ALL', sheight, swidth, work, sheight,
     $                   b( k+1,
     $                   k+ns+np ), ldb )
         END IF
         IF ( ilq ) THEN
            CALL cgemm( 'N', 'N', n, nblock, nblock, cone, q( 1,
     $                  k+1 ),
     $                  ldq, qc, ldqc, czero, work, n )
            CALL clacpy( '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 cgemm( 'N', 'N', sheight, swidth, swidth, cone,
     $                  a( istartm, k ), lda, zc, ldzc, czero, work,
     $                  sheight )
            CALL clacpy( 'ALL', sheight, swidth, work, sheight,
     $                   a( istartm, k ), lda )
            CALL cgemm( 'N', 'N', sheight, swidth, swidth, cone,
     $                  b( istartm, k ), ldb, zc, ldzc, czero, work,
     $                  sheight )
            CALL clacpy( 'ALL', sheight, swidth, work, sheight,
     $                   b( istartm, k ), ldb )
         END IF
         IF ( ilz ) THEN
            CALL cgemm( 'N', 'N', n, nblock, nblock, cone, z( 1, k ),
     $                  ldz, zc, ldzc, czero, work, n )
            CALL clacpy( '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 claset( 'FULL', ns, ns, czero, cone, qc, ldqc )
      CALL claset( '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 claqz1( .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 cgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
     $               ldqc,
     $               a( ihi-ns+1, ihi+1 ), lda, czero, work, sheight )
         CALL clacpy( 'ALL', sheight, swidth, work, sheight,
     $                a( ihi-ns+1, ihi+1 ), lda )
         CALL cgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
     $               ldqc,
     $               b( ihi-ns+1, ihi+1 ), ldb, czero, work, sheight )
         CALL clacpy( 'ALL', sheight, swidth, work, sheight,
     $                b( ihi-ns+1, ihi+1 ), ldb )
      END IF
      IF ( ilq ) THEN
         CALL cgemm( 'N', 'N', n, ns, ns, cone, q( 1, ihi-ns+1 ),
     $               ldq,
     $               qc, ldqc, czero, work, n )
         CALL clacpy( '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 cgemm( 'N', 'N', sheight, swidth, swidth, cone,
     $               a( istartm, ihi-ns ), lda, zc, ldzc, czero, work,
     $               sheight )
         CALL clacpy( 'ALL', sheight, swidth, work, sheight,
     $                a( istartm,
     $                ihi-ns ), lda )
         CALL cgemm( 'N', 'N', sheight, swidth, swidth, cone,
     $               b( istartm, ihi-ns ), ldb, zc, ldzc, czero, work,
     $               sheight )
         CALL clacpy( 'ALL', sheight, swidth, work, sheight,
     $                b( istartm,
     $                ihi-ns ), ldb )
      END IF
      IF ( ilz ) THEN
         CALL cgemm( 'N', 'N', n, ns+1, ns+1, cone, z( 1, ihi-ns ),
     $               ldz,
     $               zc, ldzc, czero, work, n )
         CALL clacpy( '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: