LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

## ◆ claqz2()

 recursive subroutine claqz2 ( 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) NW, 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, integer, intent(out) NS, integer, intent(out) ND, complex, dimension( * ), intent(inout) ALPHA, complex, dimension( * ), intent(inout) BETA, complex, dimension( ldqc, * ) QC, integer, intent(in) LDQC, complex, dimension( ldzc, * ) ZC, integer, intent(in) LDZC, complex, dimension( * ) WORK, integer, intent(in) LWORK, real, dimension( * ) RWORK, integer, intent(in) REC, integer, intent(out) INFO )

CLAQZ2

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

Purpose:
` CLAQZ2 performs AED`
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 ILO and IHI mark the rows and columns of (A,B) which are to be normalized``` [in] NW ``` NW is INTEGER The desired size of the deflation window.``` [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` [out] NS ``` NS is INTEGER The number of unconverged eigenvalues available to use as shifts.``` [out] ND ``` ND is INTEGER The number of converged eigenvalues found.``` [out] ALPHA ``` ALPHA is COMPLEX array, dimension (N) Each scalar alpha defining an eigenvalue of GNEP.``` [out] BETA ``` BETA is COMPLEX array, dimension (N) The scalars beta that define the eigenvalues of GNEP. Together, the quantities alpha = ALPHA(j) and beta = BETA(j) represent the j-th eigenvalue of the matrix pair (A,B), in one of the forms lambda = alpha/beta or mu = beta/alpha. Since either lambda or mu may overflow, they should not, in general, be computed.``` [in,out] QC ` QC is COMPLEX array, dimension (LDQC, NW)` [in] LDQC ` LDQC is INTEGER` [in,out] ZC ` ZC is COMPLEX array, dimension (LDZC, NW)` [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] RWORK ` RWORK is REAL array, dimension (N)` [in] REC ``` REC is INTEGER REC indicates the current recursion level. Should be set to 0 on first call. \param[out] INFO \verbatim INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value```
Date
May 2020

Definition at line 230 of file claqz2.f.

234  IMPLICIT NONE
235
236 * Arguments
237  LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
238  INTEGER, INTENT( IN ) :: N, ILO, IHI, NW, LDA, LDB, LDQ, LDZ,
239  \$ LDQC, LDZC, LWORK, REC
240
241  COMPLEX, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
242  \$ Z( LDZ, * ), ALPHA( * ), BETA( * )
243  INTEGER, INTENT( OUT ) :: NS, ND, INFO
244  COMPLEX :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
245  REAL :: RWORK( * )
246
247 * Parameters
248  COMPLEX CZERO, CONE
249  parameter( czero = ( 0.0, 0.0 ), cone = ( 1.0, 0.0 ) )
250  REAL :: ZERO, ONE, HALF
251  parameter( zero = 0.0, one = 1.0, half = 0.5 )
252
253 * Local Scalars
254  INTEGER :: JW, KWTOP, KWBOT, ISTOPM, ISTARTM, K, K2, CTGEXC_INFO,
255  \$ IFST, ILST, LWORKREQ, QZ_SMALL_INFO
256  REAL :: SMLNUM, ULP, SAFMIN, SAFMAX, C1, TEMPR
257  COMPLEX :: S, S1, TEMP
258
259 * External Functions
260  EXTERNAL :: xerbla, claqz0, claqz1, slabad, clacpy, claset, cgemm,
261  \$ ctgexc, clartg, crot
262  REAL, EXTERNAL :: SLAMCH
263
264  info = 0
265
266 * Set up deflation window
267  jw = min( nw, ihi-ilo+1 )
268  kwtop = ihi-jw+1
269  IF ( kwtop .EQ. ilo ) THEN
270  s = czero
271  ELSE
272  s = a( kwtop, kwtop-1 )
273  END IF
274
275 * Determine required workspace
276  ifst = 1
277  ilst = jw
278  CALL claqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
279  \$ b( kwtop, kwtop ), ldb, alpha, beta, qc, ldqc, zc,
280  \$ ldzc, work, -1, rwork, rec+1, qz_small_info )
281  lworkreq = int( work( 1 ) )+2*jw**2
282  lworkreq = max( lworkreq, n*nw, 2*nw**2+n )
283  IF ( lwork .EQ.-1 ) THEN
284 * workspace query, quick return
285  work( 1 ) = lworkreq
286  RETURN
287  ELSE IF ( lwork .LT. lworkreq ) THEN
288  info = -26
289  END IF
290
291  IF( info.NE.0 ) THEN
292  CALL xerbla( 'CLAQZ2', -info )
293  RETURN
294  END IF
295
296 * Get machine constants
297  safmin = slamch( 'SAFE MINIMUM' )
298  safmax = one/safmin
299  CALL slabad( safmin, safmax )
300  ulp = slamch( 'PRECISION' )
301  smlnum = safmin*( real( n )/ulp )
302
303  IF ( ihi .EQ. kwtop ) THEN
304 * 1 by 1 deflation window, just try a regular deflation
305  alpha( kwtop ) = a( kwtop, kwtop )
306  beta( kwtop ) = b( kwtop, kwtop )
307  ns = 1
308  nd = 0
309  IF ( abs( s ) .LE. max( smlnum, ulp*abs( a( kwtop,
310  \$ kwtop ) ) ) ) THEN
311  ns = 0
312  nd = 1
313  IF ( kwtop .GT. ilo ) THEN
314  a( kwtop, kwtop-1 ) = czero
315  END IF
316  END IF
317  END IF
318
319
320 * Store window in case of convergence failure
321  CALL clacpy( 'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
322  CALL clacpy( 'ALL', jw, jw, b( kwtop, kwtop ), ldb, work( jw**2+
323  \$ 1 ), jw )
324
325 * Transform window to real schur form
326  CALL claset( 'FULL', jw, jw, czero, cone, qc, ldqc )
327  CALL claset( 'FULL', jw, jw, czero, cone, zc, ldzc )
328  CALL claqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
329  \$ b( kwtop, kwtop ), ldb, alpha, beta, qc, ldqc, zc,
330  \$ ldzc, work( 2*jw**2+1 ), lwork-2*jw**2, rwork,
331  \$ rec+1, qz_small_info )
332
333  IF( qz_small_info .NE. 0 ) THEN
334 * Convergence failure, restore the window and exit
335  nd = 0
336  ns = jw-qz_small_info
337  CALL clacpy( 'ALL', jw, jw, work, jw, a( kwtop, kwtop ), lda )
338  CALL clacpy( 'ALL', jw, jw, work( jw**2+1 ), jw, b( kwtop,
339  \$ kwtop ), ldb )
340  RETURN
341  END IF
342
343 * Deflation detection loop
344  IF ( kwtop .EQ. ilo .OR. s .EQ. czero ) THEN
345  kwbot = kwtop-1
346  ELSE
347  kwbot = ihi
348  k = 1
349  k2 = 1
350  DO WHILE ( k .LE. jw )
351 * Try to deflate eigenvalue
352  tempr = abs( a( kwbot, kwbot ) )
353  IF( tempr .EQ. zero ) THEN
354  tempr = abs( s )
355  END IF
356  IF ( ( abs( s*qc( 1, kwbot-kwtop+1 ) ) ) .LE. max( ulp*
357  \$ tempr, smlnum ) ) THEN
358 * Deflatable
359  kwbot = kwbot-1
360  ELSE
361 * Not deflatable, move out of the way
362  ifst = kwbot-kwtop+1
363  ilst = k2
364  CALL ctgexc( .true., .true., jw, a( kwtop, kwtop ),
365  \$ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
366  \$ zc, ldzc, ifst, ilst, ctgexc_info )
367  k2 = k2+1
368  END IF
369
370  k = k+1
371  END DO
372  END IF
373
374 * Store eigenvalues
375  nd = ihi-kwbot
376  ns = jw-nd
377  k = kwtop
378  DO WHILE ( k .LE. ihi )
379  alpha( k ) = a( k, k )
380  beta( k ) = b( k, k )
381  k = k+1
382  END DO
383
384  IF ( kwtop .NE. ilo .AND. s .NE. czero ) THEN
385 * Reflect spike back, this will create optimally packed bulges
386  a( kwtop:kwbot, kwtop-1 ) = a( kwtop, kwtop-1 ) *conjg( qc( 1,
387  \$ 1:jw-nd ) )
388  DO k = kwbot-1, kwtop, -1
389  CALL clartg( a( k, kwtop-1 ), a( k+1, kwtop-1 ), c1, s1,
390  \$ temp )
391  a( k, kwtop-1 ) = temp
392  a( k+1, kwtop-1 ) = czero
393  k2 = max( kwtop, k-1 )
394  CALL crot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda, c1,
395  \$ s1 )
396  CALL crot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1, k-1 ),
397  \$ ldb, c1, s1 )
398  CALL crot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1, k+1-kwtop+1 ),
399  \$ 1, c1, conjg( s1 ) )
400  END DO
401
402 * Chase bulges down
403  istartm = kwtop
404  istopm = ihi
405  k = kwbot-1
406  DO WHILE ( k .GE. kwtop )
407
408 * Move bulge down and remove it
409  DO k2 = k, kwbot-1
410  CALL claqz1( .true., .true., k2, kwtop, kwtop+jw-1,
411  \$ kwbot, a, lda, b, ldb, jw, kwtop, qc, ldqc,
412  \$ jw, kwtop, zc, ldzc )
413  END DO
414
415  k = k-1
416  END DO
417
418  END IF
419
420 * Apply Qc and Zc to rest of the matrix
421  IF ( ilschur ) THEN
422  istartm = 1
423  istopm = n
424  ELSE
425  istartm = ilo
426  istopm = ihi
427  END IF
428
429  IF ( istopm-ihi > 0 ) THEN
430  CALL cgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
431  \$ a( kwtop, ihi+1 ), lda, czero, work, jw )
432  CALL clacpy( 'ALL', jw, istopm-ihi, work, jw, a( kwtop,
433  \$ ihi+1 ), lda )
434  CALL cgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
435  \$ b( kwtop, ihi+1 ), ldb, czero, work, jw )
436  CALL clacpy( 'ALL', jw, istopm-ihi, work, jw, b( kwtop,
437  \$ ihi+1 ), ldb )
438  END IF
439  IF ( ilq ) THEN
440  CALL cgemm( 'N', 'N', n, jw, jw, cone, q( 1, kwtop ), ldq, qc,
441  \$ ldqc, czero, work, n )
442  CALL clacpy( 'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
443  END IF
444
445  IF ( kwtop-1-istartm+1 > 0 ) THEN
446  CALL cgemm( 'N', 'N', kwtop-istartm, jw, jw, cone, a( istartm,
447  \$ kwtop ), lda, zc, ldzc, czero, work,
448  \$ kwtop-istartm )
449  CALL clacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
450  \$ a( istartm, kwtop ), lda )
451  CALL cgemm( 'N', 'N', kwtop-istartm, jw, jw, cone, b( istartm,
452  \$ kwtop ), ldb, zc, ldzc, czero, work,
453  \$ kwtop-istartm )
454  CALL clacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
455  \$ b( istartm, kwtop ), ldb )
456  END IF
457  IF ( ilz ) THEN
458  CALL cgemm( 'N', 'N', n, jw, jw, cone, z( 1, kwtop ), ldz, zc,
459  \$ ldzc, czero, work, n )
460  CALL clacpy( 'ALL', n, jw, work, n, z( 1, kwtop ), ldz )
461  END IF
462
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
subroutine clartg(f, g, c, s, r)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition: clartg.f90:118
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
subroutine claqz1(ILQ, ILZ, K, ISTARTM, ISTOPM, IHI, A, LDA, B, LDB, NQ, QSTART, Q, LDQ, NZ, ZSTART, Z, LDZ)
CLAQZ1
Definition: claqz1.f:173
subroutine ctgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, INFO)
CTGEXC
Definition: ctgexc.f:200
recursive subroutine claqz0(WANTS, WANTQ, WANTZ, N, ILO, IHI, A, LDA, B, LDB, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, REC, INFO)
CLAQZ0
Definition: claqz0.f:284
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: claset.f:106
subroutine crot(N, CX, INCX, CY, INCY, C, S)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition: crot.f:103
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: