ssyl01.f

Go to the documentation of this file.

*> \brief \b SSYL01
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE SSYL01( THRESH, NFAIL, RMAX, NINFO, KNT )
*
*      .. Scalar Arguments ..
*      INTEGER            KNT
*      REAL               THRESH
*      ..
*      .. Array Arguments ..
*      INTEGER            NFAIL( 3 ), NINFO( 2 )
*      REAL               RMAX( 2 )
*      ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SSYL01 tests STRSYL and STRSYL3, routines 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 does not exceed
*> the provided threshold:
*>
*>    norm(op(A)*X + ISGN*X*op(B) - scale*C) /
*>        (EPS*max(norm(A),norm(B))*norm(X))
*>
*> This routine complements SGET35 by testing with larger,
*> random matrices, of which some require rescaling of X to avoid overflow.
*>
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] THRESH
*> \verbatim
*>          THRESH is REAL
*>          A test will count as "failed" if the residual, computed as
*>          described above, exceeds THRESH.
*> \endverbatim
*>
*> \param[out] NFAIL
*> \verbatim
*>          NFAIL is INTEGER array, dimension (3)
*>          NFAIL(1) = No. of times residual STRSYL exceeds threshold THRESH
*>          NFAIL(2) = No. of times residual STRSYL3 exceeds threshold THRESH
*>          NFAIL(3) = No. of times STRSYL3 and STRSYL deviate
*> \endverbatim
*>
*> \param[out] RMAX
*> \verbatim
*>          RMAX is REAL, dimension (2)
*>          RMAX(1) = Value of the largest test ratio of STRSYL
*>          RMAX(2) = Value of the largest test ratio of STRSYL3
*> \endverbatim
*>
*> \param[out] NINFO
*> \verbatim
*>          NINFO is INTEGER array, dimension (2)
*>          NINFO(1) = No. of times STRSYL returns an expected INFO
*>          NINFO(2) = No. of times STRSYL3 returns an expected INFO
*> \endverbatim
*>
*> \param[out] KNT
*> \verbatim
*>          KNT is INTEGER
*>          Total number of examples tested.
*> \endverbatim
 
*
*  -- LAPACK test routine --
      SUBROUTINE ssyl01( THRESH, NFAIL, RMAX, NINFO, KNT )
      IMPLICIT NONE
*
*  -- 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
      REAL               THRESH
*     ..
*     .. Array Arguments ..
      INTEGER            NFAIL( 3 ), NINFO( 2 )
      REAL               RMAX( 2 )
*     ..
*
*  =====================================================================
*     ..
*     .. Parameters ..
      REAL               ZERO, ONE
      parameter( zero = 0.0e+0, one = 1.0e+0 )
      INTEGER            MAXM, MAXN, LDSWORK
      parameter( maxm = 101, maxn = 138, ldswork = 18 )
*     ..
*     .. Local Scalars ..
      CHARACTER          TRANA, TRANB
      INTEGER            I, INFO, IINFO, ISGN, ITRANA, ITRANB, J, KLA,
     $                   KUA, KLB, KUB, LIWORK, M, N
      REAL               ANRM, BNRM, BIGNUM, EPS, RES, RES1, RMUL,
     $                   SCALE, SCALE3, SMLNUM, TNRM, XNRM
*     ..
*     .. Local Arrays ..
      REAL               A( MAXM, MAXM ), B( MAXN, MAXN ),
     $                   C( MAXM, MAXN ), CC( MAXM, MAXN ),
     $                   X( MAXM, MAXN ),
     $                   DUML( MAXM ), DUMR( MAXN ),
     $                   D( MAX( MAXM, MAXN ) ), DUM( MAXN ),
     $                   SWORK( LDSWORK, 54 ), VM( 2 )
      INTEGER            ISEED( 4 ), IWORK( MAXM + MAXN + 2 ), IDUM( 2 )
*     ..
*     .. External Functions ..
      LOGICAL            SISNAN
      REAL               SLAMCH, SLANGE
      EXTERNAL           sisnan, slamch, slange
*     ..
*     .. External Subroutines ..
      EXTERNAL           slatmr, slacpy, sgemm, strsyl, strsyl3
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, real, max
*     ..
*     .. Executable Statements ..
*
*     Get machine parameters
*
      eps = slamch( 'P' )
      smlnum = slamch( 'S' ) / eps
      bignum = one / smlnum
*
      vm( 1 ) = one
      vm( 2 ) = 0.05e+0
*
*     Begin test loop
*
      ninfo( 1 ) = 0
      ninfo( 2 ) = 0
      nfail( 1 ) = 0
      nfail( 2 ) = 0
      nfail( 3 ) = 0
      rmax( 1 ) = zero
      rmax( 2 ) = zero
      knt = 0
      DO i = 1, 4
         iseed( i ) = 1
      END DO
      scale = one
      scale3 = one
      liwork = maxm + maxn + 2
      DO j = 1, 2
         DO isgn = -1, 1, 2
