LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

## ◆ dlaqz4()

 subroutine dlaqz4 ( 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, double precision, dimension( * ), intent(inout) SR, double precision, dimension( * ), intent(inout) SI, double precision, dimension( * ), intent(inout) SS, double precision, dimension( lda, * ), intent(inout) A, integer, intent(in) LDA, double precision, dimension( ldb, * ), intent(inout) B, integer, intent(in) LDB, double precision, dimension( ldq, * ), intent(inout) Q, integer, intent(in) LDQ, double precision, dimension( ldz, * ), intent(inout) Z, integer, intent(in) LDZ, double precision, dimension( ldqc, * ), intent(inout) QC, integer, intent(in) LDQC, double precision, dimension( ldzc, * ), intent(inout) ZC, integer, intent(in) LDZC, double precision, dimension( * ), intent(inout) WORK, integer, intent(in) LWORK, integer, intent(out) INFO )

DLAQZ4

Purpose:
` DLAQZ4 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] SR ``` SR is DOUBLE PRECISION array. SR contains the real parts of the shifts to use.``` [in] SI ``` SI is DOUBLE PRECISION array. SI contains the imaginary parts of the shifts to use.``` [in] SS ``` SS is DOUBLE PRECISION array. SS contains the scale of the shifts to use.``` [in,out] A ` A is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDQ, N)` [in] LDQ ` LDQ is INTEGER` [in,out] Z ` Z is DOUBLE PRECISION array, dimension (LDZ, N)` [in] LDZ ` LDZ is INTEGER` [in,out] QC ` QC is DOUBLE PRECISION array, dimension (LDQC, NBLOCK_DESIRED)` [in] LDQC ` LDQC is INTEGER` [in,out] ZC ` ZC is DOUBLE PRECISION array, dimension (LDZC, NBLOCK_DESIRED)` [in] LDZC ` LDZ is INTEGER` [out] WORK ``` WORK is DOUBLE PRECISION 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```
Date
May 2020

Definition at line 209 of file dlaqz4.f.

213  IMPLICIT NONE
214
215 * Function arguments
216  LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
217  INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
218  \$ NSHIFTS, NBLOCK_DESIRED, LDQC, LDZC
219
220  DOUBLE PRECISION, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ),
221  \$ Q( LDQ, * ), Z( LDZ, * ), QC( LDQC, * ),
222  \$ ZC( LDZC, * ), WORK( * ), SR( * ), SI( * ),
223  \$ SS( * )
224
225  INTEGER, INTENT( OUT ) :: INFO
226
227 * Parameters
228  DOUBLE PRECISION :: ZERO, ONE, HALF
229  parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
230
231 * Local scalars
232  INTEGER :: I, J, NS, ISTARTM, ISTOPM, SHEIGHT, SWIDTH, K, NP,
233  \$ ISTARTB, ISTOPB, ISHIFT, NBLOCK, NPOS
234  DOUBLE PRECISION :: TEMP, V( 3 ), C1, S1, C2, S2, SWAP
235 *
236 * External functions
237  EXTERNAL :: xerbla, dgemm, dlaqz1, dlaqz2, dlaset, dlartg, drot,
238  \$ dlacpy
239
240  info = 0
241  IF ( nblock_desired .LT. nshifts+1 ) THEN
242  info = -8
243  END IF
244  IF ( lwork .EQ.-1 ) THEN
245 * workspace query, quick return
246  work( 1 ) = n*nblock_desired
247  RETURN
248  ELSE IF ( lwork .LT. n*nblock_desired ) THEN
249  info = -25
250  END IF
251
252  IF( info.NE.0 ) THEN
253  CALL xerbla( 'DLAQZ4', -info )
254  RETURN
255  END IF
256
257 * Executable statements
258
259  IF ( nshifts .LT. 2 ) THEN
260  RETURN
261  END IF
262
263  IF ( ilo .GE. ihi ) THEN
264  RETURN
265  END IF
266
267  IF ( ilschur ) THEN
268  istartm = 1
269  istopm = n
270  ELSE
271  istartm = ilo
272  istopm = ihi
273  END IF
274
275 * Shuffle shifts into pairs of real shifts and pairs
276 * of complex conjugate shifts assuming complex
278 * another
279
280  DO i = 1, nshifts-2, 2
281  IF( si( i ).NE.-si( i+1 ) ) THEN
282 *
283  swap = sr( i )
284  sr( i ) = sr( i+1 )
285  sr( i+1 ) = sr( i+2 )
286  sr( i+2 ) = swap
287
288  swap = si( i )
289  si( i ) = si( i+1 )
290  si( i+1 ) = si( i+2 )
291  si( i+2 ) = swap
292
293  swap = ss( i )
294  ss( i ) = ss( i+1 )
295  ss( i+1 ) = ss( i+2 )
296  ss( i+2 ) = swap
297  END IF
298  END DO
299
300 * NSHFTS is supposed to be even, but if it is odd,
301 * then simply reduce it by one. The shuffle above
302 * ensures that the dropped shift is real and that
303 * the remaining shifts are paired.
304
305  ns = nshifts-mod( nshifts, 2 )
306  npos = max( nblock_desired-ns, 1 )
307
308 * The following block introduces the shifts and chases
309 * them down one by one just enough to make space for
310 * the other shifts. The near-the-diagonal block is
311 * of size (ns+1) x ns.
312
313  CALL dlaset( 'FULL', ns+1, ns+1, zero, one, qc, ldqc )
314  CALL dlaset( 'FULL', ns, ns, zero, one, zc, ldzc )
315
316  DO i = 1, ns, 2
317 * Introduce the shift
318  CALL dlaqz1( a( ilo, ilo ), lda, b( ilo, ilo ), ldb, sr( i ),
319  \$ sr( i+1 ), si( i ), ss( i ), ss( i+1 ), v )
320
321  temp = v( 2 )
322  CALL dlartg( temp, v( 3 ), c1, s1, v( 2 ) )
323  CALL dlartg( v( 1 ), v( 2 ), c2, s2, temp )
324
325  CALL drot( ns, a( ilo+1, ilo ), lda, a( ilo+2, ilo ), lda, c1,
326  \$ s1 )
327  CALL drot( ns, a( ilo, ilo ), lda, a( ilo+1, ilo ), lda, c2,
328  \$ s2 )
329  CALL drot( ns, b( ilo+1, ilo ), ldb, b( ilo+2, ilo ), ldb, c1,
330  \$ s1 )
331  CALL drot( ns, b( ilo, ilo ), ldb, b( ilo+1, ilo ), ldb, c2,
332  \$ s2 )
333  CALL drot( ns+1, qc( 1, 2 ), 1, qc( 1, 3 ), 1, c1, s1 )
334  CALL drot( ns+1, qc( 1, 1 ), 1, qc( 1, 2 ), 1, c2, s2 )
335
336 * Chase the shift down
337  DO j = 1, ns-1-i
338
339  CALL dlaqz2( .true., .true., j, 1, ns, ihi-ilo+1, a( ilo,
340  \$ ilo ), lda, b( ilo, ilo ), ldb, ns+1, 1, qc,
341  \$ ldqc, ns, 1, zc, ldzc )
342
343  END DO
344
345  END DO
346
347 * Update the rest of the pencil
348
349 * Update A(ilo:ilo+ns,ilo+ns:istopm) and B(ilo:ilo+ns,ilo+ns:istopm)
350 * from the left with Qc(1:ns+1,1:ns+1)'
351  sheight = ns+1
352  swidth = istopm-( ilo+ns )+1
353  IF ( swidth > 0 ) THEN
354  CALL dgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
355  \$ a( ilo, ilo+ns ), lda, zero, work, sheight )
356  CALL dlacpy( 'ALL', sheight, swidth, work, sheight, a( ilo,
357  \$ ilo+ns ), lda )
358  CALL dgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
359  \$ b( ilo, ilo+ns ), ldb, zero, work, sheight )
360  CALL dlacpy( 'ALL', sheight, swidth, work, sheight, b( ilo,
361  \$ ilo+ns ), ldb )
362  END IF
363  IF ( ilq ) THEN
364  CALL dgemm( 'N', 'N', n, sheight, sheight, one, q( 1, ilo ),
365  \$ ldq, qc, ldqc, zero, work, n )
366  CALL dlacpy( 'ALL', n, sheight, work, n, q( 1, ilo ), ldq )
367  END IF
368
369 * Update A(istartm:ilo-1,ilo:ilo+ns-1) and B(istartm:ilo-1,ilo:ilo+ns-1)
370 * from the right with Zc(1:ns,1:ns)
371  sheight = ilo-1-istartm+1
372  swidth = ns
373  IF ( sheight > 0 ) THEN
374  CALL dgemm( 'N', 'N', sheight, swidth, swidth, one, a( istartm,
375  \$ ilo ), lda, zc, ldzc, zero, work, sheight )
376  CALL dlacpy( 'ALL', sheight, swidth, work, sheight, a( istartm,
377  \$ ilo ), lda )
378  CALL dgemm( 'N', 'N', sheight, swidth, swidth, one, b( istartm,
379  \$ ilo ), ldb, zc, ldzc, zero, work, sheight )
380  CALL dlacpy( 'ALL', sheight, swidth, work, sheight, b( istartm,
381  \$ ilo ), ldb )
382  END IF
383  IF ( ilz ) THEN
384  CALL dgemm( 'N', 'N', n, swidth, swidth, one, z( 1, ilo ), ldz,
385  \$ zc, ldzc, zero, work, n )
386  CALL dlacpy( 'ALL', n, swidth, work, n, z( 1, ilo ), ldz )
387  END IF
388
389 * The following block chases the shifts down to the bottom
390 * right block. If possible, a shift is moved down npos
391 * positions at a time
392
393  k = ilo
394  DO WHILE ( k < ihi-ns )
395  np = min( ihi-ns-k, npos )
396 * Size of the near-the-diagonal block
397  nblock = ns+np
398 * istartb points to the first row we will be updating
399  istartb = k+1
400 * istopb points to the last column we will be updating
401  istopb = k+nblock-1
402
403  CALL dlaset( 'FULL', ns+np, ns+np, zero, one, qc, ldqc )
404  CALL dlaset( 'FULL', ns+np, ns+np, zero, one, zc, ldzc )
405
406 * Near the diagonal shift chase
407  DO i = ns-1, 0, -2
408  DO j = 0, np-1
409 * Move down the block with index k+i+j-1, updating
410 * the (ns+np x ns+np) block:
411 * (k:k+ns+np,k:k+ns+np-1)
412  CALL dlaqz2( .true., .true., k+i+j-1, istartb, istopb,
413  \$ ihi, a, lda, b, ldb, nblock, k+1, qc, ldqc,
414  \$ nblock, k, zc, ldzc )
415  END DO
416  END DO
417
418 * Update rest of the pencil
419
420 * Update A(k+1:k+ns+np, k+ns+np:istopm) and
421 * B(k+1:k+ns+np, k+ns+np:istopm)
422 * from the left with Qc(1:ns+np,1:ns+np)'
423  sheight = ns+np
424  swidth = istopm-( k+ns+np )+1
425  IF ( swidth > 0 ) THEN
426  CALL dgemm( 'T', 'N', sheight, swidth, sheight, one, qc,
427  \$ ldqc, a( k+1, k+ns+np ), lda, zero, work,
428  \$ sheight )
429  CALL dlacpy( 'ALL', sheight, swidth, work, sheight, a( k+1,
430  \$ k+ns+np ), lda )
431  CALL dgemm( 'T', 'N', sheight, swidth, sheight, one, qc,
432  \$ ldqc, b( k+1, k+ns+np ), ldb, zero, work,
433  \$ sheight )
434  CALL dlacpy( 'ALL', sheight, swidth, work, sheight, b( k+1,
435  \$ k+ns+np ), ldb )
436  END IF
437  IF ( ilq ) THEN
438  CALL dgemm( 'N', 'N', n, nblock, nblock, one, q( 1, k+1 ),
439  \$ ldq, qc, ldqc, zero, work, n )
440  CALL dlacpy( 'ALL', n, nblock, work, n, q( 1, k+1 ), ldq )
441  END IF
442
443 * Update A(istartm:k,k:k+ns+npos-1) and B(istartm:k,k:k+ns+npos-1)
444 * from the right with Zc(1:ns+np,1:ns+np)
445  sheight = k-istartm+1
446  swidth = nblock
447  IF ( sheight > 0 ) THEN
448  CALL dgemm( 'N', 'N', sheight, swidth, swidth, one,
449  \$ a( istartm, k ), lda, zc, ldzc, zero, work,
450  \$ sheight )
451  CALL dlacpy( 'ALL', sheight, swidth, work, sheight,
452  \$ a( istartm, k ), lda )
453  CALL dgemm( 'N', 'N', sheight, swidth, swidth, one,
454  \$ b( istartm, k ), ldb, zc, ldzc, zero, work,
455  \$ sheight )
456  CALL dlacpy( 'ALL', sheight, swidth, work, sheight,
457  \$ b( istartm, k ), ldb )
458  END IF
459  IF ( ilz ) THEN
460  CALL dgemm( 'N', 'N', n, nblock, nblock, one, z( 1, k ),
461  \$ ldz, zc, ldzc, zero, work, n )
462  CALL dlacpy( 'ALL', n, nblock, work, n, z( 1, k ), ldz )
463  END IF
464
465  k = k+np
466
467  END DO
468
469 * The following block removes the shifts from the bottom right corner
470 * one by one. Updates are initially applied to A(ihi-ns+1:ihi,ihi-ns:ihi).
471
472  CALL dlaset( 'FULL', ns, ns, zero, one, qc, ldqc )
473  CALL dlaset( 'FULL', ns+1, ns+1, zero, one, zc, ldzc )
474
475 * istartb points to the first row we will be updating
476  istartb = ihi-ns+1
477 * istopb points to the last column we will be updating
478  istopb = ihi
479
480  DO i = 1, ns, 2
481 * Chase the shift down to the bottom right corner
482  DO ishift = ihi-i-1, ihi-2
483  CALL dlaqz2( .true., .true., ishift, istartb, istopb, ihi,
484  \$ a, lda, b, ldb, ns, ihi-ns+1, qc, ldqc, ns+1,
485  \$ ihi-ns, zc, ldzc )
486  END DO
487
488  END DO
489
490 * Update rest of the pencil
491
492 * Update A(ihi-ns+1:ihi, ihi+1:istopm)
493 * from the left with Qc(1:ns,1:ns)'
494  sheight = ns
495  swidth = istopm-( ihi+1 )+1
496  IF ( swidth > 0 ) THEN
497  CALL dgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
498  \$ a( ihi-ns+1, ihi+1 ), lda, zero, work, sheight )
499  CALL dlacpy( 'ALL', sheight, swidth, work, sheight,
500  \$ a( ihi-ns+1, ihi+1 ), lda )
501  CALL dgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
502  \$ b( ihi-ns+1, ihi+1 ), ldb, zero, work, sheight )
503  CALL dlacpy( 'ALL', sheight, swidth, work, sheight,
504  \$ b( ihi-ns+1, ihi+1 ), ldb )
505  END IF
506  IF ( ilq ) THEN
507  CALL dgemm( 'N', 'N', n, ns, ns, one, q( 1, ihi-ns+1 ), ldq,
508  \$ qc, ldqc, zero, work, n )
509  CALL dlacpy( 'ALL', n, ns, work, n, q( 1, ihi-ns+1 ), ldq )
510  END IF
511
512 * Update A(istartm:ihi-ns,ihi-ns:ihi)
513 * from the right with Zc(1:ns+1,1:ns+1)
514  sheight = ihi-ns-istartm+1
515  swidth = ns+1
516  IF ( sheight > 0 ) THEN
517  CALL dgemm( 'N', 'N', sheight, swidth, swidth, one, a( istartm,
518  \$ ihi-ns ), lda, zc, ldzc, zero, work, sheight )
519  CALL dlacpy( 'ALL', sheight, swidth, work, sheight, a( istartm,
520  \$ ihi-ns ), lda )
521  CALL dgemm( 'N', 'N', sheight, swidth, swidth, one, b( istartm,
522  \$ ihi-ns ), ldb, zc, ldzc, zero, work, sheight )
523  CALL dlacpy( 'ALL', sheight, swidth, work, sheight, b( istartm,
524  \$ ihi-ns ), ldb )
525  END IF
526  IF ( ilz ) THEN
527  CALL dgemm( 'N', 'N', n, ns+1, ns+1, one, z( 1, ihi-ns ), ldz,
528  \$ zc, ldzc, zero, work, n )
529  CALL dlacpy( 'ALL', n, ns+1, work, n, z( 1, ihi-ns ), ldz )
530  END IF
531
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f90:113
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:92
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
subroutine dlaqz1(A, LDA, B, LDB, SR1, SR2, SI, BETA1, BETA2, V)
DLAQZ1
Definition: dlaqz1.f:127
subroutine dlaqz2(ILQ, ILZ, K, ISTARTM, ISTOPM, IHI, A, LDA, B, LDB, NQ, QSTART, Q, LDQ, NZ, ZSTART, Z, LDZ)
DLAQZ2
Definition: dlaqz2.f:174
Here is the call graph for this function:
Here is the caller graph for this function: