LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

## ◆ dlaqz3()

 recursive subroutine dlaqz3 ( 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, 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, integer, intent(out) NS, integer, intent(out) ND, double precision, dimension( * ), intent(inout) ALPHAR, double precision, dimension( * ), intent(inout) ALPHAI, double precision, dimension( * ), intent(inout) BETA, double precision, dimension( ldqc, * ) QC, integer, intent(in) LDQC, double precision, dimension( ldzc, * ) ZC, integer, intent(in) LDZC, double precision, dimension( * ) WORK, integer, intent(in) LWORK, integer, intent(in) REC, integer, intent(out) INFO )

DLAQZ3

Purpose:
` DLAQZ3 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 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` [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] ALPHAR ``` ALPHAR is DOUBLE PRECISION array, dimension (N) The real parts of each scalar alpha defining an eigenvalue of GNEP.``` [out] ALPHAI ``` ALPHAI is DOUBLE PRECISION array, dimension (N) The imaginary parts of each scalar alpha defining an eigenvalue of GNEP. If ALPHAI(j) is zero, then the j-th eigenvalue is real; if positive, then the j-th and (j+1)-st eigenvalues are a complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).``` [out] BETA ``` BETA is DOUBLE PRECISION array, dimension (N) The scalars beta that define the eigenvalues of GNEP. Together, the quantities alpha = (ALPHAR(j),ALPHAI(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 DOUBLE PRECISION array, dimension (LDQC, NW)` [in] LDQC ` LDQC is INTEGER` [in,out] ZC ` ZC is DOUBLE PRECISION array, dimension (LDZC, NW)` [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.``` [in] REC ``` REC is INTEGER REC indicates the current recursion level. Should be set to 0 on first call.``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value```
Date
May 2020

Definition at line 235 of file dlaqz3.f.

239  IMPLICIT NONE
240
241 * Arguments
242  LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
243  INTEGER, INTENT( IN ) :: N, ILO, IHI, NW, LDA, LDB, LDQ, LDZ,
244  \$ LDQC, LDZC, LWORK, REC
245
246  DOUBLE PRECISION, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ),
247  \$ Q( LDQ, * ), Z( LDZ, * ), ALPHAR( * ),
248  \$ ALPHAI( * ), BETA( * )
249  INTEGER, INTENT( OUT ) :: NS, ND, INFO
250  DOUBLE PRECISION :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
251
252 * Parameters
253  DOUBLE PRECISION :: ZERO, ONE, HALF
254  parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
255
256 * Local Scalars
257  LOGICAL :: BULGE
258  INTEGER :: JW, KWTOP, KWBOT, ISTOPM, ISTARTM, K, K2, DTGEXC_INFO,
259  \$ IFST, ILST, LWORKREQ, QZ_SMALL_INFO
260  DOUBLE PRECISION :: S, SMLNUM, ULP, SAFMIN, SAFMAX, C1, S1, TEMP
261
262 * External Functions
263  EXTERNAL :: xerbla, dtgexc, dlabad, dlaqz0, dlacpy, dlaset,
265  DOUBLE PRECISION, EXTERNAL :: DLAMCH
266
267  info = 0
268
269 * Set up deflation window
270  jw = min( nw, ihi-ilo+1 )
271  kwtop = ihi-jw+1
272  IF ( kwtop .EQ. ilo ) THEN
273  s = zero
274  ELSE
275  s = a( kwtop, kwtop-1 )
276  END IF
277
278 * Determine required workspace
279  ifst = 1
280  ilst = jw
281  CALL dtgexc( .true., .true., jw, a, lda, b, ldb, qc, ldqc, zc,
282  \$ ldzc, ifst, ilst, work, -1, dtgexc_info )
283  lworkreq = int( work( 1 ) )
284  CALL dlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
285  \$ b( kwtop, kwtop ), ldb, alphar, alphai, beta, qc,
286  \$ ldqc, zc, ldzc, work, -1, rec+1, qz_small_info )
287  lworkreq = max( lworkreq, int( work( 1 ) )+2*jw**2 )
288  lworkreq = max( lworkreq, n*nw, 2*nw**2+n )
289  IF ( lwork .EQ.-1 ) THEN
290 * workspace query, quick return
291  work( 1 ) = lworkreq
292  RETURN
293  ELSE IF ( lwork .LT. lworkreq ) THEN
294  info = -26
295  END IF
296
297  IF( info.NE.0 ) THEN
298  CALL xerbla( 'DLAQZ3', -info )
299  RETURN
300  END IF
301
302 * Get machine constants
303  safmin = dlamch( 'SAFE MINIMUM' )
304  safmax = one/safmin
305  CALL dlabad( safmin, safmax )
306  ulp = dlamch( 'PRECISION' )
307  smlnum = safmin*( dble( n )/ulp )
308
309  IF ( ihi .EQ. kwtop ) THEN
310 * 1 by 1 deflation window, just try a regular deflation
311  alphar( kwtop ) = a( kwtop, kwtop )
312  alphai( kwtop ) = zero
313  beta( kwtop ) = b( kwtop, kwtop )
314  ns = 1
315  nd = 0
316  IF ( abs( s ) .LE. max( smlnum, ulp*abs( a( kwtop,
317  \$ kwtop ) ) ) ) THEN
318  ns = 0
319  nd = 1
320  IF ( kwtop .GT. ilo ) THEN
321  a( kwtop, kwtop-1 ) = zero
322  END IF
323  END IF
324  END IF
325
326
327 * Store window in case of convergence failure
328  CALL dlacpy( 'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
329  CALL dlacpy( 'ALL', jw, jw, b( kwtop, kwtop ), ldb, work( jw**2+
330  \$ 1 ), jw )
331
332 * Transform window to real schur form
333  CALL dlaset( 'FULL', jw, jw, zero, one, qc, ldqc )
334  CALL dlaset( 'FULL', jw, jw, zero, one, zc, ldzc )
335  CALL dlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
336  \$ b( kwtop, kwtop ), ldb, alphar, alphai, beta, qc,
337  \$ ldqc, zc, ldzc, work( 2*jw**2+1 ), lwork-2*jw**2,
338  \$ rec+1, qz_small_info )
339
340  IF( qz_small_info .NE. 0 ) THEN
341 * Convergence failure, restore the window and exit
342  nd = 0
343  ns = jw-qz_small_info
344  CALL dlacpy( 'ALL', jw, jw, work, jw, a( kwtop, kwtop ), lda )
345  CALL dlacpy( 'ALL', jw, jw, work( jw**2+1 ), jw, b( kwtop,
346  \$ kwtop ), ldb )
347  RETURN
348  END IF
349
350 * Deflation detection loop
351  IF ( kwtop .EQ. ilo .OR. s .EQ. zero ) THEN
352  kwbot = kwtop-1
353  ELSE
354  kwbot = ihi
355  k = 1
356  k2 = 1
357  DO WHILE ( k .LE. jw )
358  bulge = .false.
359  IF ( kwbot-kwtop+1 .GE. 2 ) THEN
360  bulge = a( kwbot, kwbot-1 ) .NE. zero
361  END IF
362  IF ( bulge ) THEN
363
364 * Try to deflate complex conjugate eigenvalue pair
365  temp = abs( a( kwbot, kwbot ) )+sqrt( abs( a( kwbot,
366  \$ kwbot-1 ) ) )*sqrt( abs( a( kwbot-1, kwbot ) ) )
367  IF( temp .EQ. zero )THEN
368  temp = abs( s )
369  END IF
370  IF ( max( abs( s*qc( 1, kwbot-kwtop ) ), abs( s*qc( 1,
371  \$ kwbot-kwtop+1 ) ) ) .LE. max( smlnum,
372  \$ ulp*temp ) ) THEN
373 * Deflatable
374  kwbot = kwbot-2
375  ELSE
376 * Not deflatable, move out of the way
377  ifst = kwbot-kwtop+1
378  ilst = k2
379  CALL dtgexc( .true., .true., jw, a( kwtop, kwtop ),
380  \$ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
381  \$ zc, ldzc, ifst, ilst, work, lwork,
382  \$ dtgexc_info )
383  k2 = k2+2
384  END IF
385  k = k+2
386  ELSE
387
388 * Try to deflate real eigenvalue
389  temp = abs( a( kwbot, kwbot ) )
390  IF( temp .EQ. zero ) THEN
391  temp = abs( s )
392  END IF
393  IF ( ( abs( s*qc( 1, kwbot-kwtop+1 ) ) ) .LE. max( ulp*
394  \$ temp, smlnum ) ) THEN
395 * Deflatable
396  kwbot = kwbot-1
397  ELSE
398 * Not deflatable, move out of the way
399  ifst = kwbot-kwtop+1
400  ilst = k2
401  CALL dtgexc( .true., .true., jw, a( kwtop, kwtop ),
402  \$ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
403  \$ zc, ldzc, ifst, ilst, work, lwork,
404  \$ dtgexc_info )
405  k2 = k2+1
406  END IF
407
408  k = k+1
409
410  END IF
411  END DO
412  END IF
413
414 * Store eigenvalues
415  nd = ihi-kwbot
416  ns = jw-nd
417  k = kwtop
418  DO WHILE ( k .LE. ihi )
419  bulge = .false.
420  IF ( k .LT. ihi ) THEN
421  IF ( a( k+1, k ) .NE. zero ) THEN
422  bulge = .true.
423  END IF
424  END IF
425  IF ( bulge ) THEN
426 * 2x2 eigenvalue block
427  CALL dlag2( a( k, k ), lda, b( k, k ), ldb, safmin,
428  \$ beta( k ), beta( k+1 ), alphar( k ),
429  \$ alphar( k+1 ), alphai( k ) )
430  alphai( k+1 ) = -alphai( k )
431  k = k+2
432  ELSE
433 * 1x1 eigenvalue block
434  alphar( k ) = a( k, k )
435  alphai( k ) = zero
436  beta( k ) = b( k, k )
437  k = k+1
438  END IF
439  END DO
440
441  IF ( kwtop .NE. ilo .AND. s .NE. zero ) THEN
442 * Reflect spike back, this will create optimally packed bulges
443  a( kwtop:kwbot, kwtop-1 ) = a( kwtop, kwtop-1 )*qc( 1,
444  \$ 1:jw-nd )
445  DO k = kwbot-1, kwtop, -1
446  CALL dlartg( a( k, kwtop-1 ), a( k+1, kwtop-1 ), c1, s1,
447  \$ temp )
448  a( k, kwtop-1 ) = temp
449  a( k+1, kwtop-1 ) = zero
450  k2 = max( kwtop, k-1 )
451  CALL drot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda, c1,
452  \$ s1 )
453  CALL drot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1, k-1 ),
454  \$ ldb, c1, s1 )
455  CALL drot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1, k+1-kwtop+1 ),
456  \$ 1, c1, s1 )
457  END DO
458
459 * Chase bulges down
460  istartm = kwtop
461  istopm = ihi
462  k = kwbot-1
463  DO WHILE ( k .GE. kwtop )
464  IF ( ( k .GE. kwtop+1 ) .AND. a( k+1, k-1 ) .NE. zero ) THEN
465
466 * Move double pole block down and remove it
467  DO k2 = k-1, kwbot-2
468  CALL dlaqz2( .true., .true., k2, kwtop, kwtop+jw-1,
469  \$ kwbot, a, lda, b, ldb, jw, kwtop, qc,
470  \$ ldqc, jw, kwtop, zc, ldzc )
471  END DO
472
473  k = k-2
474  ELSE
475
476 * k points to single shift
477  DO k2 = k, kwbot-2
478
479 * Move shift down
480  CALL dlartg( b( k2+1, k2+1 ), b( k2+1, k2 ), c1, s1,
481  \$ temp )
482  b( k2+1, k2+1 ) = temp
483  b( k2+1, k2 ) = zero
484  CALL drot( k2+2-istartm+1, a( istartm, k2+1 ), 1,
485  \$ a( istartm, k2 ), 1, c1, s1 )
486  CALL drot( k2-istartm+1, b( istartm, k2+1 ), 1,
487  \$ b( istartm, k2 ), 1, c1, s1 )
488  CALL drot( jw, zc( 1, k2+1-kwtop+1 ), 1, zc( 1,
489  \$ k2-kwtop+1 ), 1, c1, s1 )
490
491  CALL dlartg( a( k2+1, k2 ), a( k2+2, k2 ), c1, s1,
492  \$ temp )
493  a( k2+1, k2 ) = temp
494  a( k2+2, k2 ) = zero
495  CALL drot( istopm-k2, a( k2+1, k2+1 ), lda, a( k2+2,
496  \$ k2+1 ), lda, c1, s1 )
497  CALL drot( istopm-k2, b( k2+1, k2+1 ), ldb, b( k2+2,
498  \$ k2+1 ), ldb, c1, s1 )
499  CALL drot( jw, qc( 1, k2+1-kwtop+1 ), 1, qc( 1,
500  \$ k2+2-kwtop+1 ), 1, c1, s1 )
501
502  END DO
503
504 * Remove the shift
505  CALL dlartg( b( kwbot, kwbot ), b( kwbot, kwbot-1 ), c1,
506  \$ s1, temp )
507  b( kwbot, kwbot ) = temp
508  b( kwbot, kwbot-1 ) = zero
509  CALL drot( kwbot-istartm, b( istartm, kwbot ), 1,
510  \$ b( istartm, kwbot-1 ), 1, c1, s1 )
511  CALL drot( kwbot-istartm+1, a( istartm, kwbot ), 1,
512  \$ a( istartm, kwbot-1 ), 1, c1, s1 )
513  CALL drot( jw, zc( 1, kwbot-kwtop+1 ), 1, zc( 1,
514  \$ kwbot-1-kwtop+1 ), 1, c1, s1 )
515
516  k = k-1
517  END IF
518  END DO
519
520  END IF
521
522 * Apply Qc and Zc to rest of the matrix
523  IF ( ilschur ) THEN
524  istartm = 1
525  istopm = n
526  ELSE
527  istartm = ilo
528  istopm = ihi
529  END IF
530
531  IF ( istopm-ihi > 0 ) THEN
532  CALL dgemm( 'T', 'N', jw, istopm-ihi, jw, one, qc, ldqc,
533  \$ a( kwtop, ihi+1 ), lda, zero, work, jw )
534  CALL dlacpy( 'ALL', jw, istopm-ihi, work, jw, a( kwtop,
535  \$ ihi+1 ), lda )
536  CALL dgemm( 'T', 'N', jw, istopm-ihi, jw, one, qc, ldqc,
537  \$ b( kwtop, ihi+1 ), ldb, zero, work, jw )
538  CALL dlacpy( 'ALL', jw, istopm-ihi, work, jw, b( kwtop,
539  \$ ihi+1 ), ldb )
540  END IF
541  IF ( ilq ) THEN
542  CALL dgemm( 'N', 'N', n, jw, jw, one, q( 1, kwtop ), ldq, qc,
543  \$ ldqc, zero, work, n )
544  CALL dlacpy( 'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
545  END IF
546
547  IF ( kwtop-1-istartm+1 > 0 ) THEN
548  CALL dgemm( 'N', 'N', kwtop-istartm, jw, jw, one, a( istartm,
549  \$ kwtop ), lda, zc, ldzc, zero, work,
550  \$ kwtop-istartm )
551  CALL dlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
552  \$ a( istartm, kwtop ), lda )
553  CALL dgemm( 'N', 'N', kwtop-istartm, jw, jw, one, b( istartm,
554  \$ kwtop ), ldb, zc, ldzc, zero, work,
555  \$ kwtop-istartm )
556  CALL dlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
557  \$ b( istartm, kwtop ), ldb )
558  END IF
559  IF ( ilz ) THEN
560  CALL dgemm( 'N', 'N', n, jw, jw, one, z( 1, kwtop ), ldz, zc,
561  \$ ldzc, zero, work, n )
562  CALL dlacpy( 'ALL', n, jw, work, n, z( 1, kwtop ), ldz )
563  END IF
564
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f90:113
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:92
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
subroutine dtgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, WORK, LWORK, INFO)
DTGEXC
Definition: dtgexc.f:220
recursive subroutine dlaqz0(WANTS, WANTQ, WANTZ, N, ILO, IHI, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, REC, INFO)
DLAQZ0
Definition: dlaqz0.f:306
subroutine dlaqz2(ILQ, ILZ, K, ISTARTM, ISTOPM, IHI, A, LDA, B, LDB, NQ, QSTART, Q, LDQ, NZ, ZSTART, Z, LDZ)
DLAQZ2
Definition: dlaqz2.f:174
subroutine dlag2(A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2, WI)
DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
Definition: dlag2.f:156
Here is the call graph for this function:
Here is the caller graph for this function: