210 SUBROUTINE slaqz4( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS,
211 $ NBLOCK_DESIRED, SR, SI, SS, A, LDA, B, LDB, Q,
212 $ LDQ, Z, LDZ, QC, LDQC, ZC, LDZC, WORK, LWORK,
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
221 REAL,
INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
222 $ Z( LDZ, * ), QC( LDQC, * ), ZC( LDZC, * ), WORK( * ), SR( * ),
225 INTEGER,
INTENT( OUT ) :: INFO
228 REAL :: ZERO, ONE, HALF
229 PARAMETER( ZERO = 0.0, one = 1.0, half = 0.5 )
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
241 IF ( nblock_desired .LT. nshifts+1 )
THEN
244 IF ( lwork .EQ.-1 )
THEN
246 work( 1 ) = n*nblock_desired
248 ELSE IF ( lwork .LT. n*nblock_desired )
THEN
253 CALL xerbla(
'SLAQZ4', -info )
259 IF ( nshifts .LT. 2 )
THEN
263 IF ( ilo .GE. ihi )
THEN
280 DO i = 1, nshifts-2, 2
281 IF( si( i ).NE.-si( i+1 ) )
THEN
285 sr( i+1 ) = sr( i+2 )
290 si( i+1 ) = si( i+2 )
295 ss( i+1 ) = ss( i+2 )
305 ns = nshifts-mod( nshifts, 2 )
306 npos = max( nblock_desired-ns, 1 )
313 CALL slaset(
'FULL', ns+1, ns+1, zero, one, qc, ldqc )
314 CALL slaset(
'FULL', ns, ns, zero, one, zc, ldzc )
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 )
322 CALL slartg( temp, v( 3 ), c1, s1, v( 2 ) )
323 CALL slartg( v( 1 ), v( 2 ), c2, s2, temp )
325 CALL srot( ns, a( ilo+1, ilo ), lda, a( ilo+2, ilo ), lda, c1,
327 CALL srot( ns, a( ilo, ilo ), lda, a( ilo+1, ilo ), lda, c2,
329 CALL srot( ns, b( ilo+1, ilo ), ldb, b( ilo+2, ilo ), ldb, c1,
331 CALL srot( ns, b( ilo, ilo ), ldb, b( ilo+1, ilo ), ldb, c2,
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 )
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 )
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,
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,
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 )
371 sheight = ilo-1-istartm+1
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,
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,
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 )
394 DO WHILE ( k < ihi-ns )
395 np = min( ihi-ns-k, npos )
403 CALL slaset(
'FULL', ns+np, ns+np, zero, one, qc, ldqc )
404 CALL slaset(
'FULL', ns+np, ns+np, zero, one, zc, ldzc )
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 )
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,
429 CALL slacpy(
'ALL', sheight, swidth, work, sheight, a( k+1,
431 CALL sgemm(
'T',
'N', sheight, swidth, sheight, one, qc,
432 $ ldqc, b( k+1, k+ns+np ), ldb, zero, work,
434 CALL slacpy(
'ALL', sheight, swidth, work, sheight, b( k+1,
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 )
445 sheight = k-istartm+1
447 IF ( sheight > 0 )
THEN
448 CALL sgemm(
'N',
'N', sheight, swidth, swidth, one,
449 $ a( istartm, k ), lda, zc, ldzc, zero, work,
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,
456 CALL slacpy(
'ALL', sheight, swidth, work, sheight,
457 $ b( istartm, k ), ldb )
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 )
472 CALL slaset(
'FULL', ns, ns, zero, one, qc, ldqc )
473 CALL slaset(
'FULL', ns+1, ns+1, zero, one, zc, ldzc )
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,
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 )
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 )
514 sheight = ihi-ns-istartm+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,
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,
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 )
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.
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine slartg(f, g, c, s, r)
SLARTG generates a plane rotation with real cosine and real sine.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine slaqz2(ILQ, ILZ, K, ISTARTM, ISTOPM, IHI, A, LDA, B, LDB, NQ, QSTART, Q, LDQ, NZ, ZSTART, Z, LDZ)
SLAQZ2
subroutine slaqz4(ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS, NBLOCK_DESIRED, SR, SI, SS, A, LDA, B, LDB, Q, LDQ, Z, LDZ, QC, LDQC, ZC, LDZC, WORK, LWORK, INFO)
SLAQZ4
subroutine slaqz1(A, LDA, B, LDB, SR1, SR2, SI, BETA1, BETA2, V)
SLAQZ1
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM