LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zlaqz2()

recursive subroutine zlaqz2 ( 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,
complex*16, dimension( lda, * ), intent(inout)  A,
integer, intent(in)  LDA,
complex*16, dimension( ldb, * ), intent(inout)  B,
integer, intent(in)  LDB,
complex*16, dimension( ldq, * ), intent(inout)  Q,
integer, intent(in)  LDQ,
complex*16, dimension( ldz, * ), intent(inout)  Z,
integer, intent(in)  LDZ,
integer, intent(out)  NS,
integer, intent(out)  ND,
complex*16, dimension( * ), intent(inout)  ALPHA,
complex*16, dimension( * ), intent(inout)  BETA,
complex*16, dimension( ldqc, * )  QC,
integer, intent(in)  LDQC,
complex*16, dimension( ldzc, * )  ZC,
integer, intent(in)  LDZC,
complex*16, dimension( * )  WORK,
integer, intent(in)  LWORK,
double precision, dimension( * )  RWORK,
integer, intent(in)  REC,
integer, intent(out)  INFO 
)

ZLAQZ2

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

Purpose:
 ZLAQZ2 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 COMPLEX*16 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 COMPLEX*16 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 COMPLEX*16 array, dimension (LDQ, N)
[in]LDQ
          LDQ is INTEGER
[in,out]Z
          Z is COMPLEX*16 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]ALPHA
          ALPHA is COMPLEX*16 array, dimension (N)
          Each scalar alpha defining an eigenvalue
          of GNEP.
