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

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```
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
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition: zlartg.f90:116
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
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
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
subroutine ztgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, INFO)
ZTGEXC
Definition: ztgexc.f:200
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
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
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: