207 SUBROUTINE dlaqz4( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS,
208 $ NBLOCK_DESIRED, SR, SI, SS, A, LDA, B, LDB, Q,
209 $ LDQ, Z, LDZ, QC, LDQC, ZC, LDZC, WORK, LWORK,
214 LOGICAL,
INTENT( IN ) :: ILSCHUR, ILQ, ILZ
215 INTEGER,
INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
216 $ NSHIFTS, NBLOCK_DESIRED, LDQC, LDZC
218 DOUBLE PRECISION,
INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ),
219 $ Q( LDQ, * ), Z( LDZ, * ), QC( LDQC, * ),
220 $ ZC( LDZC, * ), WORK( * ), SR( * ), SI( * ),
223 INTEGER,
INTENT( OUT ) :: INFO
226 DOUBLE PRECISION :: ZERO, ONE, HALF
227 PARAMETER( ZERO = 0.0d0, one = 1.0d0, half = 0.5d0 )
230 INTEGER :: I, J, NS, ISTARTM, ISTOPM, SHEIGHT, SWIDTH, K, NP,
231 $ ISTARTB, ISTOPB, ISHIFT, NBLOCK, NPOS
232 DOUBLE PRECISION :: TEMP, V( 3 ), C1, S1, C2, S2, SWAP
239 IF ( nblock_desired .LT. nshifts+1 )
THEN
242 IF ( lwork .EQ.-1 )
THEN
244 work( 1 ) = n*nblock_desired
246 ELSE IF ( lwork .LT. n*nblock_desired )
THEN
251 CALL xerbla(
'DLAQZ4', -info )
257 IF ( nshifts .LT. 2 )
THEN
261 IF ( ilo .GE. ihi )
THEN
278 DO i = 1, nshifts-2, 2
279 IF( si( i ).NE.-si( i+1 ) )
THEN
283 sr( i+1 ) = sr( i+2 )
288 si( i+1 ) = si( i+2 )
293 ss( i+1 ) = ss( i+2 )
303 ns = nshifts-mod( nshifts, 2 )
304 npos = max( nblock_desired-ns, 1 )
311 CALL dlaset(
'FULL', ns+1, ns+1, zero, one, qc, ldqc )
312 CALL dlaset(
'FULL', ns, ns, zero, one, zc, ldzc )
316 CALL dlaqz1( a( ilo, ilo ), lda, b( ilo, ilo ), ldb,
318 $ sr( i+1 ), si( i ), ss( i ), ss( i+1 ), v )
321 CALL dlartg( temp, v( 3 ), c1, s1, v( 2 ) )
322 CALL dlartg( v( 1 ), v( 2 ), c2, s2, temp )
324 CALL drot( ns, a( ilo+1, ilo ), lda, a( ilo+2, ilo ), lda,
327 CALL drot( ns, a( ilo, ilo ), lda, a( ilo+1, ilo ), lda, c2,
329 CALL drot( ns, b( ilo+1, ilo ), ldb, b( ilo+2, ilo ), ldb,
332 CALL drot( ns, b( ilo, ilo ), ldb, b( ilo+1, ilo ), ldb, c2,
334 CALL drot( ns+1, qc( 1, 2 ), 1, qc( 1, 3 ), 1, c1, s1 )
335 CALL drot( ns+1, qc( 1, 1 ), 1, qc( 1, 2 ), 1, c2, s2 )
340 CALL dlaqz2( .true., .true., j, 1, ns, ihi-ilo+1, a( ilo,
341 $ ilo ), lda, b( ilo, ilo ), ldb, ns+1, 1, qc,
342 $ ldqc, ns, 1, zc, ldzc )
353 swidth = istopm-( ilo+ns )+1
354 IF ( swidth > 0 )
THEN
355 CALL dgemm(
'T',
'N', sheight, swidth, sheight, one, qc,
357 $ a( ilo, ilo+ns ), lda, zero, work, sheight )
358 CALL dlacpy(
'ALL', sheight, swidth, work, sheight, a( ilo,
360 CALL dgemm(
'T',
'N', sheight, swidth, sheight, one, qc,
362 $ b( ilo, ilo+ns ), ldb, zero, work, sheight )
363 CALL dlacpy(
'ALL', sheight, swidth, work, sheight, b( ilo,
367 CALL dgemm(
'N',
'N', n, sheight, sheight, one, q( 1, ilo ),
368 $ ldq, qc, ldqc, zero, work, n )
369 CALL dlacpy(
'ALL', n, sheight, work, n, q( 1, ilo ), ldq )
374 sheight = ilo-1-istartm+1
376 IF ( sheight > 0 )
THEN
377 CALL dgemm(
'N',
'N', sheight, swidth, swidth, one,
379 $ ilo ), lda, zc, ldzc, zero, work, sheight )
380 CALL dlacpy(
'ALL', sheight, swidth, work, sheight,
383 CALL dgemm(
'N',
'N', sheight, swidth, swidth, one,
385 $ ilo ), ldb, zc, ldzc, zero, work, sheight )
386 CALL dlacpy(
'ALL', sheight, swidth, work, sheight,
391 CALL dgemm(
'N',
'N', n, swidth, swidth, one, z( 1, ilo ),
393 $ zc, ldzc, zero, work, n )
394 CALL dlacpy(
'ALL', n, swidth, work, n, z( 1, ilo ), ldz )
402 DO WHILE ( k < ihi-ns )
403 np = min( ihi-ns-k, npos )
411 CALL dlaset(
'FULL', ns+np, ns+np, zero, one, qc, ldqc )
412 CALL dlaset(
'FULL', ns+np, ns+np, zero, one, zc, ldzc )
420 CALL dlaqz2( .true., .true., k+i+j-1, istartb, istopb,
421 $ ihi, a, lda, b, ldb, nblock, k+1, qc, ldqc,
422 $ nblock, k, zc, ldzc )
432 swidth = istopm-( k+ns+np )+1
433 IF ( swidth > 0 )
THEN
434 CALL dgemm(
'T',
'N', sheight, swidth, sheight, one, qc,
435 $ ldqc, a( k+1, k+ns+np ), lda, zero, work,
437 CALL dlacpy(
'ALL', sheight, swidth, work, sheight,
440 CALL dgemm(
'T',
'N', sheight, swidth, sheight, one, qc,
441 $ ldqc, b( k+1, k+ns+np ), ldb, zero, work,
443 CALL dlacpy(
'ALL', sheight, swidth, work, sheight,
448 CALL dgemm(
'N',
'N', n, nblock, nblock, one, q( 1,
450 $ ldq, qc, ldqc, zero, work, n )
451 CALL dlacpy(
'ALL', n, nblock, work, n, q( 1, k+1 ),
457 sheight = k-istartm+1
459 IF ( sheight > 0 )
THEN
460 CALL dgemm(
'N',
'N', sheight, swidth, swidth, one,
461 $ a( istartm, k ), lda, zc, ldzc, zero, work,
463 CALL dlacpy(
'ALL', sheight, swidth, work, sheight,
464 $ a( istartm, k ), lda )
465 CALL dgemm(
'N',
'N', sheight, swidth, swidth, one,
466 $ b( istartm, k ), ldb, zc, ldzc, zero, work,
468 CALL dlacpy(
'ALL', sheight, swidth, work, sheight,
469 $ b( istartm, k ), ldb )
472 CALL dgemm(
'N',
'N', n, nblock, nblock, one, z( 1, k ),
473 $ ldz, zc, ldzc, zero, work, n )
474 CALL dlacpy(
'ALL', n, nblock, work, n, z( 1, k ), ldz )
484 CALL dlaset(
'FULL', ns, ns, zero, one, qc, ldqc )
485 CALL dlaset(
'FULL', ns+1, ns+1, zero, one, zc, ldzc )
494 DO ishift = ihi-i-1, ihi-2
495 CALL dlaqz2( .true., .true., ishift, istartb, istopb,
497 $ a, lda, b, ldb, ns, ihi-ns+1, qc, ldqc, ns+1,
508 swidth = istopm-( ihi+1 )+1
509 IF ( swidth > 0 )
THEN
510 CALL dgemm(
'T',
'N', sheight, swidth, sheight, one, qc,
512 $ a( ihi-ns+1, ihi+1 ), lda, zero, work, sheight )
513 CALL dlacpy(
'ALL', sheight, swidth, work, sheight,
514 $ a( ihi-ns+1, ihi+1 ), lda )
515 CALL dgemm(
'T',
'N', sheight, swidth, sheight, one, qc,
517 $ b( ihi-ns+1, ihi+1 ), ldb, zero, work, sheight )
518 CALL dlacpy(
'ALL', sheight, swidth, work, sheight,
519 $ b( ihi-ns+1, ihi+1 ), ldb )
522 CALL dgemm(
'N',
'N', n, ns, ns, one, q( 1, ihi-ns+1 ), ldq,
523 $ qc, ldqc, zero, work, n )
524 CALL dlacpy(
'ALL', n, ns, work, n, q( 1, ihi-ns+1 ), ldq )
529 sheight = ihi-ns-istartm+1
531 IF ( sheight > 0 )
THEN
532 CALL dgemm(
'N',
'N', sheight, swidth, swidth, one,
534 $ ihi-ns ), lda, zc, ldzc, zero, work, sheight )
535 CALL dlacpy(
'ALL', sheight, swidth, work, sheight,
538 CALL dgemm(
'N',
'N', sheight, swidth, swidth, one,
540 $ ihi-ns ), ldb, zc, ldzc, zero, work, sheight )
541 CALL dlacpy(
'ALL', sheight, swidth, work, sheight,
546 CALL dgemm(
'N',
'N', n, ns+1, ns+1, one, z( 1, ihi-ns ),
548 $ zc, ldzc, zero, work, n )
549 CALL dlacpy(
'ALL', n, ns+1, work, n, z( 1, ihi-ns ), ldz )