LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ claqz3()

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

CLAQZ3

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

Purpose:
 CLAQZ3 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 array. SR contains
          the alpha parts of the shifts to use.
[in]BETA
          BETA is COMPLEX array. SS contains
          the scale of the shifts to use.
[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
[in,out]QC
          QC is COMPLEX array, dimension (LDQC, NBLOCK_DESIRED)
[in]LDQC
          LDQC is INTEGER
[in,out]ZC
          ZC is COMPLEX array, dimension (LDZC, NBLOCK_DESIRED)
[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]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 203 of file claqz3.f.

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