sget35.f

Go to the documentation of this file.

*> \brief \b SGET35
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE SGET35( RMAX, LMAX, NINFO, KNT )
*
*       .. Scalar Arguments ..
*       INTEGER            KNT, LMAX, NINFO
*       REAL               RMAX
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SGET35 tests STRSYL, a routine for solving the Sylvester matrix
*> equation
*>
*>    op(A)*X + ISGN*X*op(B) = scale*C,
*>
*> A and B are assumed to be in Schur canonical form, op() represents an
*> optional transpose, and ISGN can be -1 or +1.  Scale is an output
*> less than or equal to 1, chosen to avoid overflow in X.
*>
*> The test code verifies that the following residual is order 1:
*>
*>    norm(op(A)*X + ISGN*X*op(B) - scale*C) /
*>        (EPS*max(norm(A),norm(B))*norm(X))
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[out] RMAX
*> \verbatim
*>          RMAX is REAL
*>          Value of the largest test ratio.
*> \endverbatim
*>
*> \param[out] LMAX
*> \verbatim
*>          LMAX is INTEGER
*>          Example number where largest test ratio achieved.
*> \endverbatim
*>
*> \param[out] NINFO
*> \verbatim
*>          NINFO is INTEGER
*>          Number of examples where INFO is nonzero.
*> \endverbatim
*>
*> \param[out] KNT
*> \verbatim
*>          KNT is INTEGER
*>          Total number of examples tested.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup single_eig
*
*  =====================================================================
      SUBROUTINE sget35( RMAX, LMAX, NINFO, KNT )
*
*  -- LAPACK test routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      INTEGER            KNT, LMAX, NINFO
      REAL               RMAX
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      parameter( zero = 0.0e0, one = 1.0e0 )
      REAL               TWO, FOUR
      parameter( two = 2.0e0, four = 4.0e0 )
*     ..
*     .. Local Scalars ..
      CHARACTER          TRANA, TRANB
      INTEGER            I, IMA, IMB, IMLDA1, IMLDA2, IMLDB1, IMLOFF,
     $                   INFO, ISGN, ITRANA, ITRANB, J, M, N
      REAL               BIGNUM, CNRM, EPS, RES, RES1, RMUL, SCALE,
     $                   SMLNUM, TNRM, XNRM
*     ..
*     .. Local Arrays ..
      INTEGER            IDIM( 8 ), IVAL( 6, 6, 8 )
      REAL               A( 6, 6 ), B( 6, 6 ), C( 6, 6 ), CC( 6, 6 ),
     $                   DUM( 1 ), VM1( 3 ), VM2( 3 )
*     ..
*     .. External Functions ..
      REAL               SLAMCH, SLANGE
      EXTERNAL           slamch, slange
*     ..
*     .. External Subroutines ..
      EXTERNAL           sgemm, strsyl
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, max, real, sin, sqrt
*     ..
*     .. Data statements ..
      DATA               idim / 1, 2, 3, 4, 3, 3, 6, 4 /
      DATA               ival / 1, 35*0, 1, 2, 4*0, -2, 0, 28*0, 1, 5*0,
     $                   5, 1, 2, 3*0, -8, -2, 1, 21*0, 3, 4, 4*0, -5,
     $                   3, 4*0, 1, 2, 1, 4, 2*0, -3, -9, -1, 1, 14*0,
     $                   1, 5*0, 2, 3, 4*0, 5, 6, 7, 21*0, 1, 5*0, 1, 3,
     $                   -4, 3*0, 2, 5, 2, 21*0, 1, 2, 4*0, -2, 0, 4*0,
     $                   5, 6, 3, 4, 2*0, -1, -9, -5, 2, 2*0, 4*8, 5, 6,
     $                   4*9, -7, 5, 1, 5*0, 1, 5, 2, 3*0, 2, -21, 5,
     $                   3*0, 1, 2, 3, 4, 14*0 /
*     ..
*     .. Executable Statements ..
*
*     Get machine parameters
*
      eps = slamch( 'P' )
      smlnum = slamch( 'S' )*four / eps
      bignum = one / smlnum
      CALL slabad( smlnum, bignum )
*
*     Set up test case parameters
*
      vm1( 1 ) = sqrt( smlnum )
      vm1( 2 ) = one
      vm1( 3 ) = sqrt( bignum )
      vm2( 1 ) = one
      vm2( 2 ) = one + two*eps
      vm2( 3 ) = two
*
      knt = 0
      ninfo = 0
      lmax = 0
      rmax = zero
*
*     Begin test loop
*
      DO 150 itrana = 1, 2
         DO 140 itranb = 1, 2
            DO 130 isgn = -1, 1, 2
               DO 120 ima = 1, 8
                  DO 110 imlda1 = 1, 3
                     DO 100 imlda2 = 1, 3
                        DO 90 imloff = 1, 2
                           DO 80 imb = 1, 8
                              DO 70 imldb1 = 1, 3
                                 IF( itrana.EQ.1 )
     $                              trana = 'N'
                                 IF( itrana.EQ.2 )
     $                              trana = 'T'
                                 IF( itranb.EQ.1 )
     $                              tranb = 'N'
                                 IF( itranb.EQ.2 )
     $                              tranb = 'T'
                                 m = idim( ima )
                                 n = idim( imb )
                                 tnrm = zero
                                 DO 20 i = 1, m
                                    DO 10 j = 1, m
                                       a( i, j ) = ival( i, j, ima )
                                       IF( abs( i-j ).LE.1 ) THEN
                                          a( i, j ) = a( i, j )*
     $                                                vm1( imlda1 )
                                          a( i, j ) = a( i, j )*
     $                                                vm2( imlda2 )
                                       ELSE
                                          a( i, j ) = a( i, j )*
     $                                                vm1( imloff )
                                       END IF
                                       tnrm = max( tnrm,
     $                                        abs( a( i, j ) ) )
   10                               CONTINUE
   20                            CONTINUE
                                 DO 40 i = 1, n
                                    DO 30 j = 1, n
                                       b( i, j ) = ival( i, j, imb )
                                       IF( abs( i-j ).LE.1 ) THEN
                                          b( i, j ) = b( i, j )*
     $                                                vm1( imldb1 )
                                       ELSE
                                          b( i, j ) = b( i, j )*
     $                                                vm1( imloff )
                                       END IF
                                       tnrm = max( tnrm,
     $                                        abs( b( i, j ) ) )
   30                               CONTINUE
   40                            CONTINUE
                                 cnrm = zero
                                 DO 60 i = 1, m
                                    DO 50 j = 1, n
                                       c( i, j ) = sin( real( i*j ) )
                                       cnrm = max( cnrm, c( i, j ) )
                                       cc( i, j ) = c( i, j )
   50                               CONTINUE
   60                            CONTINUE
                                 knt = knt + 1
                                 CALL strsyl( trana, tranb, isgn, m, n,
     $                                        a, 6, b, 6, c, 6, scale,
     $                                        info )
                                 IF( info.NE.0 )
     $                              ninfo = ninfo + 1
                                 xnrm = slange( 'M', m, n, c, 6, dum )
                                 rmul = one
                                 IF( xnrm.GT.one .AND. tnrm.GT.one )
     $                                THEN
                                    IF( xnrm.GT.bignum / tnrm ) THEN
                                       rmul = one / max( xnrm, tnrm )
                                    END IF
                                 END IF
                                 CALL sgemm( trana, 'N', m, n, m, rmul,
     $                                       a, 6, c, 6, -scale*rmul,
     $                                       cc, 6 )
                                 CALL sgemm( 'N', tranb, m, n, n,
     $                                       real( isgn )*rmul, c, 6, b,
     $                                       6, one, cc, 6 )
                                 res1 = slange( 'M', m, n, cc, 6, dum )
                                 res = res1 / max( smlnum, smlnum*xnrm,
     $                                 ( ( rmul*tnrm )*eps )*xnrm )
                                 IF( res.GT.rmax ) THEN
                                    lmax = knt
                                    rmax = res
                                 END IF
   70                         CONTINUE
   80                      CONTINUE
   90                   CONTINUE
  100                CONTINUE
  110             CONTINUE
  120          CONTINUE
  130       CONTINUE
  140    CONTINUE
  150 CONTINUE
*
      RETURN
*
*     End of SGET35
*
      END