LAPACK 3.11.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, 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 dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
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: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 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: