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

◆ slaqz4()

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

SLAQZ4

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

Purpose:
 SLAQZ4 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 REAL array. SR contains
          the real parts of the shifts to use.
[in]SI
          SI is REAL array. SI contains
          the imaginary parts of the shifts to use.
[in]SS
          SS is REAL array. SS contains
          the scale of the shifts to use.
[in,out]A
          A is REAL 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 REAL 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 REAL array, dimension (LDQ, N)
[in]LDQ
          LDQ is INTEGER
[in,out]Z
          Z is REAL array, dimension (LDZ, N)
[in]LDZ
          LDZ is INTEGER
[in,out]QC
          QC is REAL array, dimension (LDQC, NBLOCK_DESIRED)
[in]LDQC
          LDQC is INTEGER
[in,out]ZC
          ZC is REAL array, dimension (LDZC, NBLOCK_DESIRED)
[in]LDZC
          LDZ is INTEGER
[out]WORK
          WORK is REAL 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 210 of file slaqz4.f.

214 IMPLICIT NONE
215
216* Function arguments
217 LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
218 INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
219 $ NSHIFTS, NBLOCK_DESIRED, LDQC, LDZC
220
221 REAL, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
222 $ Z( LDZ, * ), QC( LDQC, * ), ZC( LDZC, * ), WORK( * ), SR( * ),
223 $ SI( * ), SS( * )
224
225 INTEGER, INTENT( OUT ) :: INFO
226
227* Parameters
228 REAL :: ZERO, ONE, HALF
229 parameter( zero = 0.0, one = 1.0, half = 0.5 )
230
231* Local scalars
232 INTEGER :: I, J, NS, ISTARTM, ISTOPM, SHEIGHT, SWIDTH, K, NP,
233 $ ISTARTB, ISTOPB, ISHIFT, NBLOCK, NPOS
234 REAL :: TEMP, V( 3 ), C1, S1, C2, S2, SWAP
235*
236* External functions
237 EXTERNAL :: xerbla, sgemm, slaqz1, slaqz2, slaset, slartg, srot,
238 $ slacpy
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( 'SLAQZ4', -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 slaset( 'FULL', ns+1, ns+1, zero, one, qc, ldqc )
314 CALL slaset( 'FULL', ns, ns, zero, one, zc, ldzc )
315
316 DO i = 1, ns, 2
317* Introduce the shift
318 CALL slaqz1( 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 slartg( temp, v( 3 ), c1, s1, v( 2 ) )
323 CALL slartg( v( 1 ), v( 2 ), c2, s2, temp )
324
325 CALL srot( ns, a( ilo+1, ilo ), lda, a( ilo+2, ilo ), lda, c1,
326 $ s1 )
327 CALL srot( ns, a( ilo, ilo ), lda, a( ilo+1, ilo ), lda, c2,
328 $ s2 )
329 CALL srot( ns, b( ilo+1, ilo ), ldb, b( ilo+2, ilo ), ldb, c1,
330 $ s1 )
331 CALL srot( ns, b( ilo, ilo ), ldb, b( ilo+1, ilo ), ldb, c2,
332 $ s2 )
333 CALL srot( ns+1, qc( 1, 2 ), 1, qc( 1, 3 ), 1, c1, s1 )
334 CALL srot( 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 slaqz2( .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 sgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
355 $ a( ilo, ilo+ns ), lda, zero, work, sheight )
356 CALL slacpy( 'ALL', sheight, swidth, work, sheight, a( ilo,
357 $ ilo+ns ), lda )
358 CALL sgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
359 $ b( ilo, ilo+ns ), ldb, zero, work, sheight )
360 CALL slacpy( 'ALL', sheight, swidth, work, sheight, b( ilo,
361 $ ilo+ns ), ldb )
362 END IF
363 IF ( ilq ) THEN
364 CALL sgemm( 'N', 'N', n, sheight, sheight, one, q( 1, ilo ),
365 $ ldq, qc, ldqc, zero, work, n )
366 CALL slacpy( '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 sgemm( 'N', 'N', sheight, swidth, swidth, one, a( istartm,
375 $ ilo ), lda, zc, ldzc, zero, work, sheight )
376 CALL slacpy( 'ALL', sheight, swidth, work, sheight, a( istartm,
377 $ ilo ), lda )
378 CALL sgemm( 'N', 'N', sheight, swidth, swidth, one, b( istartm,
379 $ ilo ), ldb, zc, ldzc, zero, work, sheight )
380 CALL slacpy( 'ALL', sheight, swidth, work, sheight, b( istartm,
381 $ ilo ), ldb )
382 END IF
383 IF ( ilz ) THEN
384 CALL sgemm( 'N', 'N', n, swidth, swidth, one, z( 1, ilo ), ldz,
385 $ zc, ldzc, zero, work, n )
386 CALL slacpy( '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 slaset( 'FULL', ns+np, ns+np, zero, one, qc, ldqc )
404 CALL slaset( '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 slaqz2( .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 sgemm( 'T', 'N', sheight, swidth, sheight, one, qc,
427 $ ldqc, a( k+1, k+ns+np ), lda, zero, work,
428 $ sheight )
429 CALL slacpy( 'ALL', sheight, swidth, work, sheight, a( k+1,
430 $ k+ns+np ), lda )
431 CALL sgemm( 'T', 'N', sheight, swidth, sheight, one, qc,
432 $ ldqc, b( k+1, k+ns+np ), ldb, zero, work,
433 $ sheight )
434 CALL slacpy( 'ALL', sheight, swidth, work, sheight, b( k+1,
435 $ k+ns+np ), ldb )
436 END IF
437 IF ( ilq ) THEN
438 CALL sgemm( 'N', 'N', n, nblock, nblock, one, q( 1, k+1 ),
439 $ ldq, qc, ldqc, zero, work, n )
440 CALL slacpy( '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 sgemm( 'N', 'N', sheight, swidth, swidth, one,
449 $ a( istartm, k ), lda, zc, ldzc, zero, work,
450 $ sheight )
451 CALL slacpy( 'ALL', sheight, swidth, work, sheight,
452 $ a( istartm, k ), lda )
453 CALL sgemm( 'N', 'N', sheight, swidth, swidth, one,
454 $ b( istartm, k ), ldb, zc, ldzc, zero, work,
455 $ sheight )
456 CALL slacpy( 'ALL', sheight, swidth, work, sheight,
457 $ b( istartm, k ), ldb )
458 END IF
459 IF ( ilz ) THEN
460 CALL sgemm( 'N', 'N', n, nblock, nblock, one, z( 1, k ),
461 $ ldz, zc, ldzc, zero, work, n )
462 CALL slacpy( '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 slaset( 'FULL', ns, ns, zero, one, qc, ldqc )
473 CALL slaset( '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 slaqz2( .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 sgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
498 $ a( ihi-ns+1, ihi+1 ), lda, zero, work, sheight )
499 CALL slacpy( 'ALL', sheight, swidth, work, sheight,
500 $ a( ihi-ns+1, ihi+1 ), lda )
501 CALL sgemm( 'T', 'N', sheight, swidth, sheight, one, qc, ldqc,
502 $ b( ihi-ns+1, ihi+1 ), ldb, zero, work, sheight )
503 CALL slacpy( 'ALL', sheight, swidth, work, sheight,
504 $ b( ihi-ns+1, ihi+1 ), ldb )
505 END IF
506 IF ( ilq ) THEN
507 CALL sgemm( 'N', 'N', n, ns, ns, one, q( 1, ihi-ns+1 ), ldq,
508 $ qc, ldqc, zero, work, n )
509 CALL slacpy( '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 sgemm( 'N', 'N', sheight, swidth, swidth, one, a( istartm,
518 $ ihi-ns ), lda, zc, ldzc, zero, work, sheight )
519 CALL slacpy( 'ALL', sheight, swidth, work, sheight, a( istartm,
520 $ ihi-ns ), lda )
521 CALL sgemm( 'N', 'N', sheight, swidth, swidth, one, b( istartm,
522 $ ihi-ns ), ldb, zc, ldzc, zero, work, sheight )
523 CALL slacpy( 'ALL', sheight, swidth, work, sheight, b( istartm,
524 $ ihi-ns ), ldb )
525 END IF
526 IF ( ilz ) THEN
527 CALL sgemm( 'N', 'N', n, ns+1, ns+1, one, z( 1, ihi-ns ), ldz, zc,
528 $ ldzc, zero, work, n )
529 CALL slacpy( 'ALL', n, ns+1, work, n, z( 1, ihi-ns ), ldz )
530 END IF
531
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
Definition: slartg.f90:111
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine slaqz2(ILQ, ILZ, K, ISTARTM, ISTOPM, IHI, A, LDA, B, LDB, NQ, QSTART, Q, LDQ, NZ, ZSTART, Z, LDZ)
SLAQZ2
Definition: slaqz2.f:173
subroutine slaqz1(A, LDA, B, LDB, SR1, SR2, SI, BETA1, BETA2, V)
SLAQZ1
Definition: slaqz1.f:127
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
Definition: srot.f:92
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
Here is the call graph for this function:
Here is the caller graph for this function: