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

◆ claqz2()

recursive subroutine claqz2 ( 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, dimension( lda, * ), intent(inout)  a,
integer, intent(in)  lda,
complex, dimension( ldb, * ), intent(inout)  b,
integer, intent(in)  ldb,
complex, dimension( ldq, * ), intent(inout)  q,
integer, intent(in)  ldq,
complex, dimension( ldz, * ), intent(inout)  z,
integer, intent(in)  ldz,
integer, intent(out)  ns,
integer, intent(out)  nd,
complex, dimension( * ), intent(inout)  alpha,
complex, dimension( * ), intent(inout)  beta,
complex, dimension( ldqc, * )  qc,
integer, intent(in)  ldqc,
complex, dimension( ldzc, * )  zc,
integer, intent(in)  ldzc,
complex, dimension( * )  work,
integer, intent(in)  lwork,
real, dimension( * )  rwork,
integer, intent(in)  rec,
integer, intent(out)  info 
)

CLAQZ2

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

Purpose:
 CLAQZ2 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 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 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 array, dimension (LDQ, N)
[in]LDQ
          LDQ is INTEGER
[in,out]Z
          Z is COMPLEX 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 array, dimension (N)
          Each scalar alpha defining an eigenvalue
          of GNEP.
[out]BETA
          BETA is COMPLEX 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 array, dimension (LDQC, NW)
[in]LDQC
          LDQC is INTEGER
[in,out]ZC
          ZC is COMPLEX array, dimension (LDZC, NW)
[in]LDZC
          LDZ is INTEGER
[out]WORK
          WORK is COMPLEX 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 REAL 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, KU Leuven
Date
May 2020

Definition at line 231 of file claqz2.f.

235 IMPLICIT NONE
236
237* Arguments
238 LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
239 INTEGER, INTENT( IN ) :: N, ILO, IHI, NW, LDA, LDB, LDQ, LDZ,
240 $ LDQC, LDZC, LWORK, REC
241
242 COMPLEX, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
243 $ Z( LDZ, * ), ALPHA( * ), BETA( * )
244 INTEGER, INTENT( OUT ) :: NS, ND, INFO
245 COMPLEX :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
246 REAL :: RWORK( * )
247
248* Parameters
249 COMPLEX CZERO, CONE
250 parameter( czero = ( 0.0, 0.0 ), cone = ( 1.0, 0.0 ) )
251 REAL :: ZERO, ONE, HALF
252 parameter( zero = 0.0, one = 1.0, half = 0.5 )
253
254* Local Scalars
255 INTEGER :: JW, KWTOP, KWBOT, ISTOPM, ISTARTM, K, K2, CTGEXC_INFO,
256 $ IFST, ILST, LWORKREQ, QZ_SMALL_INFO
257 REAL :: SMLNUM, ULP, SAFMIN, SAFMAX, C1, TEMPR
258 COMPLEX :: S, S1, TEMP
259
260* External Functions
261 EXTERNAL :: xerbla, claqz0, claqz1, clacpy, claset, cgemm,
262 $ ctgexc, clartg, crot
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 = 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 claqz0( '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( 'CLAQZ2', -info )
294 RETURN
295 END IF
296
297* Get machine constants
298 safmin = slamch( 'SAFE MINIMUM' )
299 safmax = one/safmin
300 ulp = slamch( 'PRECISION' )
301 smlnum = safmin*( real( 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 clacpy( 'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
322 CALL clacpy( 'ALL', jw, jw, b( kwtop, kwtop ), ldb, work( jw**2+
323 $ 1 ), jw )
324
325* Transform window to real schur form
326 CALL claset( 'FULL', jw, jw, czero, cone, qc, ldqc )
327 CALL claset( 'FULL', jw, jw, czero, cone, zc, ldzc )
328 CALL claqz0( '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 clacpy( 'ALL', jw, jw, work, jw, a( kwtop, kwtop ), lda )
338 CALL clacpy( '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 ctgexc( .true., .true., jw, a( kwtop, kwtop ),
365 $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
366 $ zc, ldzc, ifst, ilst, ctgexc_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 ) *conjg( qc( 1,
387 $ 1:jw-nd ) )
388 DO k = kwbot-1, kwtop, -1
389 CALL clartg( 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 crot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda, c1,
395 $ s1 )
396 CALL crot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1, k-1 ),
397 $ ldb, c1, s1 )
398 CALL crot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1, k+1-kwtop+1 ),
399 $ 1, c1, conjg( 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 claqz1( .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 cgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
431 $ a( kwtop, ihi+1 ), lda, czero, work, jw )
432 CALL clacpy( 'ALL', jw, istopm-ihi, work, jw, a( kwtop,
433 $ ihi+1 ), lda )
434 CALL cgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
435 $ b( kwtop, ihi+1 ), ldb, czero, work, jw )
436 CALL clacpy( 'ALL', jw, istopm-ihi, work, jw, b( kwtop,
437 $ ihi+1 ), ldb )
438 END IF
439 IF ( ilq ) THEN
440 CALL cgemm( 'N', 'N', n, jw, jw, cone, q( 1, kwtop ), ldq, qc,
441 $ ldqc, czero, work, n )
442 CALL clacpy( 'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
443 END IF
444
445 IF ( kwtop-1-istartm+1 > 0 ) THEN
446 CALL cgemm( 'N', 'N', kwtop-istartm, jw, jw, cone, a( istartm,
447 $ kwtop ), lda, zc, ldzc, czero, work,
448 $ kwtop-istartm )
449 CALL clacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
450 $ a( istartm, kwtop ), lda )
451 CALL cgemm( 'N', 'N', kwtop-istartm, jw, jw, cone, b( istartm,
452 $ kwtop ), ldb, zc, ldzc, czero, work,
453 $ kwtop-istartm )
454 CALL clacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
455 $ b( istartm, kwtop ), ldb )
456 END IF
457 IF ( ilz ) THEN
458 CALL cgemm( 'N', 'N', n, jw, jw, cone, z( 1, kwtop ), ldz, zc,
459 $ ldzc, czero, work, n )
460 CALL clacpy( 'ALL', n, jw, work, n, z( 1, kwtop ), ldz )
461 END IF
462
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:103
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
claqz0
recursive subroutine claqz0(wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb, alpha, beta, q, ldq, z, ldz, work, lwork, rwork, rec, info)
CLAQZ0
Definition claqz0.f:284
claqz1
subroutine claqz1(ilq, ilz, k, istartm, istopm, ihi, a, lda, b, ldb, nq, qstart, q, ldq, nz, zstart, z, ldz)
CLAQZ1
Definition claqz1.f:173
clartg
subroutine clartg(f, g, c, s, r)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition clartg.f90:116
claset
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition claset.f:106
crot
subroutine crot(n, cx, incx, cy, incy, c, s)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition crot.f:103
ctgexc
subroutine ctgexc(wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, ifst, ilst, info)
CTGEXC
Definition ctgexc.f:200
Here is the call graph for this function:
Here is the caller graph for this function: