298 SUBROUTINE dtgsyl( 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,
311 DOUBLE PRECISION dif, scale
315 DOUBLE PRECISION a( lda, * ),
b( ldb, * ), c( ldc, * ),
316 $ d( ldd, * ), e(
lde, * ), f( ldf, * ),
325 DOUBLE PRECISION zero, one
326 parameter( zero = 0.0d+0, one = 1.0d+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 DOUBLE PRECISION dscale, dsum, scale2, scaloc
343 INTRINSIC dble, max, 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(
'DTGSYL', -info )
400 ELSE IF( lquery )
THEN
406 IF( m.EQ.0 .OR. n.EQ.0 )
THEN
418 mb =
ilaenv( 2,
'DTGSYL', trans, m, n, -1, -1 )
419 nb =
ilaenv( 5,
'DTGSYL', trans, m, n, -1, -1 )
426 CALL
dlaset(
'F', m, n, zero, zero, c, ldc )
427 CALL
dlaset(
'F', m, n, zero, zero, f, ldf )
428 ELSE IF( ijob.GE.1 )
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
dtgsy2( 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( dble( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
450 dif = sqrt( dble( pq ) ) / ( dscale*sqrt( dsum ) )
454 IF( isolve.EQ.2 .AND. iround.EQ.1 )
THEN
459 CALL
dlacpy(
'F', m, n, c, ldc, work, m )
460 CALL
dlacpy(
'F', m, n, f, ldf, work( m*n+1 ), m )
461 CALL
dlaset(
'F', m, n, zero, zero, c, ldc )
462 CALL
dlaset(
'F', m, n, zero, zero, f, ldf )
463 ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 )
THEN
464 CALL
dlacpy(
'F', m, n, work, m, c, ldc )
465 CALL
dlacpy(
'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
dtgsy2( 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
dscal( m, scaloc, c( 1, k ), 1 )
549 CALL
dscal( m, scaloc, f( 1, k ), 1 )
552 CALL
dscal( is-1, scaloc, c( 1, k ), 1 )
553 CALL
dscal( is-1, scaloc, f( 1, k ), 1 )
556 CALL
dscal( m-ie, scaloc, c( ie+1, k ), 1 )
557 CALL
dscal( m-ie, scaloc, f( ie+1, k ), 1 )
560 CALL
dscal( m, scaloc, c( 1, k ), 1 )
561 CALL
dscal( m, scaloc, f( 1, k ), 1 )
570 CALL
dgemm(
'N',
'N', is-1, nb, mb, -one,
571 $ a( 1, is ), lda, c( is,
js ), ldc, one,
573 CALL
dgemm(
'N',
'N', is-1, nb, mb, -one,
574 $ d( 1, is ), ldd, c( is,
js ), ldc, one,
578 CALL
dgemm(
'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
dgemm(
'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( dble( 2*m*n ) ) / ( dscale*sqrt( dsum ) )
591 dif = sqrt( dble( pq ) ) / ( dscale*sqrt( dsum ) )
594 IF( isolve.EQ.2 .AND. iround.EQ.1 )
THEN
599 CALL
dlacpy(
'F', m, n, c, ldc, work, m )
600 CALL
dlacpy(
'F', m, n, f, ldf, work( m*n+1 ), m )
601 CALL
dlaset(
'F', m, n, zero, zero, c, ldc )
602 CALL
dlaset(
'F', m, n, zero, zero, f, ldf )
603 ELSE IF( isolve.EQ.2 .AND. iround.EQ.2 )
THEN
604 CALL
dlacpy(
'F', m, n, work, m, c, ldc )
605 CALL
dlacpy(
'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
dtgsy2( 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
dscal( m, scaloc, c( 1, k ), 1 )
636 CALL
dscal( m, scaloc, f( 1, k ), 1 )
639 CALL
dscal( is-1, scaloc, c( 1, k ), 1 )
640 CALL
dscal( is-1, scaloc, f( 1, k ), 1 )
643 CALL
dscal( m-ie, scaloc, c( ie+1, k ), 1 )
644 CALL
dscal( m-ie, scaloc, f( ie+1, k ), 1 )
647 CALL
dscal( m, scaloc, c( 1, k ), 1 )
648 CALL
dscal( m, scaloc, f( 1, k ), 1 )
656 CALL
dgemm(
'N',
'T', mb,
js-1, nb, one, c( is,
js ),
657 $ ldc,
b( 1,
js ), ldb, one, f( is, 1 ),
659 CALL
dgemm(
'N',
'T', mb,
js-1, nb, one, f( is,
js ),
660 $ ldf, e( 1,
js ),
lde, one, f( is, 1 ),
664 CALL
dgemm(
'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
dgemm(
'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 lsame(CA, CB)
LSAME
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dtgsy2(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL, IWORK, PQ, INFO)
DTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
LOGICAL function lde(RI, RJ, LR)
subroutine xerbla(SRNAME, INFO)
XERBLA
for consistency with libblas dev *Do not ship a copy of jquery js
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine dscal(N, DA, DX, INCX)
DSCAL
INTEGER function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine dtgsyl(TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK, IWORK, INFO)
DTGSYL