211 LOGICAL,
INTENT( IN ) :: ILSCHUR, ILQ, ILZ
212 INTEGER,
INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
213 $ NSHIFTS, NBLOCK_DESIRED, LDQC, LDZC
215 COMPLEX*16,
INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ,
216 $ * ), Z( LDZ, * ), QC( LDQC, * ), ZC( LDZC, * ), WORK( * ),
217 $ ALPHA( * ), BETA( * )
219 INTEGER,
INTENT( OUT ) :: INFO
222 COMPLEX*16 CZERO, CONE
223 parameter( czero = ( 0.0d+0, 0.0d+0 ), cone = ( 1.0d+0,
225 DOUBLE PRECISION :: ZERO, ONE, HALF
226 parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
229 INTEGER :: I, J, NS, ISTARTM, ISTOPM, SHEIGHT, SWIDTH, K, NP,
230 $ ISTARTB, ISTOPB, ISHIFT, NBLOCK, NPOS
231 DOUBLE PRECISION :: SAFMIN, SAFMAX, C, SCALE
232 COMPLEX*16 :: TEMP, TEMP2, TEMP3, S
237 DOUBLE PRECISION,
EXTERNAL :: DLAMCH
240 IF ( nblock_desired .LT. nshifts+1 )
THEN
243 IF ( lwork .EQ.-1 )
THEN
245 work( 1 ) = n*nblock_desired
247 ELSE IF ( lwork .LT. n*nblock_desired )
THEN
252 CALL xerbla(
'ZLAQZ3', -info )
261 safmin =
dlamch(
'SAFE MINIMUM' )
263 CALL dlabad( safmin, safmax )
265 IF ( ilo .GE. ihi )
THEN
278 npos = max( nblock_desired-ns, 1 )
286 CALL zlaset(
'FULL', ns+1, ns+1, czero, cone, qc, ldqc )
287 CALL zlaset(
'FULL', ns, ns, czero, cone, zc, ldzc )
291 scale = sqrt( abs( alpha( i ) ) ) * sqrt( abs( beta( i ) ) )
292 IF( scale .GE. safmin .AND. scale .LE. safmax )
THEN
293 alpha( i ) = alpha( i )/scale
294 beta( i ) = beta( i )/scale
297 temp2 = beta( i )*a( ilo, ilo )-alpha( i )*b( ilo, ilo )
298 temp3 = beta( i )*a( ilo+1, ilo )
300 IF ( abs( temp2 ) .GT. safmax .OR.
301 $ abs( temp3 ) .GT. safmax )
THEN
306 CALL zlartg( temp2, temp3, c, s, temp )
307 CALL zrot( ns, a( ilo, ilo ), lda, a( ilo+1, ilo ), lda, c,
309 CALL zrot( ns, b( ilo, ilo ), ldb, b( ilo+1, ilo ), ldb, c,
311 CALL zrot( ns+1, qc( 1, 1 ), 1, qc( 1, 2 ), 1, c,
317 CALL zlaqz1( .true., .true., j, 1, ns, ihi-ilo+1, a( ilo,
318 $ ilo ), lda, b( ilo, ilo ), ldb, ns+1, 1, qc,
319 $ ldqc, ns, 1, zc, ldzc )
330 swidth = istopm-( ilo+ns )+1
331 IF ( swidth > 0 )
THEN
332 CALL zgemm(
'C',
'N', sheight, swidth, sheight, cone, qc, ldqc,
333 $ a( ilo, ilo+ns ), lda, czero, work, sheight )
334 CALL zlacpy(
'ALL', sheight, swidth, work, sheight, a( ilo,
336 CALL zgemm(
'C',
'N', sheight, swidth, sheight, cone, qc, ldqc,
337 $ b( ilo, ilo+ns ), ldb, czero, work, sheight )
338 CALL zlacpy(
'ALL', sheight, swidth, work, sheight, b( ilo,
342 CALL zgemm(
'N',
'N', n, sheight, sheight, cone, q( 1, ilo ),
343 $ ldq, qc, ldqc, czero, work, n )
344 CALL zlacpy(
'ALL', n, sheight, work, n, q( 1, ilo ), ldq )
349 sheight = ilo-1-istartm+1
351 IF ( sheight > 0 )
THEN
352 CALL zgemm(
'N',
'N', sheight, swidth, swidth, cone,
353 $ a( istartm, ilo ), lda, zc, ldzc, czero, work,
355 CALL zlacpy(
'ALL', sheight, swidth, work, sheight, a( istartm,
357 CALL zgemm(
'N',
'N', sheight, swidth, swidth, cone,
358 $ b( istartm, ilo ), ldb, zc, ldzc, czero, work,
360 CALL zlacpy(
'ALL', sheight, swidth, work, sheight, b( istartm,
364 CALL zgemm(
'N',
'N', n, swidth, swidth, cone, z( 1, ilo ),
365 $ ldz, zc, ldzc, czero, work, n )
366 CALL zlacpy(
'ALL', n, swidth, work, n, z( 1, ilo ), ldz )
374 DO WHILE ( k < ihi-ns )
375 np = min( ihi-ns-k, npos )
383 CALL zlaset(
'FULL', ns+np, ns+np, czero, cone, qc, ldqc )
384 CALL zlaset(
'FULL', ns+np, ns+np, czero, cone, zc, ldzc )
392 CALL zlaqz1( .true., .true., k+i+j, istartb, istopb, ihi,
393 $ a, lda, b, ldb, nblock, k+1, qc, ldqc,
394 $ nblock, k, zc, ldzc )
404 swidth = istopm-( k+ns+np )+1
405 IF ( swidth > 0 )
THEN
406 CALL zgemm(
'C',
'N', sheight, swidth, sheight, cone, qc,
407 $ ldqc, a( k+1, k+ns+np ), lda, czero, work,
409 CALL zlacpy(
'ALL', sheight, swidth, work, sheight, a( k+1,
411 CALL zgemm(
'C',
'N', sheight, swidth, sheight, cone, qc,
412 $ ldqc, b( k+1, k+ns+np ), ldb, czero, work,
414 CALL zlacpy(
'ALL', sheight, swidth, work, sheight, b( k+1,
418 CALL zgemm(
'N',
'N', n, nblock, nblock, cone, q( 1, k+1 ),
419 $ ldq, qc, ldqc, czero, work, n )
420 CALL zlacpy(
'ALL', n, nblock, work, n, q( 1, k+1 ), ldq )
425 sheight = k-istartm+1
427 IF ( sheight > 0 )
THEN
428 CALL zgemm(
'N',
'N', sheight, swidth, swidth, cone,
429 $ a( istartm, k ), lda, zc, ldzc, czero, work,
431 CALL zlacpy(
'ALL', sheight, swidth, work, sheight,
432 $ a( istartm, k ), lda )
433 CALL zgemm(
'N',
'N', sheight, swidth, swidth, cone,
434 $ b( istartm, k ), ldb, zc, ldzc, czero, work,
436 CALL zlacpy(
'ALL', sheight, swidth, work, sheight,
437 $ b( istartm, k ), ldb )
440 CALL zgemm(
'N',
'N', n, nblock, nblock, cone, z( 1, k ),
441 $ ldz, zc, ldzc, czero, work, n )
442 CALL zlacpy(
'ALL', n, nblock, work, n, z( 1, k ), ldz )
452 CALL zlaset(
'FULL', ns, ns, czero, cone, qc, ldqc )
453 CALL zlaset(
'FULL', ns+1, ns+1, czero, cone, zc, ldzc )
462 DO ishift = ihi-i, ihi-1
463 CALL zlaqz1( .true., .true., ishift, istartb, istopb, ihi,
464 $ a, lda, b, ldb, ns, ihi-ns+1, qc, ldqc, ns+1,
475 swidth = istopm-( ihi+1 )+1
476 IF ( swidth > 0 )
THEN
477 CALL zgemm(
'C',
'N', sheight, swidth, sheight, cone, qc, ldqc,
478 $ a( ihi-ns+1, ihi+1 ), lda, czero, work, sheight )
479 CALL zlacpy(
'ALL', sheight, swidth, work, sheight,
480 $ a( ihi-ns+1, ihi+1 ), lda )
481 CALL zgemm(
'C',
'N', sheight, swidth, sheight, cone, qc, ldqc,
482 $ b( ihi-ns+1, ihi+1 ), ldb, czero, work, sheight )
483 CALL zlacpy(
'ALL', sheight, swidth, work, sheight,
484 $ b( ihi-ns+1, ihi+1 ), ldb )
487 CALL zgemm(
'N',
'N', n, ns, ns, cone, q( 1, ihi-ns+1 ), ldq,
488 $ qc, ldqc, czero, work, n )
489 CALL zlacpy(
'ALL', n, ns, work, n, q( 1, ihi-ns+1 ), ldq )
494 sheight = ihi-ns-istartm+1
496 IF ( sheight > 0 )
THEN
497 CALL zgemm(
'N',
'N', sheight, swidth, swidth, cone,
498 $ a( istartm, ihi-ns ), lda, zc, ldzc, czero, work,
500 CALL zlacpy(
'ALL', sheight, swidth, work, sheight, a( istartm,
502 CALL zgemm(
'N',
'N', sheight, swidth, swidth, cone,
503 $ b( istartm, ihi-ns ), ldb, zc, ldzc, czero, work,
505 CALL zlacpy(
'ALL', sheight, swidth, work, sheight, b( istartm,
509 CALL zgemm(
'N',
'N', n, ns+1, ns+1, cone, z( 1, ihi-ns ), ldz,
510 $ zc, ldzc, czero, work, n )
511 CALL zlacpy(
'ALL', n, ns+1, work, n, z( 1, ihi-ns ), ldz )
double precision function dlamch(CMACH)
DLAMCH
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zlaqz1(ILQ, ILZ, K, ISTARTM, ISTOPM, IHI, A, LDA, B, LDB, NQ, QSTART, Q, LDQ, NZ, ZSTART, Z, LDZ)
ZLAQZ1
subroutine zrot(N, CX, INCX, CY, INCY, C, S)
ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.