LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
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

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

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
Author
Thijs Steel, KU Leuven
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 slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
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: