82
   83
   84
   85
   86
   87
   88      INTEGER            KNT, LMAX, NINFO
   89      REAL               RMAX
   90
   91
   92
   93
   94
   95      REAL               ZERO, ONE
   96      parameter( zero = 0.0e0, one = 1.0e0 )
   97      REAL               TWO, FOUR, EIGHT
   98      parameter( two = 2.0e0, four = 4.0e0, eight = 8.0e0 )
   99
  100
  101      LOGICAL            LTRANL, LTRANR
  102      INTEGER            IB, IB1, IB2, IB3, INFO, ISGN, ITL, ITLSCL,
  103     $                   ITR, ITRANL, ITRANR, ITRSCL, N1, N2
  104      REAL               BIGNUM, DEN, EPS, RES, SCALE, SGN, SMLNUM, TMP,
  105     $                   TNRM, XNORM, XNRM
  106
  107
  108      INTEGER            ITVAL( 2, 2, 8 )
  109      REAL               B( 2, 2 ), TL( 2, 2 ), TR( 2, 2 ), VAL( 3 ),
  110     $                   X( 2, 2 )
  111
  112
  113      REAL               SLAMCH
  115
  116
  118
  119
  120      INTRINSIC          abs, max, min, sqrt
  121
  122
  123      DATA               itval / 8, 4, 2, 1, 4, 8, 1, 2, 2, 1, 8, 4, 1,
  124     $                   2, 4, 8, 9, 4, 2, 1, 4, 9, 1, 2, 2, 1, 9, 4, 1,
  125     $                   2, 4, 9 /
  126
  127
  128
  129
  130
  132      smlnum = 
