LAPACK 3.12.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, 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 ulp = dlamch( 'PRECISION' )
301 smlnum = safmin*( dble( n )/ulp )
302
303 IF ( ihi .EQ. kwtop ) THEN
304* 1 by 1 deflation window, just try a regular deflation
305 alpha( kwtop ) = a( kwtop, kwtop )
306 beta( kwtop ) = b( kwtop, kwtop )
307 ns = 1
308 nd = 0
309 IF ( abs( s ) .LE. max( smlnum, ulp*abs( a( kwtop,
310 $ kwtop ) ) ) ) THEN
311 ns = 0
312 nd = 1
313 IF ( kwtop .GT. ilo ) THEN
314 a( kwtop, kwtop-1 ) = czero
315 END IF
316 END IF
317 END IF
318
319
320* Store window in case of convergence failure
321 CALL zlacpy( 'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
322 CALL zlacpy( 'ALL', jw, jw, b( kwtop, kwtop ), ldb, work( jw**2+
323 $ 1 ), jw )
324
325* Transform window to real schur form
326 CALL zlaset( 'FULL', jw, jw, czero, cone, qc, ldqc )
327 CALL zlaset( 'FULL', jw, jw, czero, cone, zc, ldzc )
328 CALL zlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
329 $ b( kwtop, kwtop ), ldb, alpha, beta, qc, ldqc, zc,
330 $ ldzc, work( 2*jw**2+1 ), lwork-2*jw**2, rwork,
331 $ rec+1, qz_small_info )
332
333 IF( qz_small_info .NE. 0 ) THEN
334* Convergence failure, restore the window and exit
335 nd = 0
336 ns = jw-qz_small_info
337 CALL zlacpy( 'ALL', jw, jw, work, jw, a( kwtop, kwtop ), lda )
338 CALL zlacpy( 'ALL', jw, jw, work( jw**2+1 ), jw, b( kwtop,
339 $ kwtop ), ldb )
340 RETURN
341 END IF
342
343* Deflation detection loop
344 IF ( kwtop .EQ. ilo .OR. s .EQ. czero ) THEN
345 kwbot = kwtop-1
346 ELSE
347 kwbot = ihi
348 k = 1
349 k2 = 1
350 DO WHILE ( k .LE. jw )
351* Try to deflate eigenvalue
352 tempr = abs( a( kwbot, kwbot ) )
353 IF( tempr .EQ. zero ) THEN
354 tempr = abs( s )
355 END IF
356 IF ( ( abs( s*qc( 1, kwbot-kwtop+1 ) ) ) .LE. max( ulp*
357 $ tempr, smlnum ) ) THEN
358* Deflatable
359 kwbot = kwbot-1
360 ELSE
361* Not deflatable, move out of the way
362 ifst = kwbot-kwtop+1
363 ilst = k2
364 CALL ztgexc( .true., .true., jw, a( kwtop, kwtop ),
365 $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
366 $ zc, ldzc, ifst, ilst, ztgexc_info )
367 k2 = k2+1
368 END IF
369
370 k = k+1
371 END DO
372 END IF
373
374* Store eigenvalues
375 nd = ihi-kwbot
376 ns = jw-nd
377 k = kwtop
378 DO WHILE ( k .LE. ihi )
379 alpha( k ) = a( k, k )
380 beta( k ) = b( k, k )
381 k = k+1
382 END DO
383
384 IF ( kwtop .NE. ilo .AND. s .NE. czero ) THEN
385* Reflect spike back, this will create optimally packed bulges
386 a( kwtop:kwbot, kwtop-1 ) = a( kwtop, kwtop-1 ) *dconjg( qc( 1,
387 $ 1:jw-nd ) )
388 DO k = kwbot-1, kwtop, -1
389 CALL zlartg( a( k, kwtop-1 ), a( k+1, kwtop-1 ), c1, s1,
390 $ temp )
391 a( k, kwtop-1 ) = temp
392 a( k+1, kwtop-1 ) = czero
393 k2 = max( kwtop, k-1 )
394 CALL zrot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda, c1,
395 $ s1 )
396 CALL zrot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1, k-1 ),
397 $ ldb, c1, s1 )
398 CALL zrot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1, k+1-kwtop+1 ),
399 $ 1, c1, dconjg( s1 ) )
400 END DO
401
402* Chase bulges down
403 istartm = kwtop
404 istopm = ihi
405 k = kwbot-1
406 DO WHILE ( k .GE. kwtop )
407
408* Move bulge down and remove it
409 DO k2 = k, kwbot-1
410 CALL zlaqz1( .true., .true., k2, kwtop, kwtop+jw-1,
411 $ kwbot, a, lda, b, ldb, jw, kwtop, qc, ldqc,
412 $ jw, kwtop, zc, ldzc )
413 END DO
414
415 k = k-1
416 END DO
417
418 END IF
419
420* Apply Qc and Zc to rest of the matrix
421 IF ( ilschur ) THEN
422 istartm = 1
423 istopm = n
424 ELSE
425 istartm = ilo
426 istopm = ihi
427 END IF
428
429 IF ( istopm-ihi > 0 ) THEN
430 CALL zgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
431 $ a( kwtop, ihi+1 ), lda, czero, work, jw )
432 CALL zlacpy( 'ALL', jw, istopm-ihi, work, jw, a( kwtop,
433 $ ihi+1 ), lda )
434 CALL zgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
435 $ b( kwtop, ihi+1 ), ldb, czero, work, jw )
436 CALL zlacpy( 'ALL', jw, istopm-ihi, work, jw, b( kwtop,
437 $ ihi+1 ), ldb )
438 END IF
439 IF ( ilq ) THEN
440 CALL zgemm( 'N', 'N', n, jw, jw, cone, q( 1, kwtop ), ldq, qc,
441 $ ldqc, czero, work, n )
442 CALL zlacpy( 'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
443 END IF
444
445 IF ( kwtop-1-istartm+1 > 0 ) THEN
446 CALL zgemm( 'N', 'N', kwtop-istartm, jw, jw, cone, a( istartm,
447 $ kwtop ), lda, zc, ldzc, czero, work,
448 $ kwtop-istartm )
449 CALL zlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
450 $ a( istartm, kwtop ), lda )
451 CALL zgemm( 'N', 'N', kwtop-istartm, jw, jw, cone, b( istartm,
452 $ kwtop ), ldb, zc, ldzc, czero, work,
453 $ kwtop-istartm )
454 CALL zlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
455 $ b( istartm, kwtop ), ldb )
456 END IF
457 IF ( ilz ) THEN
458 CALL zgemm( 'N', 'N', n, jw, jw, cone, z( 1, kwtop ), ldz, zc,
459 $ ldzc, czero, work, n )
460 CALL zlacpy( 'ALL', n, jw, work, n, z( 1, kwtop ), ldz )
461 END IF
462
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
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
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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
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
zlartg
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition zlartg.f90:116
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
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
ztgexc
subroutine ztgexc(wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, ifst, ilst, info)
ZTGEXC
Definition ztgexc.f:200
Here is the call graph for this function:
Here is the caller graph for this function: