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