slamch( 
'S' ) / eps
 
  133      bignum = one / smlnum
  134
  135
  136
  137      val( 1 ) = sqrt( smlnum )
  138      val( 2 ) = one
  139      val( 3 ) = sqrt( bignum )
  140
  141      knt = 0
  142      ninfo = 0
  143      lmax = 0
  144      rmax = zero
  145
  146
  147
  148      DO 230 itranl = 0, 1
  149         DO 220 itranr = 0, 1
  150            DO 210 isgn = -1, 1, 2
  151               sgn = isgn
  152               ltranl = itranl.EQ.1
  153               ltranr = itranr.EQ.1
  154
  155               n1 = 1
  156               n2 = 1
  157               DO 30 itl = 1, 3
  158                  DO 20 itr = 1, 3
  159                     DO 10 ib = 1, 3
  160                        tl( 1, 1 ) = val( itl )
  161                        tr( 1, 1 ) = val( itr )
  162                        b( 1, 1 ) = val( ib )
  163                        knt = knt + 1
  164                        CALL slasy2( ltranl, ltranr, isgn, n1, n2, tl,
 
  165     $                               2, tr, 2, b, 2, scale, x, 2, xnorm,
  166     $                               info )
  167                        IF( info.NE.0 )
  168     $                     ninfo = ninfo + 1
  169                        res = abs( ( tl( 1, 1 )+sgn*tr( 1, 1 ) )*
  170     $                        x( 1, 1 )-scale*b( 1, 1 ) )
  171                        IF( info.EQ.0 ) THEN
  172                           den = max( eps*( ( abs( tr( 1,
  173     $                           1 ) )+abs( tl( 1, 1 ) ) )*abs( x( 1,
  174     $                           1 ) ) ), smlnum )
  175                        ELSE
  176                           den = smlnum*max( abs( x( 1, 1 ) ), one )
  177                        END IF
  178                        res = res / den
  179                        IF( scale.GT.one )
  180     $                     res = res + one / eps
  181                        res = res + abs( xnorm-abs( x( 1, 1 ) ) ) /
  182     $                        max( smlnum, xnorm ) / eps
  183                        IF( info.NE.0 .AND. info.NE.1 )
  184     $                     res = res + one / eps
  185                        IF( res.GT.rmax ) THEN
  186                           lmax = knt
  187                           rmax = res
  188                        END IF
  189   10                CONTINUE
  190   20             CONTINUE
  191   30          CONTINUE
  192
  193               n1 = 2
  194               n2 = 1
  195               DO 80 itl = 1, 8
  196                  DO 70 itlscl = 1, 3
  197                     DO 60 itr = 1, 3
  198                        DO 50 ib1 = 1, 3
  199                           DO 40 ib2 = 1, 3
  200                              b( 1, 1 ) = val( ib1 )
  201                              b( 2, 1 ) = -four*val( ib2 )
  202                              tl( 1, 1 ) = itval( 1, 1, itl )*
  203     $                                     val( itlscl )
  204                              tl( 2, 1 ) = itval( 2, 1, itl )*
  205     $                                     val( itlscl )
  206                              tl( 1, 2 ) = itval( 1, 2, itl )*
  207     $                                     val( itlscl )
  208                              tl( 2, 2 ) = itval( 2, 2, itl )*
  209     $                                     val( itlscl )
  210                              tr( 1, 1 ) = val( itr )
  211                              knt = knt + 1
  212                              CALL slasy2( ltranl, ltranr, isgn, n1, n2,
 
  213     $                                     tl, 2, tr, 2, b, 2, scale, x,
  214     $                                     2, xnorm, info )
  215                              IF( info.NE.0 )
  216     $                           ninfo = ninfo + 1
  217                              IF( ltranl ) THEN
  218                                 tmp = tl( 1, 2 )
  219                                 tl( 1, 2 ) = tl( 2, 1 )
  220                                 tl( 2, 1 ) = tmp
  221                              END IF
  222                              res = abs( ( tl( 1, 1 )+sgn*tr( 1, 1 ) )*
  223     $                              x( 1, 1 )+tl( 1, 2 )*x( 2, 1 )-
  224     $                              scale*b( 1, 1 ) )
  225                              res = res + abs( ( tl( 2, 2 )+sgn*tr( 1,
  226     $                              1 ) )*x( 2, 1 )+tl( 2, 1 )*
  227     $                              x( 1, 1 )-scale*b( 2, 1 ) )
  228                              tnrm = abs( tr( 1, 1 ) ) +
  229     $                               abs( tl( 1, 1 ) ) +
  230     $                               abs( tl( 1, 2 ) ) +
  231     $                               abs( tl( 2, 1 ) ) +
  232     $                               abs( tl( 2, 2 ) )
  233                              xnrm = max( abs( x( 1, 1 ) ),
  234     $                               abs( x( 2, 1 ) ) )
  235                              den = max( smlnum, smlnum*xnrm,
  236     $                              ( tnrm*eps )*xnrm )
  237                              res = res / den
  238                              IF( scale.GT.one )
  239     $                           res = res + one / eps
  240                              res = res + abs( xnorm-xnrm ) /
  241     $                              max( smlnum, xnorm ) / eps
  242                              IF( res.GT.rmax ) THEN
  243                                 lmax = knt
  244                                 rmax = res
  245                              END IF
  246   40                      CONTINUE
  247   50                   CONTINUE
  248   60                CONTINUE
  249   70             CONTINUE
  250   80          CONTINUE
  251
  252               n1 = 1
  253               n2 = 2
  254               DO 130 itr = 1, 8
  255                  DO 120 itrscl = 1, 3
  256                     DO 110 itl = 1, 3
  257                        DO 100 ib1 = 1, 3
  258                           DO 90 ib2 = 1, 3
  259                              b( 1, 1 ) = val( ib1 )
  260                              b( 1, 2 ) = -two*val( ib2 )
  261                              tr( 1, 1 ) = itval( 1, 1, itr )*
  262     $                                     val( itrscl )
  263                              tr( 2, 1 ) = itval( 2, 1, itr )*
  264     $                                     val( itrscl )
  265                              tr( 1, 2 ) = itval( 1, 2, itr )*
  266     $                                     val( itrscl )
  267                              tr( 2, 2 ) = itval( 2, 2, itr )*
  268     $                                     val( itrscl )
  269                              tl( 1, 1 ) = val( itl )
  270                              knt = knt + 1
  271                              CALL slasy2( ltranl, ltranr, isgn, n1, n2,
 
  272     $                                     tl, 2, tr, 2, b, 2, scale, x,
  273     $                                     2, xnorm, info )
  274                              IF( info.NE.0 )
  275     $                           ninfo = ninfo + 1
  276                              IF( ltranr ) THEN
  277                                 tmp = tr( 1, 2 )
  278                                 tr( 1, 2 ) = tr( 2, 1 )
  279                                 tr( 2, 1 ) = tmp
  280                              END IF
  281                              tnrm = abs( tl( 1, 1 ) ) +
  282     $                               abs( tr( 1, 1 ) ) +
  283     $                               abs( tr( 1, 2 ) ) +
  284     $                               abs( tr( 2, 2 ) ) +
  285     $                               abs( tr( 2, 1 ) )
  286                              xnrm = abs( x( 1, 1 ) ) + abs( x( 1, 2 ) )
  287                              res = abs( ( ( tl( 1, 1 )+sgn*tr( 1,
  288     $                              1 ) ) )*( x( 1, 1 ) )+
  289     $                              ( sgn*tr( 2, 1 ) )*( x( 1, 2 ) )-
  290     $                              ( scale*b( 1, 1 ) ) )
  291                              res = res + abs( ( ( tl( 1, 1 )+sgn*tr( 2,
  292     $                              2 ) ) )*( x( 1, 2 ) )+
  293     $                              ( sgn*tr( 1, 2 ) )*( x( 1, 1 ) )-
  294     $                              ( scale*b( 1, 2 ) ) )
  295                              den = max( smlnum, smlnum*xnrm,
  296     $                              ( tnrm*eps )*xnrm )
  297                              res = res / den
  298                              IF( scale.GT.one )
  299     $                           res = res + one / eps
  300                              res = res + abs( xnorm-xnrm ) /
  301     $                              max( smlnum, xnorm ) / eps
  302                              IF( res.GT.rmax ) THEN
  303                                 lmax = knt
  304                                 rmax = res
  305                              END IF
  306   90                      CONTINUE
  307  100                   CONTINUE
  308  110                CONTINUE
  309  120             CONTINUE
  310  130          CONTINUE
  311
  312               n1 = 2
  313               n2 = 2
  314               DO 200 itr = 1, 8
  315                  DO 190 itrscl = 1, 3
  316                     DO 180 itl = 1, 8
  317                        DO 170 itlscl = 1, 3
  318                           DO 160 ib1 = 1, 3
  319                              DO 150 ib2 = 1, 3
  320                                 DO 140 ib3 = 1, 3
  321                                    b( 1, 1 ) = val( ib1 )
  322                                    b( 2, 1 ) = -four*val( ib2 )
  323                                    b( 1, 2 ) = -two*val( ib3 )
  324                                    b( 2, 2 ) = eight*
  325     $                                          min( val( ib1 ), val
  326     $                                          ( ib2 ), val( ib3 ) )
  327                                    tr( 1, 1 ) = itval( 1, 1, itr )*
  328     $                                           val( itrscl )
  329                                    tr( 2, 1 ) = itval( 2, 1, itr )*
  330     $                                           val( itrscl )
  331                                    tr( 1, 2 ) = itval( 1, 2, itr )*
  332     $                                           val( itrscl )
  333                                    tr( 2, 2 ) = itval( 2, 2, itr )*
  334     $                                           val( itrscl )
  335                                    tl( 1, 1 ) = itval( 1, 1, itl )*
  336     $                                           val( itlscl )
  337                                    tl( 2, 1 ) = itval( 2, 1, itl )*
  338     $                                           val( itlscl )
  339                                    tl( 1, 2 ) = itval( 1, 2, itl )*
  340     $                                           val( itlscl )
  341                                    tl( 2, 2 ) = itval( 2, 2, itl )*
  342     $                                           val( itlscl )
  343                                    knt = knt + 1
  344                                    CALL slasy2( ltranl, ltranr, isgn,
 
  345     $                                           n1, n2, tl, 2, tr, 2,
  346     $                                           b, 2, scale, x, 2,
  347     $                                           xnorm, info )
  348                                    IF( info.NE.0 )
  349     $                                 ninfo = ninfo + 1
  350                                    IF( ltranr ) THEN
  351                                       tmp = tr( 1, 2 )
  352                                       tr( 1, 2 ) = tr( 2, 1 )
  353                                       tr( 2, 1 ) = tmp
  354                                    END IF
  355                                    IF( ltranl ) THEN
  356                                       tmp = tl( 1, 2 )
  357                                       tl( 1, 2 ) = tl( 2, 1 )
  358                                       tl( 2, 1 ) = tmp
  359                                    END IF
  360                                    tnrm = abs( tr( 1, 1 ) ) +
  361     $                                     abs( tr( 2, 1 ) ) +
  362     $                                     abs( tr( 1, 2 ) ) +
  363     $                                     abs( tr( 2, 2 ) ) +
  364     $                                     abs( tl( 1, 1 ) ) +
  365     $                                     abs( tl( 2, 1 ) ) +
  366     $                                     abs( tl( 1, 2 ) ) +
  367     $                                     abs( tl( 2, 2 ) )
  368                                    xnrm = max( abs( x( 1, 1 ) )+
  369     $                                     abs( x( 1, 2 ) ),
  370     $                                     abs( x( 2, 1 ) )+
  371     $                                     abs( x( 2, 2 ) ) )
  372                                    res = abs( ( ( tl( 1, 1 )+sgn*tr( 1,
  373     $                                    1 ) ) )*( x( 1, 1 ) )+
  374     $                                    ( sgn*tr( 2, 1 ) )*
  375     $                                    ( x( 1, 2 ) )+( tl( 1, 2 ) )*
  376     $                                    ( x( 2, 1 ) )-
  377     $                                    ( scale*b( 1, 1 ) ) )
  378                                    res = res + abs( ( tl( 1, 1 ) )*
  379     $                                    ( x( 1, 2 ) )+
  380     $                                    ( sgn*tr( 1, 2 ) )*
  381     $                                    ( x( 1, 1 ) )+
  382     $                                    ( sgn*tr( 2, 2 ) )*
  383     $                                    ( x( 1, 2 ) )+( tl( 1, 2 ) )*
  384     $                                    ( x( 2, 2 ) )-
  385     $                                    ( scale*b( 1, 2 ) ) )
  386                                    res = res + abs( ( tl( 2, 1 ) )*
  387     $                                    ( x( 1, 1 ) )+
  388     $                                    ( sgn*tr( 1, 1 ) )*
  389     $                                    ( x( 2, 1 ) )+
  390     $                                    ( sgn*tr( 2, 1 ) )*
  391     $                                    ( x( 2, 2 ) )+( tl( 2, 2 ) )*
  392     $                                    ( x( 2, 1 ) )-
  393     $                                    ( scale*b( 2, 1 ) ) )
  394                                    res = res + abs( ( ( tl( 2,
  395     $                                    2 )+sgn*tr( 2, 2 ) ) )*
  396     $                                    ( x( 2, 2 ) )+
  397     $                                    ( sgn*tr( 1, 2 ) )*
  398     $                                    ( x( 2, 1 ) )+( tl( 2, 1 ) )*
  399     $                                    ( x( 1, 2 ) )-
  400     $                                    ( scale*b( 2, 2 ) ) )
  401                                    den = max( smlnum, smlnum*xnrm,
  402     $                                    ( tnrm*eps )*xnrm )
  403                                    res = res / den
  404                                    IF( scale.GT.one )
  405     $                                 res = res + one / eps
  406                                    res = res + abs( xnorm-xnrm ) /
  407     $                                    max( smlnum, xnorm ) / eps
  408                                    IF( res.GT.rmax ) THEN
  409                                       lmax = knt
  410                                       rmax = res
  411                                    END IF
  412  140                            CONTINUE
  413  150                         CONTINUE
  414  160                      CONTINUE
  415  170                   CONTINUE
  416  180                CONTINUE
  417  190             CONTINUE
  418  200          CONTINUE
  419  210       CONTINUE
  420  220    CONTINUE
  421  230 CONTINUE
  422
  423      RETURN
  424
  425
  426
real function slamch(cmach)
SLAMCH
subroutine slasy2(ltranl, ltranr, isgn, n1, n2, tl, ldtl, tr, ldtr, b, ldb, scale, x, ldx, xnorm, info)
SLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2.