1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
|
SUBROUTINE DGET35( RMAX, LMAX, NINFO, KNT )
*
* -- LAPACK test routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* February 29, 1992
*
* .. Scalar Arguments ..
INTEGER KNT, LMAX, NINFO
DOUBLE PRECISION RMAX
* ..
*
* Purpose
* =======
*
* DGET35 tests DTRSYL, 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))
*
* Arguments
* ==========
*
* RMAX (output) DOUBLE PRECISION
* Value of the largest test ratio.
*
* LMAX (output) INTEGER
* Example number where largest test ratio achieved.
*
* NINFO (output) INTEGER
* Number of examples where INFO is nonzero.
*
* KNT (output) INTEGER
* Total number of examples tested.
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
DOUBLE PRECISION TWO, FOUR
PARAMETER ( TWO = 2.0D0, FOUR = 4.0D0 )
* ..
* .. Local Scalars ..
CHARACTER TRANA, TRANB
INTEGER I, IMA, IMB, IMLDA1, IMLDA2, IMLDB1, IMLOFF,
$ INFO, ISGN, ITRANA, ITRANB, J, M, N
DOUBLE PRECISION BIGNUM, CNRM, EPS, RES, RES1, RMUL, SCALE,
$ SMLNUM, TNRM, XNRM
* ..
* .. Local Arrays ..
INTEGER IDIM( 8 ), IVAL( 6, 6, 8 )
DOUBLE PRECISION A( 6, 6 ), B( 6, 6 ), C( 6, 6 ), CC( 6, 6 ),
$ DUM( 1 ), VM1( 3 ), VM2( 3 )
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH, DLANGE
EXTERNAL DLAMCH, DLANGE
* ..
* .. External Subroutines ..
EXTERNAL DGEMM, DLABAD, DTRSYL
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, MAX, 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 = DLAMCH( 'P' )
SMLNUM = DLAMCH( 'S' )*FOUR / EPS
BIGNUM = ONE / SMLNUM
CALL DLABAD( 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( DBLE( I*J ) )
CNRM = MAX( CNRM, C( I, J ) )
CC( I, J ) = C( I, J )
50 CONTINUE
60 CONTINUE
KNT = KNT + 1
CALL DTRSYL( TRANA, TRANB, ISGN, M, N,
$ A, 6, B, 6, C, 6, SCALE,
$ INFO )
IF( INFO.NE.0 )
$ NINFO = NINFO + 1
XNRM = DLANGE( '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 DGEMM( TRANA, 'N', M, N, M, RMUL,
$ A, 6, C, 6, -SCALE*RMUL,
$ CC, 6 )
CALL DGEMM( 'N', TRANB, M, N, N,
$ DBLE( ISGN )*RMUL, C, 6, B,
$ 6, ONE, CC, 6 )
RES1 = DLANGE( '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 DGET35
*
END
|
