298 SUBROUTINE stgsyl( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
299 $ ldd, e,
lde, f, ldf, scale, dif, work, lwork,
309 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF,
315 REAL A( lda, * ), B( ldb, * ), C( ldc, * ),
316 $ d( ldd, * ), e( lde, * ), f( ldf, * ),
326 parameter ( zero = 0.0e+0, one = 1.0e+0 )
329 LOGICAL LQUERY, NOTRAN
330 INTEGER I, IE, IFUNC, IROUND, IS, ISOLVE, J, JE, JS, K,
331 $ linfo, lwmin, mb, nb, p, ppqq, pq, q
332 REAL DSCALE, DSUM, SCALE2, SCALOC
337 EXTERNAL lsame, ilaenv
343 INTRINSIC max,
REAL, SQRT
350 notran = lsame( trans,
'N' )
351 lquery = ( lwork.EQ.-1 )
353 IF( .NOT.notran .AND. .NOT.lsame( trans,
'T' ) )
THEN
355 ELSE IF( notran )
THEN
356 IF( ( ijob.LT.0 ) .OR. ( ijob.GT.4 ) )
THEN
363 ELSE IF( n.LE.0 )
THEN
365 ELSE IF( lda.LT.max( 1, m ) )
THEN
367 ELSE IF( ldb.LT.max( 1, n ) )
THEN
369 ELSE IF( ldc.LT.max( 1, m ) )
THEN
371 ELSE IF( ldd.LT.max( 1, m ) )
THEN
373 ELSE IF( lde.LT.max( 1, n ) )
THEN
375 ELSE IF( ldf.LT.max( 1, m ) )
THEN
382 IF( ijob.EQ.1 .OR. ijob.EQ.2 )
THEN
383 lwmin = max( 1, 2*m*n )
392 IF( lwork.LT.lwmin .AND. .NOT.lquery )
THEN
398 CALL xerbla(
'STGSYL', -info )
400 ELSE IF( lquery )
THEN
406 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
418 mb = ilaenv( 2,
'STGSYL', trans, m, n, -1, -1 )
419 nb = ilaenv( 5,
'STGSYL', trans, m, n, -1, -1 )
426 CALL slaset(
'F', m, n, zero, zero, c, ldc )
427 CALL slaset(
'F', m, n, zero, zero, f, ldf )
428 ELSE IF( ijob.GE.1 .AND. notran )
THEN
433 IF( ( mb.LE.1 .AND. nb.LE.1 ) .OR. ( mb.GE.m .AND. nb.GE.n ) )
436 DO 30 iround = 1, isolve
443 CALL stgsy2( trans, ifunc, m, n, a, lda, b, ldb, c, ldc, d,
444 $ ldd, e, lde, f, ldf, scale, dsum, dscale,
446 IF( dscale.NE.zero )
THEN
447 IF( ijob.EQ.1 .OR. ijob.EQ.3 )
THEN
448 dif = sqrt(
REAL( 2*M*N ) ) / ( DSCALE*SQRT( dsum ) )
450 dif = sqrt(
REAL( PQ ) ) / ( DSCALE*SQRT( dsum ) )
454 IF( isolve.EQ.2 .AND. iround.EQ.1 )
THEN
459 CALL slacpy(
'F', m, n, c, ldc, work, m )
460 CALL slacpy(
'F', m, n, f, ldf, work( m*n+1 ), m )
461 CALL slaset(
'F', m, n, zero, zero, c, ldc )
462 CALL slaset(
'F', m, n, zero, zero, f, ldf )
463 ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 )
THEN
464 CALL slacpy(
'F', m, n, work, m, c, ldc )
465 CALL slacpy(
'F', m, n, work( m*n+1 ), m, f, ldf )
485 IF( a( i, i-1 ).NE.zero )
491 IF( iwork( p ).EQ.iwork( p+1 ) )
506 IF( b( j, j-1 ).NE.zero )
512 IF( iwork( q ).EQ.iwork( q+1 ) )
517 DO 150 iround = 1, isolve
530 je = iwork( j+1 ) - 1
534 ie = iwork( i+1 ) - 1
537 CALL stgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
538 $ b( js, js ), ldb, c( is, js ), ldc,
539 $ d( is, is ), ldd, e( js, js ), lde,
540 $ f( is, js ), ldf, scaloc, dsum, dscale,
541 $ iwork( q+2 ), ppqq, linfo )
546 IF( scaloc.NE.one )
THEN
548 CALL sscal( m, scaloc, c( 1, k ), 1 )
549 CALL sscal( m, scaloc, f( 1, k ), 1 )
552 CALL sscal( is-1, scaloc, c( 1, k ), 1 )
553 CALL sscal( is-1, scaloc, f( 1, k ), 1 )
556 CALL sscal( m-ie, scaloc, c( ie+1, k ), 1 )
557 CALL sscal( m-ie, scaloc, f( ie+1, k ), 1 )
560 CALL sscal( m, scaloc, c( 1, k ), 1 )
561 CALL sscal( m, scaloc, f( 1, k ), 1 )
570 CALL sgemm(
'N',
'N', is-1, nb, mb, -one,
571 $ a( 1, is ), lda, c( is, js ), ldc, one,
573 CALL sgemm(
'N',
'N', is-1, nb, mb, -one,
574 $ d( 1, is ), ldd, c( is, js ), ldc, one,
578 CALL sgemm(
'N',
'N', mb, n-je, nb, one,
579 $ f( is, js ), ldf, b( js, je+1 ), ldb,
580 $ one, c( is, je+1 ), ldc )
581 CALL sgemm(
'N',
'N', mb, n-je, nb, one,
582 $ f( is, js ), ldf, e( js, je+1 ), lde,
583 $ one, f( is, je+1 ), ldf )
587 IF( dscale.NE.zero )
THEN
588 IF( ijob.EQ.1 .OR. ijob.EQ.3 )
THEN
589 dif = sqrt(
REAL( 2*M*N ) ) / ( DSCALE*SQRT( dsum ) )
591 dif = sqrt(
REAL( PQ ) ) / ( DSCALE*SQRT( dsum ) )
594 IF( isolve.EQ.2 .AND. iround.EQ.1 )
THEN
599 CALL slacpy(
'F', m, n, c, ldc, work, m )
600 CALL slacpy(
'F', m, n, f, ldf, work( m*n+1 ), m )
601 CALL slaset(
'F', m, n, zero, zero, c, ldc )
602 CALL slaset(
'F', m, n, zero, zero, f, ldf )
603 ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 )
THEN
604 CALL slacpy(
'F', m, n, work, m, c, ldc )
605 CALL slacpy(
'F', m, n, work( m*n+1 ), m, f, ldf )
620 ie = iwork( i+1 ) - 1
622 DO 200 j = q, p + 2, -1
624 je = iwork( j+1 ) - 1
626 CALL stgsy2( trans, ifunc, mb, nb, a( is, is ), lda,
627 $ b( js, js ), ldb, c( is, js ), ldc,
628 $ d( is, is ), ldd, e( js, js ), lde,
629 $ f( is, js ), ldf, scaloc, dsum, dscale,
630 $ iwork( q+2 ), ppqq, linfo )
633 IF( scaloc.NE.one )
THEN
635 CALL sscal( m, scaloc, c( 1, k ), 1 )
636 CALL sscal( m, scaloc, f( 1, k ), 1 )
639 CALL sscal( is-1, scaloc, c( 1, k ), 1 )
640 CALL sscal( is-1, scaloc, f( 1, k ), 1 )
643 CALL sscal( m-ie, scaloc, c( ie+1, k ), 1 )
644 CALL sscal( m-ie, scaloc, f( ie+1, k ), 1 )
647 CALL sscal( m, scaloc, c( 1, k ), 1 )
648 CALL sscal( m, scaloc, f( 1, k ), 1 )
656 CALL sgemm(
'N',
'T', mb, js-1, nb, one, c( is, js ),
657 $ ldc, b( 1, js ), ldb, one, f( is, 1 ),
659 CALL sgemm(
'N',
'T', mb, js-1, nb, one, f( is, js ),
660 $ ldf, e( 1, js ), lde, one, f( is, 1 ),
664 CALL sgemm(
'T',
'N', m-ie, nb, mb, -one,
665 $ a( is, ie+1 ), lda, c( is, js ), ldc, one,
666 $ c( ie+1, js ), ldc )
667 CALL sgemm(
'T',
'N', m-ie, nb, mb, -one,
668 $ d( is, ie+1 ), ldd, f( is, js ), ldf, one,
669 $ c( ie+1, js ), ldc )
logical function lde(RI, RJ, LR)
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine stgsyl(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO)
STGSYL
subroutine xerbla(SRNAME, INFO)
XERBLA
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 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 sscal(N, SA, SX, INCX)
SSCAL
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.