LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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

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

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
Author
Thijs Steel, KU Leuven
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, 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 ulp = dlamch( 'PRECISION' )
306 smlnum = safmin*( dble( n )/ulp )
307
308 IF ( ihi .EQ. kwtop ) THEN
309* 1 by 1 deflation window, just try a regular deflation
310 alphar( kwtop ) = a( kwtop, kwtop )
311 alphai( kwtop ) = zero
312 beta( kwtop ) = b( kwtop, kwtop )
313 ns = 1
314 nd = 0
315 IF ( abs( s ) .LE. max( smlnum, ulp*abs( a( kwtop,
316 $ kwtop ) ) ) ) THEN
317 ns = 0
318 nd = 1
319 IF ( kwtop .GT. ilo ) THEN
320 a( kwtop, kwtop-1 ) = zero
321 END IF
322 END IF
323 END IF
324
325
326* Store window in case of convergence failure
327 CALL dlacpy( 'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
328 CALL dlacpy( 'ALL', jw, jw, b( kwtop, kwtop ), ldb, work( jw**2+
329 $ 1 ), jw )
330
331* Transform window to real schur form
332 CALL dlaset( 'FULL', jw, jw, zero, one, qc, ldqc )
333 CALL dlaset( 'FULL', jw, jw, zero, one, zc, ldzc )
334 CALL dlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
335 $ b( kwtop, kwtop ), ldb, alphar, alphai, beta, qc,
336 $ ldqc, zc, ldzc, work( 2*jw**2+1 ), lwork-2*jw**2,
337 $ rec+1, qz_small_info )
338
339 IF( qz_small_info .NE. 0 ) THEN
340* Convergence failure, restore the window and exit
341 nd = 0
342 ns = jw-qz_small_info
343 CALL dlacpy( 'ALL', jw, jw, work, jw, a( kwtop, kwtop ), lda )
344 CALL dlacpy( 'ALL', jw, jw, work( jw**2+1 ), jw, b( kwtop,
345 $ kwtop ), ldb )
346 RETURN
347 END IF
348
349* Deflation detection loop
350 IF ( kwtop .EQ. ilo .OR. s .EQ. zero ) THEN
351 kwbot = kwtop-1
352 ELSE
353 kwbot = ihi
354 k = 1
355 k2 = 1
356 DO WHILE ( k .LE. jw )
357 bulge = .false.
358 IF ( kwbot-kwtop+1 .GE. 2 ) THEN
359 bulge = a( kwbot, kwbot-1 ) .NE. zero
360 END IF
361 IF ( bulge ) THEN
362
363* Try to deflate complex conjugate eigenvalue pair
364 temp = abs( a( kwbot, kwbot ) )+sqrt( abs( a( kwbot,
365 $ kwbot-1 ) ) )*sqrt( abs( a( kwbot-1, kwbot ) ) )
366 IF( temp .EQ. zero )THEN
367 temp = abs( s )
368 END IF
369 IF ( max( abs( s*qc( 1, kwbot-kwtop ) ), abs( s*qc( 1,
370 $ kwbot-kwtop+1 ) ) ) .LE. max( smlnum,
371 $ ulp*temp ) ) THEN
372* Deflatable
373 kwbot = kwbot-2
374 ELSE
375* Not deflatable, move out of the way
376 ifst = kwbot-kwtop+1
377 ilst = k2
378 CALL dtgexc( .true., .true., jw, a( kwtop, kwtop ),
379 $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
380 $ zc, ldzc, ifst, ilst, work, lwork,
381 $ dtgexc_info )
382 k2 = k2+2
383 END IF
384 k = k+2
385 ELSE
386
387* Try to deflate real eigenvalue
388 temp = abs( a( kwbot, kwbot ) )
389 IF( temp .EQ. zero ) THEN
390 temp = abs( s )
391 END IF
392 IF ( ( abs( s*qc( 1, kwbot-kwtop+1 ) ) ) .LE. max( ulp*
393 $ temp, smlnum ) ) THEN
394* Deflatable
395 kwbot = kwbot-1
396 ELSE
397* Not deflatable, move out of the way
398 ifst = kwbot-kwtop+1
399 ilst = k2
400 CALL dtgexc( .true., .true., jw, a( kwtop, kwtop ),
401 $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
402 $ zc, ldzc, ifst, ilst, work, lwork,
403 $ dtgexc_info )
404 k2 = k2+1
405 END IF
406
407 k = k+1
408
409 END IF
410 END DO
411 END IF
412
413* Store eigenvalues
414 nd = ihi-kwbot
415 ns = jw-nd
416 k = kwtop
417 DO WHILE ( k .LE. ihi )
418 bulge = .false.
419 IF ( k .LT. ihi ) THEN
420 IF ( a( k+1, k ) .NE. zero ) THEN
421 bulge = .true.
422 END IF
423 END IF
424 IF ( bulge ) THEN
425* 2x2 eigenvalue block
426 CALL dlag2( a( k, k ), lda, b( k, k ), ldb, safmin,
427 $ beta( k ), beta( k+1 ), alphar( k ),
428 $ alphar( k+1 ), alphai( k ) )
429 alphai( k+1 ) = -alphai( k )
430 k = k+2
431 ELSE
432* 1x1 eigenvalue block
433 alphar( k ) = a( k, k )
434 alphai( k ) = zero
435 beta( k ) = b( k, k )
436 k = k+1
437 END IF
438 END DO
439
440 IF ( kwtop .NE. ilo .AND. s .NE. zero ) THEN
441* Reflect spike back, this will create optimally packed bulges
442 a( kwtop:kwbot, kwtop-1 ) = a( kwtop, kwtop-1 )*qc( 1,
443 $ 1:jw-nd )
444 DO k = kwbot-1, kwtop, -1
445 CALL dlartg( a( k, kwtop-1 ), a( k+1, kwtop-1 ), c1, s1,
446 $ temp )
447 a( k, kwtop-1 ) = temp
448 a( k+1, kwtop-1 ) = zero
449 k2 = max( kwtop, k-1 )
450 CALL drot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda, c1,
451 $ s1 )
452 CALL drot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1, k-1 ),
453 $ ldb, c1, s1 )
454 CALL drot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1, k+1-kwtop+1 ),
455 $ 1, c1, s1 )
456 END DO
457
458* Chase bulges down
459 istartm = kwtop
460 istopm = ihi
461 k = kwbot-1
462 DO WHILE ( k .GE. kwtop )
463 IF ( ( k .GE. kwtop+1 ) .AND. a( k+1, k-1 ) .NE. zero ) THEN
464
465* Move double pole block down and remove it
466 DO k2 = k-1, kwbot-2
467 CALL dlaqz2( .true., .true., k2, kwtop, kwtop+jw-1,
468 $ kwbot, a, lda, b, ldb, jw, kwtop, qc,
469 $ ldqc, jw, kwtop, zc, ldzc )
470 END DO
471
472 k = k-2
473 ELSE
474
475* k points to single shift
476 DO k2 = k, kwbot-2
477
478* Move shift down
479 CALL dlartg( b( k2+1, k2+1 ), b( k2+1, k2 ), c1, s1,
480 $ temp )
481 b( k2+1, k2+1 ) = temp
482 b( k2+1, k2 ) = zero
483 CALL drot( k2+2-istartm+1, a( istartm, k2+1 ), 1,
484 $ a( istartm, k2 ), 1, c1, s1 )
485 CALL drot( k2-istartm+1, b( istartm, k2+1 ), 1,
486 $ b( istartm, k2 ), 1, c1, s1 )
487 CALL drot( jw, zc( 1, k2+1-kwtop+1 ), 1, zc( 1,
488 $ k2-kwtop+1 ), 1, c1, s1 )
489
490 CALL dlartg( a( k2+1, k2 ), a( k2+2, k2 ), c1, s1,
491 $ temp )
492 a( k2+1, k2 ) = temp
493 a( k2+2, k2 ) = zero
494 CALL drot( istopm-k2, a( k2+1, k2+1 ), lda, a( k2+2,
495 $ k2+1 ), lda, c1, s1 )
496 CALL drot( istopm-k2, b( k2+1, k2+1 ), ldb, b( k2+2,
497 $ k2+1 ), ldb, c1, s1 )
498 CALL drot( jw, qc( 1, k2+1-kwtop+1 ), 1, qc( 1,
499 $ k2+2-kwtop+1 ), 1, c1, s1 )
500
501 END DO
502
503* Remove the shift
504 CALL dlartg( b( kwbot, kwbot ), b( kwbot, kwbot-1 ), c1,
505 $ s1, temp )
506 b( kwbot, kwbot ) = temp
507 b( kwbot, kwbot-1 ) = zero
508 CALL drot( kwbot-istartm, b( istartm, kwbot ), 1,
509 $ b( istartm, kwbot-1 ), 1, c1, s1 )
510 CALL drot( kwbot-istartm+1, a( istartm, kwbot ), 1,
511 $ a( istartm, kwbot-1 ), 1, c1, s1 )
512 CALL drot( jw, zc( 1, kwbot-kwtop+1 ), 1, zc( 1,
513 $ kwbot-1-kwtop+1 ), 1, c1, s1 )
514
515 k = k-1
516 END IF
517 END DO
518
519 END IF
520
521* Apply Qc and Zc to rest of the matrix
522 IF ( ilschur ) THEN
523 istartm = 1
524 istopm = n
525 ELSE
526 istartm = ilo
527 istopm = ihi
528 END IF
529
530 IF ( istopm-ihi > 0 ) THEN
531 CALL dgemm( 'T', 'N', jw, istopm-ihi, jw, one, qc, ldqc,
532 $ a( kwtop, ihi+1 ), lda, zero, work, jw )
533 CALL dlacpy( 'ALL', jw, istopm-ihi, work, jw, a( kwtop,
534 $ ihi+1 ), lda )
535 CALL dgemm( 'T', 'N', jw, istopm-ihi, jw, one, qc, ldqc,
536 $ b( kwtop, ihi+1 ), ldb, zero, work, jw )
537 CALL dlacpy( 'ALL', jw, istopm-ihi, work, jw, b( kwtop,
538 $ ihi+1 ), ldb )
539 END IF
540 IF ( ilq ) THEN
541 CALL dgemm( 'N', 'N', n, jw, jw, one, q( 1, kwtop ), ldq, qc,
542 $ ldqc, zero, work, n )
543 CALL dlacpy( 'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
544 END IF
545
546 IF ( kwtop-1-istartm+1 > 0 ) THEN
547 CALL dgemm( 'N', 'N', kwtop-istartm, jw, jw, one, a( istartm,
548 $ kwtop ), lda, zc, ldzc, zero, work,
549 $ kwtop-istartm )
550 CALL dlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
551 $ a( istartm, kwtop ), lda )
552 CALL dgemm( 'N', 'N', kwtop-istartm, jw, jw, one, b( istartm,
553 $ kwtop ), ldb, zc, ldzc, zero, work,
554 $ kwtop-istartm )
555 CALL dlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
556 $ b( istartm, kwtop ), ldb )
557 END IF
558 IF ( ilz ) THEN
559 CALL dgemm( 'N', 'N', n, jw, jw, one, z( 1, kwtop ), ldz, zc,
560 $ ldzc, zero, work, n )
561 CALL dlacpy( 'ALL', n, jw, work, n, z( 1, kwtop ), ldz )
562 END IF
563
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
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 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
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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 dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
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 drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
subroutine dtgexc(wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, ifst, ilst, work, lwork, info)
DTGEXC
Definition dtgexc.f:220
Here is the call graph for this function:
Here is the caller graph for this function: