LAPACK: TESTING/EIG/dget35.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DGET35( RMAX, LMAX, NINFO, KNT )
00002 *
00003 *  -- LAPACK test routine (version 3.1) --
00004 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00005 *     November 2006
00006 *
00007 *     .. Scalar Arguments ..
00008       INTEGER            KNT, LMAX, NINFO
00009       DOUBLE PRECISION   RMAX
00010 *     ..
00011 *
00012 *  Purpose
00013 *  =======
00014 *
00015 *  DGET35 tests DTRSYL, a routine for solving the Sylvester matrix
00016 *  equation
00017 *
00018 *     op(A)*X + ISGN*X*op(B) = scale*C,
00019 *
00020 *  A and B are assumed to be in Schur canonical form, op() represents an
00021 *  optional transpose, and ISGN can be -1 or +1.  Scale is an output
00022 *  less than or equal to 1, chosen to avoid overflow in X.
00023 *
00024 *  The test code verifies that the following residual is order 1:
00025 *
00026 *     norm(op(A)*X + ISGN*X*op(B) - scale*C) /
00027 *         (EPS*max(norm(A),norm(B))*norm(X))
00028 *
00029 *  Arguments
00030 *  ==========
00031 *
00032 *  RMAX    (output) DOUBLE PRECISION
00033 *          Value of the largest test ratio.
00034 *
00035 *  LMAX    (output) INTEGER
00036 *          Example number where largest test ratio achieved.
00037 *
00038 *  NINFO   (output) INTEGER
00039 *          Number of examples where INFO is nonzero.
00040 *
00041 *  KNT     (output) INTEGER
00042 *          Total number of examples tested.
00043 *
00044 *  =====================================================================
00045 *
00046 *     .. Parameters ..
00047       DOUBLE PRECISION   ZERO, ONE
00048       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00049       DOUBLE PRECISION   TWO, FOUR
00050       PARAMETER          ( TWO = 2.0D0, FOUR = 4.0D0 )
00051 *     ..
00052 *     .. Local Scalars ..
00053       CHARACTER          TRANA, TRANB
00054       INTEGER            I, IMA, IMB, IMLDA1, IMLDA2, IMLDB1, IMLOFF,
00055      $                   INFO, ISGN, ITRANA, ITRANB, J, M, N
00056       DOUBLE PRECISION   BIGNUM, CNRM, EPS, RES, RES1, RMUL, SCALE,
00057      $                   SMLNUM, TNRM, XNRM
00058 *     ..
00059 *     .. Local Arrays ..
00060       INTEGER            IDIM( 8 ), IVAL( 6, 6, 8 )
00061       DOUBLE PRECISION   A( 6, 6 ), B( 6, 6 ), C( 6, 6 ), CC( 6, 6 ),
00062      $                   DUM( 1 ), VM1( 3 ), VM2( 3 )
00063 *     ..
00064 *     .. External Functions ..
00065       DOUBLE PRECISION   DLAMCH, DLANGE
00066       EXTERNAL           DLAMCH, DLANGE
00067 *     ..
00068 *     .. External Subroutines ..
00069       EXTERNAL           DGEMM, DLABAD, DTRSYL
00070 *     ..
00071 *     .. Intrinsic Functions ..
00072       INTRINSIC          ABS, DBLE, MAX, SIN, SQRT
00073 *     ..
00074 *     .. Data statements ..
00075       DATA               IDIM / 1, 2, 3, 4, 3, 3, 6, 4 /
00076       DATA               IVAL / 1, 35*0, 1, 2, 4*0, -2, 0, 28*0, 1, 5*0,
00077      $                   5, 1, 2, 3*0, -8, -2, 1, 21*0, 3, 4, 4*0, -5,
00078      $                   3, 4*0, 1, 2, 1, 4, 2*0, -3, -9, -1, 1, 14*0,
00079      $                   1, 5*0, 2, 3, 4*0, 5, 6, 7, 21*0, 1, 5*0, 1, 3,
00080      $                   -4, 3*0, 2, 5, 2, 21*0, 1, 2, 4*0, -2, 0, 4*0,
00081      $                   5, 6, 3, 4, 2*0, -1, -9, -5, 2, 2*0, 4*8, 5, 6,
00082      $                   4*9, -7, 5, 1, 5*0, 1, 5, 2, 3*0, 2, -21, 5,
00083      $                   3*0, 1, 2, 3, 4, 14*0 /
00084 *     ..
00085 *     .. Executable Statements ..
00086 *
00087 *     Get machine parameters
00088 *
00089       EPS = DLAMCH( 'P' )
00090       SMLNUM = DLAMCH( 'S' )*FOUR / EPS
00091       BIGNUM = ONE / SMLNUM
00092       CALL DLABAD( SMLNUM, BIGNUM )
00093 *
00094 *     Set up test case parameters
00095 *
00096       VM1( 1 ) = SQRT( SMLNUM )
00097       VM1( 2 ) = ONE
00098       VM1( 3 ) = SQRT( BIGNUM )
00099       VM2( 1 ) = ONE
00100       VM2( 2 ) = ONE + TWO*EPS
00101       VM2( 3 ) = TWO
00102 *
00103       KNT = 0
00104       NINFO = 0
00105       LMAX = 0
00106       RMAX = ZERO
00107 *
00108 *     Begin test loop
00109 *
00110       DO 150 ITRANA = 1, 2
00111          DO 140 ITRANB = 1, 2
00112             DO 130 ISGN = -1, 1, 2
00113                DO 120 IMA = 1, 8
00114                   DO 110 IMLDA1 = 1, 3
00115                      DO 100 IMLDA2 = 1, 3
00116                         DO 90 IMLOFF = 1, 2
00117                            DO 80 IMB = 1, 8
00118                               DO 70 IMLDB1 = 1, 3
00119                                  IF( ITRANA.EQ.1 )
00120      $                              TRANA = 'N'
00121                                  IF( ITRANA.EQ.2 )
00122      $                              TRANA = 'T'
00123                                  IF( ITRANB.EQ.1 )
00124      $                              TRANB = 'N'
00125                                  IF( ITRANB.EQ.2 )
00126      $                              TRANB = 'T'
00127                                  M = IDIM( IMA )
00128                                  N = IDIM( IMB )
00129                                  TNRM = ZERO
00130                                  DO 20 I = 1, M
00131                                     DO 10 J = 1, M
00132                                        A( I, J ) = IVAL( I, J, IMA )
00133                                        IF( ABS( I-J ).LE.1 ) THEN
00134                                           A( I, J ) = A( I, J )*
00135      $                                                VM1( IMLDA1 )
00136                                           A( I, J ) = A( I, J )*
00137      $                                                VM2( IMLDA2 )
00138                                        ELSE
00139                                           A( I, J ) = A( I, J )*
00140      $                                                VM1( IMLOFF )
00141                                        END IF
00142                                        TNRM = MAX( TNRM,
00143      $                                        ABS( A( I, J ) ) )
00144    10                               CONTINUE
00145    20                            CONTINUE
00146                                  DO 40 I = 1, N
00147                                     DO 30 J = 1, N
00148                                        B( I, J ) = IVAL( I, J, IMB )
00149                                        IF( ABS( I-J ).LE.1 ) THEN
00150                                           B( I, J ) = B( I, J )*
00151      $                                                VM1( IMLDB1 )
00152                                        ELSE
00153                                           B( I, J ) = B( I, J )*
00154      $                                                VM1( IMLOFF )
00155                                        END IF
00156                                        TNRM = MAX( TNRM,
00157      $                                        ABS( B( I, J ) ) )
00158    30                               CONTINUE
00159    40                            CONTINUE
00160                                  CNRM = ZERO
00161                                  DO 60 I = 1, M
00162                                     DO 50 J = 1, N
00163                                        C( I, J ) = SIN( DBLE( I*J ) )
00164                                        CNRM = MAX( CNRM, C( I, J ) )
00165                                        CC( I, J ) = C( I, J )
00166    50                               CONTINUE
00167    60                            CONTINUE
00168                                  KNT = KNT + 1
00169                                  CALL DTRSYL( TRANA, TRANB, ISGN, M, N,
00170      $                                        A, 6, B, 6, C, 6, SCALE,
00171      $                                        INFO )
00172                                  IF( INFO.NE.0 )
00173      $                              NINFO = NINFO + 1
00174                                  XNRM = DLANGE( 'M', M, N, C, 6, DUM )
00175                                  RMUL = ONE
00176                                  IF( XNRM.GT.ONE .AND. TNRM.GT.ONE )
00177      $                                THEN
00178                                     IF( XNRM.GT.BIGNUM / TNRM ) THEN
00179                                        RMUL = ONE / MAX( XNRM, TNRM )
00180                                     END IF
00181                                  END IF
00182                                  CALL DGEMM( TRANA, 'N', M, N, M, RMUL,
00183      $                                       A, 6, C, 6, -SCALE*RMUL,
00184      $                                       CC, 6 )
00185                                  CALL DGEMM( 'N', TRANB, M, N, N,
00186      $                                       DBLE( ISGN )*RMUL, C, 6, B,
00187      $                                       6, ONE, CC, 6 )
00188                                  RES1 = DLANGE( 'M', M, N, CC, 6, DUM )
00189                                  RES = RES1 / MAX( SMLNUM, SMLNUM*XNRM,
00190      $                                 ( ( RMUL*TNRM )*EPS )*XNRM )
00191                                  IF( RES.GT.RMAX ) THEN
00192                                     LMAX = KNT
00193                                     RMAX = RES
00194                                  END IF
00195    70                         CONTINUE
00196    80                      CONTINUE
00197    90                   CONTINUE
00198   100                CONTINUE
00199   110             CONTINUE
00200   120          CONTINUE
00201   130       CONTINUE
00202   140    CONTINUE
00203   150 CONTINUE
00204 *
00205       RETURN
00206 *
00207 *     End of DGET35
00208 *
00209       END