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 203 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 ) = 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: