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

◆ 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

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

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
Author
Thijs Steel, KU Leuven
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
277* conjugate shifts are already adjacent to one
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:111
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: