LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ zlaqz3()

subroutine zlaqz3 ( 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)  NSHIFTS,
integer, intent(in)  NBLOCK_DESIRED,
complex*16, dimension( * ), intent(inout)  ALPHA,
complex*16, dimension( * ), intent(inout)  BETA,
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,
complex*16, dimension( ldqc, * ), intent(inout)  QC,
integer, intent(in)  LDQC,
complex*16, dimension( ldzc, * ), intent(inout)  ZC,
integer, intent(in)  LDZC,
complex*16, dimension( * ), intent(inout)  WORK,
integer, intent(in)  LWORK,
integer, intent(out)  INFO 
)

ZLAQZ3

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

Purpose:
 ZLAQZ3 Executes a single multishift QZ sweep
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
[in]NSHIFTS
          NSHIFTS is INTEGER
          The desired number of shifts to use
[in]NBLOCK_DESIRED
          NBLOCK_DESIRED is INTEGER
          The desired size of the computational windows
[in]ALPHA
          ALPHA is COMPLEX*16 array. SR contains
          the alpha parts of the shifts to use.
[in]BETA
          BETA is COMPLEX*16 array. SS contains
          the scale of the shifts to use.
[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
[in,out]QC
          QC is COMPLEX*16 array, dimension (LDQC, NBLOCK_DESIRED)
[in]LDQC
          LDQC is INTEGER
[in,out]ZC
          ZC is COMPLEX*16 array, dimension (LDZC, NBLOCK_DESIRED)
[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]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 204 of file zlaqz3.f.

208  IMPLICIT NONE
209 
210 * Function arguments
211  LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
212  INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
213  $ NSHIFTS, NBLOCK_DESIRED, LDQC, LDZC
214 
215  COMPLEX*16, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ,
216  $ * ), Z( LDZ, * ), QC( LDQC, * ), ZC( LDZC, * ), WORK( * ),
217  $ ALPHA( * ), BETA( * )
218 
219  INTEGER, INTENT( OUT ) :: INFO
220 
221 * Parameters
222  COMPLEX*16 CZERO, CONE
223  parameter( czero = ( 0.0d+0, 0.0d+0 ), cone = ( 1.0d+0,
224  $ 0.0d+0 ) )
225  DOUBLE PRECISION :: ZERO, ONE, HALF
226  parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
227 
228 * Local scalars
229  INTEGER :: I, J, NS, ISTARTM, ISTOPM, SHEIGHT, SWIDTH, K, NP,
230  $ ISTARTB, ISTOPB, ISHIFT, NBLOCK, NPOS
231  DOUBLE PRECISION :: SAFMIN, SAFMAX, C, SCALE
232  COMPLEX*16 :: TEMP, TEMP2, TEMP3, S
233 
234 * External Functions
235  EXTERNAL :: xerbla, dlabad, zlaset, zlartg, zrot, zlaqz1, zgemm,
236  $ zlacpy
237  DOUBLE PRECISION, EXTERNAL :: DLAMCH
238 
239  info = 0
240  IF ( nblock_desired .LT. nshifts+1 ) THEN
241  info = -8
242  END IF
243  IF ( lwork .EQ.-1 ) THEN
244 * workspace query, quick return
245  work( 1 ) = n*nblock_desired
246  RETURN
247  ELSE IF ( lwork .LT. n*nblock_desired ) THEN
248  info = -25
249  END IF
250 
251  IF( info.NE.0 ) THEN
252  CALL xerbla( 'ZLAQZ3', -info )
253  RETURN
254  END IF
255 
256 *
257 * Executable statements
258 *
259 
260 * Get machine constants
261  safmin = dlamch( 'SAFE MINIMUM' )
262  safmax = one/safmin
263  CALL dlabad( safmin, safmax )
264 
265  IF ( ilo .GE. ihi ) THEN
266  RETURN
267  END IF
268 
269  IF ( ilschur ) THEN
270  istartm = 1
271  istopm = n
272  ELSE
273  istartm = ilo
274  istopm = ihi
275  END IF
276 
277  ns = nshifts
278  npos = max( nblock_desired-ns, 1 )
279 
280 
281 * The following block introduces the shifts and chases
282 * them down one by one just enough to make space for
283 * the other shifts. The near-the-diagonal block is
284 * of size (ns+1) x ns.
285 
286  CALL zlaset( 'FULL', ns+1, ns+1, czero, cone, qc, ldqc )
287  CALL zlaset( 'FULL', ns, ns, czero, cone, zc, ldzc )
288 
289  DO i = 1, ns
290 * Introduce the shift
291  scale = sqrt( abs( alpha( i ) ) ) * sqrt( abs( beta( i ) ) )
292  IF( scale .GE. safmin .AND. scale .LE. safmax ) THEN
293  alpha( i ) = alpha( i )/scale
294  beta( i ) = beta( i )/scale
295  END IF
296 
297  temp2 = beta( i )*a( ilo, ilo )-alpha( i )*b( ilo, ilo )
298  temp3 = beta( i )*a( ilo+1, ilo )
299 
300  IF ( abs( temp2 ) .GT. safmax .OR.
301  $ abs( temp3 ) .GT. safmax ) THEN
302  temp2 = cone
303  temp3 = czero
304  END IF
305 
306  CALL zlartg( temp2, temp3, c, s, temp )
307  CALL zrot( ns, a( ilo, ilo ), lda, a( ilo+1, ilo ), lda, c,
308  $ s )
309  CALL zrot( ns, b( ilo, ilo ), ldb, b( ilo+1, ilo ), ldb, c,
310  $ s )
311  CALL zrot( ns+1, qc( 1, 1 ), 1, qc( 1, 2 ), 1, c,
312  $ dconjg( s ) )
313 
314 * Chase the shift down
315  DO j = 1, ns-i
316 
317  CALL zlaqz1( .true., .true., j, 1, ns, ihi-ilo+1, a( ilo,
318  $ ilo ), lda, b( ilo, ilo ), ldb, ns+1, 1, qc,
319  $ ldqc, ns, 1, zc, ldzc )
320 
321  END DO
322 
323  END DO
324 
325 * Update the rest of the pencil
326 
327 * Update A(ilo:ilo+ns,ilo+ns:istopm) and B(ilo:ilo+ns,ilo+ns:istopm)
328 * from the left with Qc(1:ns+1,1:ns+1)'
329  sheight = ns+1
330  swidth = istopm-( ilo+ns )+1
331  IF ( swidth > 0 ) THEN
332  CALL zgemm( 'C', 'N', sheight, swidth, sheight, cone, qc, ldqc,
333  $ a( ilo, ilo+ns ), lda, czero, work, sheight )
334  CALL zlacpy( 'ALL', sheight, swidth, work, sheight, a( ilo,
335  $ ilo+ns ), lda )
336  CALL zgemm( 'C', 'N', sheight, swidth, sheight, cone, qc, ldqc,
337  $ b( ilo, ilo+ns ), ldb, czero, work, sheight )
338  CALL zlacpy( 'ALL', sheight, swidth, work, sheight, b( ilo,
339  $ ilo+ns ), ldb )
340  END IF
341  IF ( ilq ) THEN
342  CALL zgemm( 'N', 'N', n, sheight, sheight, cone, q( 1, ilo ),
343  $ ldq, qc, ldqc, czero, work, n )
344  CALL zlacpy( 'ALL', n, sheight, work, n, q( 1, ilo ), ldq )
345  END IF
346 
347 * Update A(istartm:ilo-1,ilo:ilo+ns-1) and B(istartm:ilo-1,ilo:ilo+ns-1)
348 * from the right with Zc(1:ns,1:ns)
349  sheight = ilo-1-istartm+1
350  swidth = ns
351  IF ( sheight > 0 ) THEN
352  CALL zgemm( 'N', 'N', sheight, swidth, swidth, cone,
353  $ a( istartm, ilo ), lda, zc, ldzc, czero, work,
354  $ sheight )
355  CALL zlacpy( 'ALL', sheight, swidth, work, sheight, a( istartm,
356  $ ilo ), lda )
357  CALL zgemm( 'N', 'N', sheight, swidth, swidth, cone,
358  $ b( istartm, ilo ), ldb, zc, ldzc, czero, work,
359  $ sheight )
360  CALL zlacpy( 'ALL', sheight, swidth, work, sheight, b( istartm,
361  $ ilo ), ldb )
362  END IF
363  IF ( ilz ) THEN
364  CALL zgemm( 'N', 'N', n, swidth, swidth, cone, z( 1, ilo ),
365  $ ldz, zc, ldzc, czero, work, n )
366  CALL zlacpy( 'ALL', n, swidth, work, n, z( 1, ilo ), ldz )
367  END IF
368 
369 * The following block chases the shifts down to the bottom
370 * right block. If possible, a shift is moved down npos
371 * positions at a time
372 
373  k = ilo
374  DO WHILE ( k < ihi-ns )
375  np = min( ihi-ns-k, npos )
376 * Size of the near-the-diagonal block
377  nblock = ns+np
378 * istartb points to the first row we will be updating
379  istartb = k+1
380 * istopb points to the last column we will be updating
381  istopb = k+nblock-1
382 
383  CALL zlaset( 'FULL', ns+np, ns+np, czero, cone, qc, ldqc )
384  CALL zlaset( 'FULL', ns+np, ns+np, czero, cone, zc, ldzc )
385 
386 * Near the diagonal shift chase
387  DO i = ns-1, 0, -1
388  DO j = 0, np-1
389 * Move down the block with index k+i+j, updating
390 * the (ns+np x ns+np) block:
391 * (k:k+ns+np,k:k+ns+np-1)
392  CALL zlaqz1( .true., .true., k+i+j, istartb, istopb, ihi,
393  $ a, lda, b, ldb, nblock, k+1, qc, ldqc,
394  $ nblock, k, zc, ldzc )
395  END DO
396  END DO
397 
398 * Update rest of the pencil
399 
400 * Update A(k+1:k+ns+np, k+ns+np:istopm) and
401 * B(k+1:k+ns+np, k+ns+np:istopm)
402 * from the left with Qc(1:ns+np,1:ns+np)'
403  sheight = ns+np
404  swidth = istopm-( k+ns+np )+1
405  IF ( swidth > 0 ) THEN
406  CALL zgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
407  $ ldqc, a( k+1, k+ns+np ), lda, czero, work,
408  $ sheight )
409  CALL zlacpy( 'ALL', sheight, swidth, work, sheight, a( k+1,
410  $ k+ns+np ), lda )
411  CALL zgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
412  $ ldqc, b( k+1, k+ns+np ), ldb, czero, work,
413  $ sheight )
414  CALL zlacpy( 'ALL', sheight, swidth, work, sheight, b( k+1,
415  $ k+ns+np ), ldb )
416  END IF
417  IF ( ilq ) THEN
418  CALL zgemm( 'N', 'N', n, nblock, nblock, cone, q( 1, k+1 ),
419  $ ldq, qc, ldqc, czero, work, n )
420  CALL zlacpy( 'ALL', n, nblock, work, n, q( 1, k+1 ), ldq )
421  END IF
422 
423 * Update A(istartm:k,k:k+ns+npos-1) and B(istartm:k,k:k+ns+npos-1)
424 * from the right with Zc(1:ns+np,1:ns+np)
425  sheight = k-istartm+1
426  swidth = nblock
427  IF ( sheight > 0 ) THEN
428  CALL zgemm( 'N', 'N', sheight, swidth, swidth, cone,
429  $ a( istartm, k ), lda, zc, ldzc, czero, work,
430  $ sheight )
431  CALL zlacpy( 'ALL', sheight, swidth, work, sheight,
432  $ a( istartm, k ), lda )
433  CALL zgemm( 'N', 'N', sheight, swidth, swidth, cone,
434  $ b( istartm, k ), ldb, zc, ldzc, czero, work,
435  $ sheight )
436  CALL zlacpy( 'ALL', sheight, swidth, work, sheight,
437  $ b( istartm, k ), ldb )
438  END IF
439  IF ( ilz ) THEN
440  CALL zgemm( 'N', 'N', n, nblock, nblock, cone, z( 1, k ),
441  $ ldz, zc, ldzc, czero, work, n )
442  CALL zlacpy( 'ALL', n, nblock, work, n, z( 1, k ), ldz )
443  END IF
444 
445  k = k+np
446 
447  END DO
448 
449 * The following block removes the shifts from the bottom right corner
450 * one by one. Updates are initially applied to A(ihi-ns+1:ihi,ihi-ns:ihi).
451 
452  CALL zlaset( 'FULL', ns, ns, czero, cone, qc, ldqc )
453  CALL zlaset( 'FULL', ns+1, ns+1, czero, cone, zc, ldzc )
454 
455 * istartb points to the first row we will be updating
456  istartb = ihi-ns+1
457 * istopb points to the last column we will be updating
458  istopb = ihi
459 
460  DO i = 1, ns
461 * Chase the shift down to the bottom right corner
462  DO ishift = ihi-i, ihi-1
463  CALL zlaqz1( .true., .true., ishift, istartb, istopb, ihi,
464  $ a, lda, b, ldb, ns, ihi-ns+1, qc, ldqc, ns+1,
465  $ ihi-ns, zc, ldzc )
466  END DO
467 
468  END DO
469 
470 * Update rest of the pencil
471 
472 * Update A(ihi-ns+1:ihi, ihi+1:istopm)
473 * from the left with Qc(1:ns,1:ns)'
474  sheight = ns
475  swidth = istopm-( ihi+1 )+1
476  IF ( swidth > 0 ) THEN
477  CALL zgemm( 'C', 'N', sheight, swidth, sheight, cone, qc, ldqc,
478  $ a( ihi-ns+1, ihi+1 ), lda, czero, work, sheight )
479  CALL zlacpy( 'ALL', sheight, swidth, work, sheight,
480  $ a( ihi-ns+1, ihi+1 ), lda )
481  CALL zgemm( 'C', 'N', sheight, swidth, sheight, cone, qc, ldqc,
482  $ b( ihi-ns+1, ihi+1 ), ldb, czero, work, sheight )
483  CALL zlacpy( 'ALL', sheight, swidth, work, sheight,
484  $ b( ihi-ns+1, ihi+1 ), ldb )
485  END IF
486  IF ( ilq ) THEN
487  CALL zgemm( 'N', 'N', n, ns, ns, cone, q( 1, ihi-ns+1 ), ldq,
488  $ qc, ldqc, czero, work, n )
489  CALL zlacpy( 'ALL', n, ns, work, n, q( 1, ihi-ns+1 ), ldq )
490  END IF
491 
492 * Update A(istartm:ihi-ns,ihi-ns:ihi)
493 * from the right with Zc(1:ns+1,1:ns+1)
494  sheight = ihi-ns-istartm+1
495  swidth = ns+1
496  IF ( sheight > 0 ) THEN
497  CALL zgemm( 'N', 'N', sheight, swidth, swidth, cone,
498  $ a( istartm, ihi-ns ), lda, zc, ldzc, czero, work,
499  $ sheight )
500  CALL zlacpy( 'ALL', sheight, swidth, work, sheight, a( istartm,
501  $ ihi-ns ), lda )
502  CALL zgemm( 'N', 'N', sheight, swidth, swidth, cone,
503  $ b( istartm, ihi-ns ), ldb, zc, ldzc, czero, work,
504  $ sheight )
505  CALL zlacpy( 'ALL', sheight, swidth, work, sheight, b( istartm,
506  $ ihi-ns ), ldb )
507  END IF
508  IF ( ilz ) THEN
509  CALL zgemm( 'N', 'N', n, ns+1, ns+1, cone, z( 1, ihi-ns ), ldz,
510  $ zc, ldzc, czero, work, n )
511  CALL zlacpy( 'ALL', n, ns+1, work, n, z( 1, ihi-ns ), ldz )
512  END IF
513 
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:118
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
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
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: