LAPACK: TESTING/EIG/zget35.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZGET35( RMAX, LMAX, NINFO, KNT, NIN )
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, NIN, NINFO
00009       DOUBLE PRECISION   RMAX
00010 *     ..
00011 *
00012 *  Purpose
00013 *  =======
00014 *
00015 *  ZGET35 tests ZTRSYL, 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 *  NIN     (input) INTEGER
00045 *          Input logical unit number.
00046 *
00047 *  =====================================================================
00048 *
00049 *     .. Parameters ..
00050       INTEGER            LDT
00051       PARAMETER          ( LDT = 10 )
00052       DOUBLE PRECISION   ZERO, ONE, TWO
00053       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
00054       DOUBLE PRECISION   LARGE
00055       PARAMETER          ( LARGE = 1.0D6 )
00056       COMPLEX*16         CONE
00057       PARAMETER          ( CONE = 1.0D0 )
00058 *     ..
00059 *     .. Local Scalars ..
00060       CHARACTER          TRANA, TRANB
00061       INTEGER            I, IMLA, IMLAD, IMLB, IMLC, INFO, ISGN, ITRANA,
00062      $                   ITRANB, J, M, N
00063       DOUBLE PRECISION   BIGNUM, EPS, RES, RES1, SCALE, SMLNUM, TNRM,
00064      $                   XNRM
00065       COMPLEX*16         RMUL
00066 *     ..
00067 *     .. Local Arrays ..
00068       DOUBLE PRECISION   DUM( 1 ), VM1( 3 ), VM2( 3 )
00069       COMPLEX*16         A( LDT, LDT ), ATMP( LDT, LDT ), B( LDT, LDT ),
00070      $                   BTMP( LDT, LDT ), C( LDT, LDT ),
00071      $                   CSAV( LDT, LDT ), CTMP( LDT, LDT )
00072 *     ..
00073 *     .. External Functions ..
00074       DOUBLE PRECISION   DLAMCH, ZLANGE
00075       EXTERNAL           DLAMCH, ZLANGE
00076 *     ..
00077 *     .. External Subroutines ..
00078       EXTERNAL           DLABAD, ZGEMM, ZTRSYL
00079 *     ..
00080 *     .. Intrinsic Functions ..
00081       INTRINSIC          ABS, DBLE, MAX, SQRT
00082 *     ..
00083 *     .. Executable Statements ..
00084 *
00085 *     Get machine parameters
00086 *
00087       EPS = DLAMCH( 'P' )
00088       SMLNUM = DLAMCH( 'S' ) / EPS
00089       BIGNUM = ONE / SMLNUM
00090       CALL DLABAD( SMLNUM, BIGNUM )
00091 *
00092 *     Set up test case parameters
00093 *
00094       VM1( 1 ) = SQRT( SMLNUM )
00095       VM1( 2 ) = ONE
00096       VM1( 3 ) = LARGE
00097       VM2( 1 ) = ONE
00098       VM2( 2 ) = ONE + TWO*EPS
00099       VM2( 3 ) = TWO
00100 *
00101       KNT = 0
00102       NINFO = 0
00103       LMAX = 0
00104       RMAX = ZERO
00105 *
00106 *     Begin test loop
00107 *
00108    10 CONTINUE
00109       READ( NIN, FMT = * )M, N
00110       IF( N.EQ.0 )
00111      $   RETURN
00112       DO 20 I = 1, M
00113          READ( NIN, FMT = * )( ATMP( I, J ), J = 1, M )
00114    20 CONTINUE
00115       DO 30 I = 1, N
00116          READ( NIN, FMT = * )( BTMP( I, J ), J = 1, N )
00117    30 CONTINUE
00118       DO 40 I = 1, M
00119          READ( NIN, FMT = * )( CTMP( I, J ), J = 1, N )
00120    40 CONTINUE
00121       DO 170 IMLA = 1, 3
00122          DO 160 IMLAD = 1, 3
00123             DO 150 IMLB = 1, 3
00124                DO 140 IMLC = 1, 3
00125                   DO 130 ITRANA = 1, 2
00126                      DO 120 ITRANB = 1, 2
00127                         DO 110 ISGN = -1, 1, 2
00128                            IF( ITRANA.EQ.1 )
00129      $                        TRANA = 'N'
00130                            IF( ITRANA.EQ.2 )
00131      $                        TRANA = 'C'
00132                            IF( ITRANB.EQ.1 )
00133      $                        TRANB = 'N'
00134                            IF( ITRANB.EQ.2 )
00135      $                        TRANB = 'C'
00136                            TNRM = ZERO
00137                            DO 60 I = 1, M
00138                               DO 50 J = 1, M
00139                                  A( I, J ) = ATMP( I, J )*VM1( IMLA )
00140                                  TNRM = MAX( TNRM, ABS( A( I, J ) ) )
00141    50                         CONTINUE
00142                               A( I, I ) = A( I, I )*VM2( IMLAD )
00143                               TNRM = MAX( TNRM, ABS( A( I, I ) ) )
00144    60                      CONTINUE
00145                            DO 80 I = 1, N
00146                               DO 70 J = 1, N
00147                                  B( I, J ) = BTMP( I, J )*VM1( IMLB )
00148                                  TNRM = MAX( TNRM, ABS( B( I, J ) ) )
00149    70                         CONTINUE
00150    80                      CONTINUE
00151                            IF( TNRM.EQ.ZERO )
00152      $                        TNRM = ONE
00153                            DO 100 I = 1, M
00154                               DO 90 J = 1, N
00155                                  C( I, J ) = CTMP( I, J )*VM1( IMLC )
00156                                  CSAV( I, J ) = C( I, J )
00157    90                         CONTINUE
00158   100                      CONTINUE
00159                            KNT = KNT + 1
00160                            CALL ZTRSYL( TRANA, TRANB, ISGN, M, N, A,
00161      $                                  LDT, B, LDT, C, LDT, SCALE,
00162      $                                  INFO )
00163                            IF( INFO.NE.0 )
00164      $                        NINFO = NINFO + 1
00165                            XNRM = ZLANGE( 'M', M, N, C, LDT, DUM )
00166                            RMUL = CONE
00167                            IF( XNRM.GT.ONE .AND. TNRM.GT.ONE ) THEN
00168                               IF( XNRM.GT.BIGNUM / TNRM ) THEN
00169                                  RMUL = MAX( XNRM, TNRM )
00170                                  RMUL = CONE / RMUL
00171                               END IF
00172                            END IF
00173                            CALL ZGEMM( TRANA, 'N', M, N, M, RMUL, A,
00174      $                                 LDT, C, LDT, -SCALE*RMUL, CSAV,
00175      $                                 LDT )
00176                            CALL ZGEMM( 'N', TRANB, M, N, N,
00177      $                                 DBLE( ISGN )*RMUL, C, LDT, B,
00178      $                                 LDT, CONE, CSAV, LDT )
00179                            RES1 = ZLANGE( 'M', M, N, CSAV, LDT, DUM )
00180                            RES = RES1 / MAX( SMLNUM, SMLNUM*XNRM,
00181      $                           ( ( ABS( RMUL )*TNRM )*EPS )*XNRM )
00182                            IF( RES.GT.RMAX ) THEN
00183                               LMAX = KNT
00184                               RMAX = RES
00185                            END IF
00186   110                   CONTINUE
00187   120                CONTINUE
00188   130             CONTINUE
00189   140          CONTINUE
00190   150       CONTINUE
00191   160    CONTINUE
00192   170 CONTINUE
00193       GO TO 10
00194 *
00195 *     End of ZGET35
00196 *
00197       END