[out]BETA
          BETA is COMPLEX*16 array, dimension (N)
          The scalars beta that define the eigenvalues of GNEP.
          Together, the quantities alpha = ALPHA(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 COMPLEX*16 array, dimension (LDQC, NW)
[in]LDQC
          LDQC is INTEGER
[in,out]ZC
          ZC is COMPLEX*16 array, dimension (LDZC, NW)
[in]LDZC
          LDZ is INTEGER
[out]WORK
          WORK is COMPLEX*16 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.
[out]RWORK
          RWORK is DOUBLE PRECISION array, dimension (N)
[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 230 of file zlaqz2.f.

234 IMPLICIT NONE
235
236* Arguments
237 LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
238 INTEGER, INTENT( IN ) :: N, ILO, IHI, NW, LDA, LDB, LDQ, LDZ,
239 $ LDQC, LDZC, LWORK, REC
240
241 COMPLEX*16, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ,
242 $ * ), Z( LDZ, * ), ALPHA( * ), BETA( * )
243 INTEGER, INTENT( OUT ) :: NS, ND, INFO
244 COMPLEX*16 :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
245 DOUBLE PRECISION :: RWORK( * )
246
247* Parameters
248 COMPLEX*16 CZERO, CONE
249 parameter( czero = ( 0.0d+0, 0.0d+0 ), cone = ( 1.0d+0,
250 $ 0.0d+0 ) )
251 DOUBLE PRECISION :: ZERO, ONE, HALF
252 parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
253
254* Local Scalars
255 INTEGER :: JW, KWTOP, KWBOT, ISTOPM, ISTARTM, K, K2, ZTGEXC_INFO,
256 $ IFST, ILST, LWORKREQ, QZ_SMALL_INFO
257 DOUBLE PRECISION ::SMLNUM, ULP, SAFMIN, SAFMAX, C1, TEMPR
258 COMPLEX*16 :: S, S1, TEMP
259
260* External Functions
261 EXTERNAL :: xerbla, zlaqz0, zlaqz1, dlabad, zlacpy, zlaset, zgemm,
262 $ ztgexc, zlartg, zrot
263 DOUBLE PRECISION, EXTERNAL :: DLAMCH
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 = czero
272 ELSE
273 s = a( kwtop, kwtop-1 )
274 END IF
275
276* Determine required workspace
277 ifst = 1
278 ilst = jw
279 CALL zlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
280 $ b( kwtop, kwtop ), ldb, alpha, beta, qc, ldqc, zc,
281 $ ldzc, work, -1, rwork, rec+1, qz_small_info )
282 lworkreq = int( work( 1 ) )+2*jw**2
283 lworkreq = max( lworkreq, n*nw, 2*nw**2+n )
284 IF ( lwork .EQ.-1 ) THEN
285* workspace query, quick return
286 work( 1 ) = lworkreq
287 RETURN
288 ELSE IF ( lwork .LT. lworkreq ) THEN
289 info = -26
290 END IF
291
292 IF( info.NE.0 ) THEN
293 CALL xerbla( 'ZLAQZ2', -info )
294 RETURN
295 END IF
296
297* Get machine constants
298 safmin = dlamch( 'SAFE MINIMUM' )
299 safmax = one/safmin
300 CALL dlabad( safmin, safmax )
301 ulp = dlamch( 'PRECISION' )
302 smlnum = safmin*( dble( n )/ulp )
303
304 IF ( ihi .EQ. kwtop ) THEN
305* 1 by 1 deflation window, just try a regular deflation
306 alpha( kwtop ) = a( kwtop, kwtop )
307 beta( kwtop ) = b( kwtop, kwtop )
308 ns = 1
309 nd = 0
310 IF ( abs( s ) .LE. max( smlnum, ulp*abs( a( kwtop,
311 $ kwtop ) ) ) ) THEN
312 ns = 0
313 nd = 1
314 IF ( kwtop .GT. ilo ) THEN
315 a( kwtop, kwtop-1 ) = czero
316 END IF
317 END IF
318 END IF
319
320
321* Store window in case of convergence failure
322 CALL zlacpy( 'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
323 CALL zlacpy( 'ALL', jw, jw, b( kwtop, kwtop ), ldb, work( jw**2+
324 $ 1 ), jw )
325
326* Transform window to real schur form
327 CALL zlaset( 'FULL', jw, jw, czero, cone, qc, ldqc )
328 CALL zlaset( 'FULL', jw, jw, czero, cone, zc, ldzc )
329 CALL zlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
330 $ b( kwtop, kwtop ), ldb, alpha, beta, qc, ldqc, zc,
331 $ ldzc, work( 2*jw**2+1 ), lwork-2*jw**2, rwork,
332 $ rec+1, qz_small_info )
333
334 IF( qz_small_info .NE. 0 ) THEN
335* Convergence failure, restore the window and exit
336 nd = 0
337 ns = jw-qz_small_info
338 CALL zlacpy( 'ALL', jw, jw, work, jw, a( kwtop, kwtop ), lda )
339 CALL zlacpy( 'ALL', jw, jw, work( jw**2+1 ), jw, b( kwtop,
340 $ kwtop ), ldb )
341 RETURN
342 END IF
343
344* Deflation detection loop
345 IF ( kwtop .EQ. ilo .OR. s .EQ. czero ) THEN
346 kwbot = kwtop-1
347 ELSE
348 kwbot = ihi
349 k = 1
350 k2 = 1
351 DO WHILE ( k .LE. jw )
352* Try to deflate eigenvalue
353 tempr = abs( a( kwbot, kwbot ) )
354 IF( tempr .EQ. zero ) THEN
355 tempr = abs( s )
356 END IF
357 IF ( ( abs( s*qc( 1, kwbot-kwtop+1 ) ) ) .LE. max( ulp*
358 $ tempr, smlnum ) ) THEN
359* Deflatable
360 kwbot = kwbot-1
361 ELSE
362* Not deflatable, move out of the way
363 ifst = kwbot-kwtop+1
364 ilst = k2
365 CALL ztgexc( .true., .true., jw, a( kwtop, kwtop ),
366 $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
367 $ zc, ldzc, ifst, ilst, ztgexc_info )
368 k2 = k2+1
369 END IF
370
371 k = k+1
372 END DO
373 END IF
374
375* Store eigenvalues
376 nd = ihi-kwbot
377 ns = jw-nd
378 k = kwtop
379 DO WHILE ( k .LE. ihi )
380 alpha( k ) = a( k, k )
381 beta( k ) = b( k, k )
382 k = k+1
383 END DO
384
385 IF ( kwtop .NE. ilo .AND. s .NE. czero ) THEN
386* Reflect spike back, this will create optimally packed bulges
387 a( kwtop:kwbot, kwtop-1 ) = a( kwtop, kwtop-1 ) *dconjg( qc( 1,
388 $ 1:jw-nd ) )
389 DO k = kwbot-1, kwtop, -1
390 CALL zlartg( a( k, kwtop-1 ), a( k+1, kwtop-1 ), c1, s1,
391 $ temp )
392 a( k, kwtop-1 ) = temp
393 a( k+1, kwtop-1 ) = czero
394 k2 = max( kwtop, k-1 )
395 CALL zrot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda, c1,
396 $ s1 )
397 CALL zrot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1, k-1 ),
398 $ ldb, c1, s1 )
399 CALL zrot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1, k+1-kwtop+1 ),
400 $ 1, c1, dconjg( s1 ) )
401 END DO
402
403* Chase bulges down
404 istartm = kwtop
405 istopm = ihi
406 k = kwbot-1
407 DO WHILE ( k .GE. kwtop )
408
409* Move bulge down and remove it
410 DO k2 = k, kwbot-1
411 CALL zlaqz1( .true., .true., k2, kwtop, kwtop+jw-1,
412 $ kwbot, a, lda, b, ldb, jw, kwtop, qc, ldqc,
413 $ jw, kwtop, zc, ldzc )
414 END DO
415
416 k = k-1
417 END DO
418
419 END IF
420
421* Apply Qc and Zc to rest of the matrix
422 IF ( ilschur ) THEN
423 istartm = 1
424 istopm = n
425 ELSE
426 istartm = ilo
427 istopm = ihi
428 END IF
429
430 IF ( istopm-ihi > 0 ) THEN
431 CALL zgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
432 $ a( kwtop, ihi+1 ), lda, czero, work, jw )
433 CALL zlacpy( 'ALL', jw, istopm-ihi, work, jw, a( kwtop,
434 $ ihi+1 ), lda )
435 CALL zgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
436 $ b( kwtop, ihi+1 ), ldb, czero, work, jw )
437 CALL zlacpy( 'ALL', jw, istopm-ihi, work, jw, b( kwtop,
438 $ ihi+1 ), ldb )
439 END IF
440 IF ( ilq ) THEN
441 CALL zgemm( 'N', 'N', n, jw, jw, cone, q( 1, kwtop ), ldq, qc,
442 $ ldqc, czero, work, n )
443 CALL zlacpy( 'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
444 END IF
445
446 IF ( kwtop-1-istartm+1 > 0 ) THEN
447 CALL zgemm( 'N', 'N', kwtop-istartm, jw, jw, cone, a( istartm,
448 $ kwtop ), lda, zc, ldzc, czero, work,
449 $ kwtop-istartm )
450 CALL zlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
451 $ a( istartm, kwtop ), lda )
452 CALL zgemm( 'N', 'N', kwtop-istartm, jw, jw, cone, b( istartm,
453 $ kwtop ), ldb, zc, ldzc, czero, work,
454 $ kwtop-istartm )
455 CALL zlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
456 $ b( istartm, kwtop ), ldb )
457 END IF
458 IF ( ilz ) THEN
459 CALL zgemm( 'N', 'N', n, jw, jw, cone, z( 1, kwtop ), ldz, zc,
460 $ ldzc, czero, work, n )
461 CALL zlacpy( 'ALL', n, jw, work, n, z( 1, kwtop ), ldz )
462 END IF
463
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
zlartg
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition: zlartg.f90:116
dlabad
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
zgemm
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
zlaqz1
subroutine zlaqz1(ILQ, ILZ, K, ISTARTM, ISTOPM, IHI, A, LDA, B, LDB, NQ, QSTART, Q, LDQ, NZ, ZSTART, Z, LDZ)
ZLAQZ1
Definition: zlaqz1.f:173
zlaqz0
recursive subroutine zlaqz0(WANTS, WANTQ, WANTZ, N, ILO, IHI, A, LDA, B, LDB, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, REC, INFO)
ZLAQZ0
Definition: zlaqz0.f:284
ztgexc
subroutine ztgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, INFO)
ZTGEXC
Definition: ztgexc.f:200
zrot
subroutine zrot(N, CX, INCX, CY, INCY, C, S)
ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition: zrot.f:103
zlacpy
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
zlaset
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:106
Here is the call graph for this function:
Here is the caller graph for this function: