zget35.f

Go to the documentation of this file.

*> \brief \b ZGET35
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at 
*            http://www.netlib.org/lapack/explore-html/ 
*
*  Definition:
*  ===========
*
*       SUBROUTINE ZGET35( RMAX, LMAX, NINFO, KNT, NIN )
* 
*       .. Scalar Arguments ..
*       INTEGER            KNT, LMAX, NIN, NINFO
*       DOUBLE PRECISION   RMAX
*       ..
*  
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> ZGET35 tests ZTRSYL, 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 DOUBLE PRECISION
*>          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
*>
*> \param[in] NIN
*> \verbatim
*>          NIN is INTEGER
*>          Input logical unit number.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee 
*> \author Univ. of California Berkeley 
*> \author Univ. of Colorado Denver 
*> \author NAG Ltd. 
*
*> \date November 2011
*
*> \ingroup complex16_eig
*
*  =====================================================================
      SUBROUTINE zget35( RMAX, LMAX, NINFO, KNT, NIN )
*
*  -- LAPACK test routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      INTEGER            knt, lmax, nin, ninfo
      DOUBLE PRECISION   rmax
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      INTEGER            ldt
      parameter( ldt = 10 )
      DOUBLE PRECISION   zero, one, two
      parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
      DOUBLE PRECISION   large
      parameter( large = 1.0d6 )
      COMPLEX*16         cone
      parameter( cone = 1.0d0 )
*     ..
*     .. Local Scalars ..
      CHARACTER          trana, tranb
      INTEGER            i, imla, imlad, imlb, imlc, info, isgn, itrana,
     $                   itranb, j, m, n
      DOUBLE PRECISION   bignum, eps, res, res1, scale, smlnum, tnrm,
     $                   xnrm
      COMPLEX*16         rmul
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   dum( 1 ), vm1( 3 ), vm2( 3 )
      COMPLEX*16         a( ldt, ldt ), atmp( ldt, ldt ), b( ldt, ldt ),
     $                   btmp( ldt, ldt ), c( ldt, ldt ),
     $                   csav( ldt, ldt ), ctmp( ldt, ldt )
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   dlamch, zlange
      EXTERNAL           dlamch, zlange
*     ..
*     .. External Subroutines ..
      EXTERNAL           dlabad, zgemm, ztrsyl
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, max, sqrt
*     ..
*     .. Executable Statements ..
*
*     Get machine parameters
*
      eps = dlamch( 'P' )
      smlnum = dlamch( 'S' ) / eps
      bignum = one / smlnum
      CALL dlabad( smlnum, bignum )
*
*     Set up test case parameters
*
      vm1( 1 ) = sqrt( smlnum )
      vm1( 2 ) = one
      vm1( 3 ) = large
      vm2( 1 ) = one
      vm2( 2 ) = one + two*eps
      vm2( 3 ) = two
*
      knt = 0
      ninfo = 0
      lmax = 0
      rmax = zero
*
*     Begin test loop
*
   10 continue
      READ( nin, fmt = * )m, n
      IF( n.EQ.0 )
     $   return
      DO 20 i = 1, m
         READ( nin, fmt = * )( atmp( i, j ), j = 1, m )
   20 continue
      DO 30 i = 1, n
         READ( nin, fmt = * )( btmp( i, j ), j = 1, n )
   30 continue
      DO 40 i = 1, m
         READ( nin, fmt = * )( ctmp( i, j ), j = 1, n )
   40 continue
      DO 170 imla = 1, 3
         DO 160 imlad = 1, 3
            DO 150 imlb = 1, 3
               DO 140 imlc = 1, 3
                  DO 130 itrana = 1, 2
                     DO 120 itranb = 1, 2
                        DO 110 isgn = -1, 1, 2
                           IF( itrana.EQ.1 )
     $                        trana = 'N'
                           IF( itrana.EQ.2 )
     $                        trana = 'C'
                           IF( itranb.EQ.1 )
     $                        tranb = 'N'
                           IF( itranb.EQ.2 )
     $                        tranb = 'C'
                           tnrm = zero
                           DO 60 i = 1, m
                              DO 50 j = 1, m
                                 a( i, j ) = atmp( i, j )*vm1( imla )
                                 tnrm = max( tnrm, abs( a( i, j ) ) )
   50                         continue
                              a( i, i ) = a( i, i )*vm2( imlad )
                              tnrm = max( tnrm, abs( a( i, i ) ) )
   60                      continue
                           DO 80 i = 1, n
                              DO 70 j = 1, n
                                 b( i, j ) = btmp( i, j )*vm1( imlb )
                                 tnrm = max( tnrm, abs( b( i, j ) ) )
   70                         continue
   80                      continue
                           IF( tnrm.EQ.zero )
     $                        tnrm = one
                           DO 100 i = 1, m
                              DO 90 j = 1, n
                                 c( i, j ) = ctmp( i, j )*vm1( imlc )
                                 csav( i, j ) = c( i, j )
   90                         continue
  100                      continue
                           knt = knt + 1
                           CALL ztrsyl( trana, tranb, isgn, m, n, a,
     $                                  ldt, b, ldt, c, ldt, scale,
     $                                  info )
                           IF( info.NE.0 )
     $                        ninfo = ninfo + 1
                           xnrm = zlange( 'M', m, n, c, ldt, dum )
                           rmul = cone
                           IF( xnrm.GT.one .AND. tnrm.GT.one ) THEN
                              IF( xnrm.GT.bignum / tnrm ) THEN
                                 rmul = max( xnrm, tnrm )
                                 rmul = cone / rmul
                              END IF
                           END IF
                           CALL zgemm( trana, 'N', m, n, m, rmul, a,
     $                                 ldt, c, ldt, -scale*rmul, csav,
     $                                 ldt )
                           CALL zgemm( 'N', tranb, m, n, n,
     $                                 dble( isgn )*rmul, c, ldt, b,
     $                                 ldt, cone, csav, ldt )
                           res1 = zlange( 'M', m, n, csav, ldt, dum )
                           res = res1 / max( smlnum, smlnum*xnrm,
     $                           ( ( abs( rmul )*tnrm )*eps )*xnrm )
                           IF( res.GT.rmax ) THEN
                              lmax = knt
                              rmax = res
                           END IF
  110                   continue
  120                continue
  130             continue
  140          continue
  150       continue
  160    continue
  170 continue
      go to 10
*
*     End of ZGET35
*
      END