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 )
598 CALL slaset(
'Full', m, m, zero, zero, work, m )
601 idum = lwork - m*m - 2
603 CALL slagv2( a( j1, j1 ), lda, b( j1, j1 ), ldb, ar, ai, be,
604 $ work( 1 ), work( 2 ), t( 1, 1 ), t( 2, 1 ) )
605 work( m+1 ) = -work( 2 )
606 work( m+2 ) = work( 1 )
607 t( n2, n2 ) = t( 1, 1 )
608 t( 1, 2 ) = -t( 2, 1 )
614 CALL slagv2( a( j1+n2, j1+n2 ), lda, b( j1+n2, j1+n2 ), ldb,
615 $ taur, taul, work( m*m+1 ), work( n2*m+n2+1 ),
616 $ work( n2*m+n2+2 ), t( n2+1, n2+1 ),
618 work( m*m ) = work( n2*m+n2+1 )
619 work( m*m-1 ) = -work( n2*m+n2+2 )
620 t( m, m ) = t( n2+1, n2+1 )
621 t( m-1, m ) = -t( m, m-1 )
623 CALL sgemm(
'T',
'N', n2, n1, n2, one, work, m, a( j1, j1+n2 ),
624 $ lda, zero, work( m*m+1 ), n2 )
625 CALL slacpy(
'Full', n2, n1, work( m*m+1 ), n2, a( j1, j1+n2 ),
627 CALL sgemm(
'T',
'N', n2, n1, n2, one, work, m, b( j1, j1+n2 ),
628 $ ldb, zero, work( m*m+1 ), n2 )
629 CALL slacpy(
'Full', n2, n1, work( m*m+1 ), n2, b( j1, j1+n2 ),
631 CALL sgemm(
'N',
'N', m, m, m, one, li, ldst, work, m, zero,
633 CALL slacpy(
'Full', m, m, work( m*m+1 ), m, li, ldst )
634 CALL sgemm(
'N',
'N', n2, n1, n1, one, a( j1, j1+n2 ), lda,
635 $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
636 CALL slacpy(
'Full', n2, n1, work, n2, a( j1, j1+n2 ), lda )
637 CALL sgemm(
'N',
'N', n2, n1, n1, one, b( j1, j1+n2 ), ldb,
638 $ t( n2+1, n2+1 ), ldst, zero, work, n2 )
639 CALL slacpy(
'Full', n2, n1, work, n2, b( j1, j1+n2 ), ldb )
640 CALL sgemm(
'T',
'N', m, m, m, one, ir, ldst, t, ldst, zero,
642 CALL slacpy(
'Full', m, m, work, m, ir, ldst )
647 CALL sgemm(
'N',
'N', n, m, m, one, q( 1, j1 ), ldq, li,
648 $ ldst, zero, work, n )
649 CALL slacpy(
'Full', n, m, work, n, q( 1, j1 ), ldq )
654 CALL sgemm(
'N',
'N', n, m, m, one, z( 1, j1 ), ldz, ir,
655 $ ldst, zero, work, n )
656 CALL slacpy(
'Full', n, m, work, n, z( 1, j1 ), ldz )
665 CALL sgemm(
'T',
'N', m, n-i+1, m, one, li, ldst,
666 $ a( j1, i ), lda, zero, work, m )
667 CALL slacpy(
'Full', m, n-i+1, work, m, a( j1, i ), lda )
668 CALL sgemm(
'T',
'N', m, n-i+1, m, one, li, ldst,
669 $ b( j1, i ), ldb, zero, work, m )
670 CALL slacpy(
'Full', m, n-i+1, work, m, b( j1, i ), ldb )
674 CALL sgemm(
'N',
'N', i, m, m, one, a( 1, j1 ), lda, ir,
675 $ ldst, zero, work, i )
676 CALL slacpy(
'Full', i, m, work, i, a( 1, j1 ), lda )
677 CALL sgemm(
'N',
'N', i, m, m, one, b( 1, j1 ), ldb, ir,
678 $ ldst, zero, work, i )
679 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 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 sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine srot(N, SX, INCX, SY, INCY, C, S)
SROT
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 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 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 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 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 sgerq2(M, N, A, LDA, TAU, WORK, INFO)
SGERQ2 computes the RQ factorization of a general rectangular matrix using an unblocked algorithm...
subroutine sscal(N, SA, SX, INCX)
SSCAL
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...
real function slamch(CMACH)
SLAMCH
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...