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

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```
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 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: