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

## ◆ slaqz3()

 recursive subroutine slaqz3 ( 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, real, dimension( lda, * ), intent(inout) a, integer, intent(in) lda, real, dimension( ldb, * ), intent(inout) b, integer, intent(in) ldb, real, dimension( ldq, * ), intent(inout) q, integer, intent(in) ldq, real, dimension( ldz, * ), intent(inout) z, integer, intent(in) ldz, integer, intent(out) ns, integer, intent(out) nd, real, dimension( * ), intent(inout) alphar, real, dimension( * ), intent(inout) alphai, real, dimension( * ), intent(inout) beta, real, dimension( ldqc, * ) qc, integer, intent(in) ldqc, real, dimension( ldzc, * ) zc, integer, intent(in) ldzc, real, dimension( * ) work, integer, intent(in) lwork, integer, intent(in) rec, integer, intent(out) info )

SLAQZ3

Purpose:
` SLAQZ3 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 REAL 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 REAL 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 REAL array, dimension (LDQ, N)` [in] LDQ ` LDQ is INTEGER` [in,out] Z ` Z is REAL 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 REAL array, dimension (N) The real parts of each scalar alpha defining an eigenvalue of GNEP.``` [out] ALPHAI ``` ALPHAI is REAL 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 REAL 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 REAL array, dimension (LDQC, NW)` [in] LDQC ` LDQC is INTEGER` [in,out] ZC ` ZC is REAL array, dimension (LDZC, NW)` [in] LDZC ` LDZ is INTEGER` [out] WORK ``` WORK is REAL 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 234 of file slaqz3.f.

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