*           Reset seed (overwritten by LATMR)
            DO i = 1, 4
               iseed( i ) = 1
            END DO
            DO m = 32, maxm, 71
               kla = 0
               kua = m - 1
               CALL slatmr( m, m, 'S', iseed, 'N', d,
     $                      6, one, one, 'T', 'N',
     $                      duml, 1, one, dumr, 1, one,
     $                      'N', iwork, kla, kua, zero,
     $                      one, 'NO', a, maxm, iwork, iinfo )
               DO i = 1, m
                  a( i, i ) = a( i, i ) * vm( j )
               END DO
               anrm = slange( 'M', m, m, a, maxm, dum )
               DO n = 51, maxn, 47
                  klb = 0
                  kub = n - 1
                  CALL slatmr( n, n, 'S', iseed, 'N', d,
     $                         6, one, one, 'T', 'N',
     $                         duml, 1, one, dumr, 1, one,
     $                         'N', iwork, klb, kub, zero,
     $                         one, 'NO', b, maxn, iwork, iinfo )
                  bnrm = slange( 'M', n, n, b, maxn, dum )
                  tnrm = max( anrm, bnrm )
                  CALL slatmr( m, n, 'S', iseed, 'N', d,
     $                         6, one, one, 'T', 'N',
     $                         duml, 1, one, dumr, 1, one,
     $                         'N', iwork, m, n, zero, one,
     $                         'NO', c, maxm, iwork, iinfo )
                  DO itrana = 1, 2
                     IF( itrana.EQ.1 ) THEN
                        trana = 'N'
                     END IF
                     IF( itrana.EQ.2 ) THEN
                        trana = 'T'
                     END IF
                     DO itranb = 1, 2
                        IF( itranb.EQ.1 ) THEN
                           tranb = 'N'
                        END IF
                        IF( itranb.EQ.2 ) THEN
                           tranb = 'T'
                        END IF
                        knt = knt + 1
*
                        CALL slacpy( 'All', m, n, c, maxm, x, maxm)
                        CALL slacpy( 'All', m, n, c, maxm, cc, maxm)
                        CALL strsyl( trana, tranb, isgn, m, n, 
     $                               a, maxm, b, maxn, x, maxm,
     $                               scale, iinfo )
                        IF( iinfo.NE.0 )
     $                     ninfo( 1 ) = ninfo( 1 ) + 1
                        xnrm = slange( 'M', m, n, x, maxm, 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, maxm, x, maxm, -scale*rmul,
     $                              c, maxm )
                        CALL sgemm( 'N', tranb, m, n, n,
     $                               real( isgn )*rmul, x, maxm, b,
     $                               maxn, one, c, maxm )
                        res1 = slange( 'M', m, n, c, maxm, dum )
                        res = res1 / max( smlnum, smlnum*xnrm,
     $                              ( ( rmul*tnrm )*eps )*xnrm )
                        IF( res.GT.thresh )
     $                     nfail( 1 ) = nfail( 1 ) + 1
                        IF( res.GT.rmax( 1 ) )
     $                     rmax( 1 ) = res
*
                        CALL slacpy( 'All', m, n, c, maxm, x, maxm )
                        CALL slacpy( 'All', m, n, c, maxm, cc, maxm )
                        CALL strsyl3( trana, tranb, isgn, m, n,
     $                                a, maxm, b, maxn, x, maxm,
     $                                scale3, iwork, liwork,
     $                                swork, ldswork, info)
                        IF( info.NE.0 )
     $                     ninfo( 2 ) = ninfo( 2 ) + 1
                        xnrm = slange( 'M', m, n, x, maxm, 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, maxm, x, maxm, -scale3*rmul,
     $                              cc, maxm )
                        CALL sgemm( 'N', tranb, m, n, n,
     $                              real( isgn )*rmul, x, maxm, b,
     $                              maxn, one, cc, maxm )
                        res1 = slange( 'M', m, n, cc, maxm, dum )
                        res = res1 / max( smlnum, smlnum*xnrm,
     $                             ( ( rmul*tnrm )*eps )*xnrm )
*                       Verify that TRSYL3 only flushes if TRSYL flushes (but
*                       there may be cases where TRSYL3 avoid flushing).
                        IF( scale3.EQ.zero .AND. scale.GT.zero .OR. 
     $                      iinfo.NE.info ) THEN
                           nfail( 3 ) = nfail( 3 ) + 1
                        END IF
                        IF( res.GT.thresh .OR. sisnan( res ) )
     $                     nfail( 2 ) = nfail( 2 ) + 1
                        IF( res.GT.rmax( 2 ) )
     $                     rmax( 2 ) = res
                     END DO
                  END DO
               END DO
            END DO
         END DO
      END DO
*
      RETURN
*
*     End of SSYL01
*
      END