231 INTEGER info, j1, lda, ldb, ldq, ldz, lwork, n, n1, n2
234 REAL a( lda, * ), b( ldb, * ), q( ldq, * ),
235 $ work( * ), z( ldz, * )
244 parameter( zero = 0.0e+0, one = 1.0e+0 )
246 parameter( twenty = 2.0e+01 )
248 parameter( ldst = 4 )
250 parameter( wands = .true. )
254 INTEGER i, idum, linfo, m
255 REAL bqra21, brqa21, ddum, dnorm, dscale, dsum, eps,
256 $ f, g, sa, sb, scale, smlnum, ss, thresh, ws
259 INTEGER iwork( ldst )
260 REAL ai( 2 ), ar( 2 ), be( 2 ), ir( ldst, ldst ),
261 $ ircop( ldst, ldst ), li( ldst, ldst ),
262 $ licop( ldst, ldst ), s( ldst, ldst ),
263 $ scpy( ldst, ldst ), t( ldst, ldst ),
264 $ taul( ldst ), taur( ldst ), tcpy( ldst, ldst )
276 INTRINSIC abs, max, sqrt
284 IF( n.LE.1 .OR. n1.LE.0 .OR. n2.LE.0 )
286 IF( n1.GT.n .OR. ( j1+n1 ).GT.n )
289 IF( lwork.LT.max( n*m, m*m*2 ) )
THEN
291 work( 1 ) = max( n*m, m*m*2 )
300 CALL slaset(
'Full', ldst, ldst, zero, zero, li, ldst )
301 CALL slaset(
'Full', ldst, ldst, zero, zero, ir, ldst )
302 CALL slacpy(
'Full', m, m, a( j1, j1 ), lda, s, ldst )
303 CALL slacpy(
'Full', m, m, b( j1, j1 ), ldb, t, ldst )
308 smlnum =
slamch(
'S' ) / eps
311 CALL slacpy(
'Full', m, m, s, ldst, work, m )
312 CALL slassq( m*m, work, 1, dscale, dsum )
313 CALL slacpy(
'Full', m, m, t, ldst, work, m )
314 CALL slassq( m*m, work, 1, dscale, dsum )
315 dnorm = dscale*sqrt( dsum )
325 thresh = max( twenty*eps*dnorm, smlnum )
334 f = s( 2, 2 )*t( 1, 1 ) - t( 2, 2 )*s( 1, 1 )
335 g = s( 2, 2 )*t( 1, 2 ) - t( 2, 2 )*s( 1, 2 )
336 sb = abs( t( 2, 2 ) )
337 sa = abs( s( 2, 2 ) )
338 CALL slartg( f, g, ir( 1, 2 ), ir( 1, 1 ), ddum )
339 ir( 2, 1 ) = -ir( 1, 2 )
340 ir( 2, 2 ) = ir( 1, 1 )
341 CALL srot( 2, s( 1, 1 ), 1, s( 1, 2 ), 1, ir( 1, 1 ),
343 CALL srot( 2, t( 1, 1 ), 1, t( 1, 2 ), 1, ir( 1, 1 ),
346 CALL slartg( s( 1, 1 ), s( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
349 CALL slartg( t( 1, 1 ), t( 2, 1 ), li( 1, 1 ), li( 2, 1 ),
352 CALL srot( 2, s( 1, 1 ), ldst, s( 2, 1 ), ldst, li( 1, 1 ),
354 CALL srot( 2, t( 1, 1 ), ldst, t( 2, 1 ), ldst, li( 1, 1 ),
356 li( 2, 2 ) = li( 1, 1 )
357 li( 1, 2 ) = -li( 2, 1 )
362 ws = abs( s( 2, 1 ) ) + abs( t( 2, 1 ) )
372 CALL slacpy(
'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
374 CALL sgemm(
'N',
'N', m, m, m, one, li, ldst, s, ldst, zero,
376 CALL sgemm(
'N',
'T', m, m, m, -one, work, m, ir, ldst, one,
380 CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
382 CALL slacpy(
'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
384 CALL sgemm(
'N',
'N', m, m, m, one, li, ldst, t, ldst, zero,
386 CALL sgemm(
'N',
'T', m, m, m, -one, work, m, ir, ldst, one,
388 CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
389 ss = dscale*sqrt( dsum )
390 strong = ss.LE.thresh
398 CALL srot( j1+1, a( 1, j1 ), 1, a( 1, j1+1 ), 1, ir( 1, 1 ),
400 CALL srot( j1+1, b( 1, j1 ), 1, b( 1, j1+1 ), 1, ir( 1, 1 ),
402 CALL srot( n-j1+1, a( j1, j1 ), lda, a( j1+1, j1 ), lda,
403 $ li( 1, 1 ), li( 2, 1 ) )
404 CALL srot( n-j1+1, b( j1, j1 ), ldb, b( j1+1, j1 ), ldb,
405 $ li( 1, 1 ), li( 2, 1 ) )
415 $
CALL srot( n, z( 1, j1 ), 1, z( 1, j1+1 ), 1, ir( 1, 1 ),
418 $
CALL srot( n, q( 1, j1 ), 1, q( 1, j1+1 ), 1, li( 1, 1 ),
435 CALL slacpy(
'Full', n1, n2, t( 1, n1+1 ), ldst, li, ldst )
436 CALL slacpy(
'Full', n1, n2, s( 1, n1+1 ), ldst,
437 $ ir( n2+1, n1+1 ), ldst )
438 CALL stgsy2(
'N', 0, n1, n2, s, ldst, s( n1+1, n1+1 ), ldst,
439 $ ir( n2+1, n1+1 ), ldst, t, ldst, t( n1+1, n1+1 ),
440 $ ldst, li, ldst, scale, dsum, dscale, iwork, idum,
452 CALL sscal( n1, -one, li( 1, i ), 1 )
453 li( n1+i, i ) = scale
455 CALL sgeqr2( m, n2, li, ldst, taul, work, linfo )
458 CALL sorg2r( m, m, n2, li, ldst, taul, work, linfo )
469 ir( n2+i, i ) = scale
471 CALL sgerq2( n1, m, ir( n2+1, 1 ), ldst, taur, work, linfo )
474 CALL sorgr2( m, m, n1, ir, ldst, taur, work, linfo )
480 CALL sgemm(
'T',
'N', m, m, m, one, li, ldst, s, ldst, zero,
482 CALL sgemm(
'N',
'T', m, m, m, one, work, m, ir, ldst, zero, s,
484 CALL sgemm(
'T',
'N', m, m, m, one, li, ldst, t, ldst, zero,
486 CALL sgemm(
'N',
'T', m, m, m, one, work, m, ir, ldst, zero, t,
488 CALL slacpy(
'F', m, m, s, ldst, scpy, ldst )
489 CALL slacpy(
'F', m, m, t, ldst, tcpy, ldst )
490 CALL slacpy(
'F', m, m, ir, ldst, ircop, ldst )
491 CALL slacpy(
'F', m, m, li, ldst, licop, ldst )
496 CALL sgerq2( m, m, t, ldst, taur, work, linfo )
499 CALL sormr2(
'R',
'T', m, m, m, t, ldst, taur, s, ldst, work,
503 CALL sormr2(
'L',
'N', m, m, m, t, ldst, taur, ir, ldst, work,
513 CALL slassq( n1, s( n2+1, i ), 1, dscale, dsum )
515 brqa21 = dscale*sqrt( dsum )
520 CALL sgeqr2( m, m, tcpy, ldst, taul, work, linfo )
523 CALL sorm2r(
'L',
'T', m, m, m, tcpy, ldst, taul, scpy, ldst,
525 CALL sorm2r(
'R',
'N', m, m, m, tcpy, ldst, taul, licop, ldst,
535 CALL slassq( n1, scpy( n2+1, i ), 1, dscale, dsum )
537 bqra21 = dscale*sqrt( dsum )
543 IF( bqra21.LE.brqa21 .AND. bqra21.LE.thresh )
THEN
544 CALL slacpy(
'F', m, m, scpy, ldst, s, ldst )
545 CALL slacpy(
'F', m, m, tcpy, ldst, t, ldst )
546 CALL slacpy(
'F', m, m, ircop, ldst, ir, ldst )
547 CALL slacpy(
'F', m, m, licop, ldst, li, ldst )
548 ELSE IF( brqa21.GE.thresh )
THEN
554 CALL slaset(
'Lower', m-1, m-1, zero, zero, t(2,1), ldst )
561 CALL slacpy(
'Full', m, m, a( j1, j1 ), lda, work( m*m+1 ),
563 CALL sgemm(
'N',
'N', m, m, m, one, li, ldst, s, ldst, zero,
565 CALL sgemm(
'N',
'N', m, m, m, -one, work, m, ir, ldst, one,
569 CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
571 CALL slacpy(
'Full', m, m, b( j1, j1 ), ldb, work( m*m+1 ),
573 CALL sgemm(
'N',
'N', m, m, m, one, li, ldst, t, ldst, zero,
575 CALL sgemm(
'N',
'N', m, m, m, -one, work, m, ir, ldst, one,
577 CALL slassq( m*m, work( m*m+1 ), 1, dscale, dsum )
578 ss = dscale*sqrt( dsum )
579 strong = ( ss.LE.thresh )
588 CALL slaset(
'Full', n1, n2, zero, zero, s(n2+1,1), ldst )
592 CALL slacpy(
'F', m, m, s, ldst, a( j1, j1 ), lda )
593 CALL slacpy(
'F', m, m, t, ldst, b( j1, j1 ), ldb )
594 CALL slaset(
'Full', ldst, ldst, zero, zero, t, ldst )
603 idum = lwork - m*m - 2
605 CALL slagv2( a( j1, j1 ), lda, b( j1, j1 ), ldb, ar, ai, be,
606 $ work( 1 ), work( 2 ), t( 1, 1 ), t( 2, 1 ) )
607 work( m+1 ) = -work( 2 )
608 work( m+2 ) = work( 1 )
609 t( n2, n2 ) = t( 1, 1 )
610 t( 1, 2 ) = -t( 2, 1 )
616 CALL slagv2( a( j1+n2, j1+n2 ), lda, b( j1+n2, j1+n2 ), ldb,
617 $ taur, taul, work( m*m+1 ), work( n2*m+n2+1 ),
618 $ work( n2*m+n2+2 ), t( n2+1, n2+1 ),
620 work( m*m ) = work( n2*m+n2+1 )
621 work( m*m-1 ) = -work( n2*m+n2+2 )
622 t( m, m ) = t( n2+1, n2+1 )
623 t( m-1, m ) = -t( m, m-1 )
625 CALL sgemm(
'T',
'N', n2, n1, n2, one, work, m, a( j1, j1+n2 ),
626 $ lda, zero, work( m*m+1 ), n2 )
627 CALL slacpy(
'Full', n2, n1, work( m*m+1 ), n2, a( j1, j1+n2 ),
629 CALL sgemm(
'T',
'N', n2, n1, n2, one, work, m, b( j1, j1+n2 ),
630 $ ldb, zero, work( m*m+1 ), n2 )
631 CALL slacpy(
'Full', n2, n1, work( m*m+1 ), n2, b( j1, j1+n2 ),
633 CALL sgemm(
'N',
'N', m, m, m, one, li, ldst, work, m, zero,
635 CALL slacpy(
'Full', m, m, work( m*m+1 ), m, li, ldst )
636 CALL sgemm(
'N',
'N', n2, n1, n1, one, a( j1, j1+n2 ), lda,
637 $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
638 CALL slacpy(
'Full', n2, n1, work, n2, a( j1, j1+n2 ), lda )
639 CALL sgemm(
'N',
'N', n2, n1, n1, one, b( j1, j1+n2 ), ldb,
640 $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
641 CALL slacpy(
'Full', n2, n1, work, n2, b( j1, j1+n2 ), ldb )
642 CALL sgemm(
'T',
'N', m, m, m, one, ir, ldst, t, ldst, zero,
644 CALL slacpy(
'Full', m, m, work, m, ir, ldst )
649 CALL sgemm(
'N',
'N', n, m, m, one, q( 1, j1 ), ldq, li,
650 $ ldst, zero, work, n )
651 CALL slacpy(
'Full', n, m, work, n, q( 1, j1 ), ldq )
656 CALL sgemm(
'N',
'N', n, m, m, one, z( 1, j1 ), ldz, ir,
657 $ ldst, zero, work, n )
658 CALL slacpy(
'Full', n, m, work, n, z( 1, j1 ), ldz )
667 CALL sgemm(
'T',
'N', m, n-i+1, m, one, li, ldst,
668 $ a( j1, i ), lda, zero, work, m )
669 CALL slacpy(
'Full', m, n-i+1, work, m, a( j1, i ), lda )
670 CALL sgemm(
'T',
'N', m, n-i+1, m, one, li, ldst,
671 $ b( j1, i ), ldb, zero, work, m )
672 CALL slacpy(
'Full', m, n-i+1, work, m, b( j1, i ), ldb )
676 CALL sgemm(
'N',
'N', i, m, m, one, a( 1, j1 ), lda, ir,
677 $ ldst, zero, work, i )
678 CALL slacpy(
'Full', i, m, work, i, a( 1, j1 ), lda )
679 CALL sgemm(
'N',
'N', i, m, m, one, b( 1, j1 ), ldb, ir,
680 $ ldst, zero, work, i )
681 CALL slacpy(
'Full', i, m, work, i, b( 1, j1 ), ldb )
subroutine slassq(N, X, INCX, SCALE, SUMSQ)
SLASSQ updates a sum of squares represented in scaled form.
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 sorgr2(M, N, K, A, LDA, TAU, WORK, INFO)
SORGR2 generates all or part of the orthogonal matrix Q from an RQ factorization determined by sgerqf...
subroutine sormr2(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
SORMR2 multiplies a general matrix by the orthogonal matrix from a RQ factorization determined by sge...
subroutine stgsy2(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, IWORK, PQ, INFO)
STGSY2 solves the generalized Sylvester equation (unblocked algorithm).
subroutine slagv2(A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, CSL, SNL, CSR, SNR)
SLAGV2 computes the Generalized Schur factorization of a real 2-by-2 matrix pencil (A...
subroutine sorg2r(M, N, K, A, LDA, TAU, WORK, INFO)
SORG2R generates all or part of the orthogonal matrix Q from a QR factorization determined by sgeqrf ...
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine sorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
SORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
real function slamch(CMACH)
SLAMCH
subroutine sgerq2(M, N, A, LDA, TAU, WORK, INFO)
SGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm...
subroutine slartg(F, G, CS, SN, R)
SLARTG generates a plane rotation with real cosine and real sine.
subroutine sgeqr2(M, N, A, LDA, TAU, WORK, INFO)
SGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm...
subroutine sscal(N, SA, SX, INCX)
SSCAL
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT