LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
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

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, 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: