LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ssyl01.f
Go to the documentation of this file.
1*> \brief \b SSYL01
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE SSYL01( THRESH, NFAIL, RMAX, NINFO, KNT )
12*
13* .. Scalar Arguments ..
14* INTEGER KNT
15* REAL THRESH
16* ..
17* .. Array Arguments ..
18* INTEGER NFAIL( 3 ), NINFO( 2 )
19* REAL RMAX( 2 )
20* ..
21*
22*
23*> \par Purpose:
24* =============
25*>
26*> \verbatim
27*>
28*> SSYL01 tests STRSYL and STRSYL3, routines for solving the Sylvester matrix
29*> equation
30*>
31*> op(A)*X + ISGN*X*op(B) = scale*C,
32*>
33*> A and B are assumed to be in Schur canonical form, op() represents an
34*> optional transpose, and ISGN can be -1 or +1. Scale is an output
35*> less than or equal to 1, chosen to avoid overflow in X.
36*>
37*> The test code verifies that the following residual does not exceed
38*> the provided threshold:
39*>
40*> norm(op(A)*X + ISGN*X*op(B) - scale*C) /
41*> (EPS*max(norm(A),norm(B))*norm(X))
42*>
43*> This routine complements SGET35 by testing with larger,
44*> random matrices, of which some require rescaling of X to avoid overflow.
45*>
46*> \endverbatim
47*
48* Arguments:
49* ==========
50*
51*> \param[in] THRESH
52*> \verbatim
53*> THRESH is REAL
54*> A test will count as "failed" if the residual, computed as
55*> described above, exceeds THRESH.
56*> \endverbatim
57*>
58*> \param[out] NFAIL
59*> \verbatim
60*> NFAIL is INTEGER array, dimension (3)
61*> NFAIL(1) = No. of times residual STRSYL exceeds threshold THRESH
62*> NFAIL(2) = No. of times residual STRSYL3 exceeds threshold THRESH
63*> NFAIL(3) = No. of times STRSYL3 and STRSYL deviate
64*> \endverbatim
65*>
66*> \param[out] RMAX
67*> \verbatim
68*> RMAX is REAL, dimension (2)
69*> RMAX(1) = Value of the largest test ratio of STRSYL
70*> RMAX(2) = Value of the largest test ratio of STRSYL3
71*> \endverbatim
72*>
73*> \param[out] NINFO
74*> \verbatim
75*> NINFO is INTEGER array, dimension (2)
76*> NINFO(1) = No. of times STRSYL returns an expected INFO
77*> NINFO(2) = No. of times STRSYL3 returns an expected INFO
78*> \endverbatim
79*>
80*> \param[out] KNT
81*> \verbatim
82*> KNT is INTEGER
83*> Total number of examples tested.
84*> \endverbatim
85
86*
87* -- LAPACK test routine --
88 SUBROUTINE ssyl01( THRESH, NFAIL, RMAX, NINFO, KNT )
89 IMPLICIT NONE
90*
91* -- LAPACK test routine --
92* -- LAPACK is a software package provided by Univ. of Tennessee, --
93* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
94*
95* .. Scalar Arguments ..
96 INTEGER KNT
97 REAL THRESH
98* ..
99* .. Array Arguments ..
100 INTEGER NFAIL( 3 ), NINFO( 2 )
101 REAL RMAX( 2 )
102* ..
103*
104* =====================================================================
105* ..
106* .. Parameters ..
107 REAL ZERO, ONE
108 parameter( zero = 0.0e+0, one = 1.0e+0 )
109 INTEGER MAXM, MAXN, LDSWORK
110 parameter( maxm = 101, maxn = 138, ldswork = 18 )
111* ..
112* .. Local Scalars ..
113 CHARACTER TRANA, TRANB
114 INTEGER I, INFO, IINFO, ISGN, ITRANA, ITRANB, J, KLA,
115 $ KUA, KLB, KUB, LIWORK, M, N
116 REAL ANRM, BNRM, BIGNUM, EPS, RES, RES1, RMUL,
117 $ SCALE, SCALE3, SMLNUM, TNRM, XNRM
118* ..
119* .. Local Arrays ..
120 REAL A( MAXM, MAXM ), B( MAXN, MAXN ),
121 $ C( MAXM, MAXN ), CC( MAXM, MAXN ),
122 $ X( MAXM, MAXN ),
123 $ DUML( MAXM ), DUMR( MAXN ),
124 $ D( MAX( MAXM, MAXN ) ), DUM( MAXN ),
125 $ SWORK( LDSWORK, 54 ), VM( 2 )
126 INTEGER ISEED( 4 ), IWORK( MAXM + MAXN + 2 ), IDUM( 2 )
127* ..
128* .. External Functions ..
129 LOGICAL SISNAN
130 REAL SLAMCH, SLANGE
131 EXTERNAL sisnan, slamch, slange
132* ..
133* .. External Subroutines ..
134 EXTERNAL slatmr, slacpy, sgemm, strsyl, strsyl3
135* ..
136* .. Intrinsic Functions ..
137 INTRINSIC abs, real, max
138* ..
139* .. Executable Statements ..
140*
141* Get machine parameters
142*
143 eps = slamch( 'P' )
144 smlnum = slamch( 'S' ) / eps
145 bignum = one / smlnum
146*
147 vm( 1 ) = one
148 vm( 2 ) = 0.05e+0
149*
150* Begin test loop
151*
152 ninfo( 1 ) = 0
153 ninfo( 2 ) = 0
154 nfail( 1 ) = 0
155 nfail( 2 ) = 0
156 nfail( 3 ) = 0
157 rmax( 1 ) = zero
158 rmax( 2 ) = zero
159 knt = 0
160 DO i = 1, 4
161 iseed( i ) = 1
162 END DO
163 scale = one
164 scale3 = one
165 liwork = maxm + maxn + 2
166 DO j = 1, 2
167 DO isgn = -1, 1, 2
168* Reset seed (overwritten by LATMR)
169 DO i = 1, 4
170 iseed( i ) = 1
171 END DO
172 DO m = 32, maxm, 71
173 kla = 0
174 kua = m - 1
175 CALL slatmr( m, m, 'S', iseed, 'N', d,
176 $ 6, one, one, 'T', 'N',
177 $ duml, 1, one, dumr, 1, one,
178 $ 'N', iwork, kla, kua, zero,
179 $ one, 'NO', a, maxm, iwork, iinfo )
180 DO i = 1, m
181 a( i, i ) = a( i, i ) * vm( j )
182 END DO
183 anrm = slange( 'M', m, m, a, maxm, dum )
184 DO n = 51, maxn, 47
185 klb = 0
186 kub = n - 1
187 CALL slatmr( n, n, 'S', iseed, 'N', d,
188 $ 6, one, one, 'T', 'N',
189 $ duml, 1, one, dumr, 1, one,
190 $ 'N', iwork, klb, kub, zero,
191 $ one, 'NO', b, maxn, iwork, iinfo )
192 bnrm = slange( 'M', n, n, b, maxn, dum )
193 tnrm = max( anrm, bnrm )
194 CALL slatmr( m, n, 'S', iseed, 'N', d,
195 $ 6, one, one, 'T', 'N',
196 $ duml, 1, one, dumr, 1, one,
197 $ 'N', iwork, m, n, zero, one,
198 $ 'NO', c, maxm, iwork, iinfo )
199 DO itrana = 1, 2
200 IF( itrana.EQ.1 ) THEN
201 trana = 'N'
202 END IF
203 IF( itrana.EQ.2 ) THEN
204 trana = 'T'
205 END IF
206 DO itranb = 1, 2
207 IF( itranb.EQ.1 ) THEN
208 tranb = 'N'
209 END IF
210 IF( itranb.EQ.2 ) THEN
211 tranb = 'T'
212 END IF
213 knt = knt + 1
214*
215 CALL slacpy( 'All', m, n, c, maxm, x, maxm)
216 CALL slacpy( 'All', m, n, c, maxm, cc, maxm)
217 CALL strsyl( trana, tranb, isgn, m, n,
218 $ a, maxm, b, maxn, x, maxm,
219 $ scale, iinfo )
220 IF( iinfo.NE.0 )
221 $ ninfo( 1 ) = ninfo( 1 ) + 1
222 xnrm = slange( 'M', m, n, x, maxm, dum )
223 rmul = one
224 IF( xnrm.GT.one .AND. tnrm.GT.one ) THEN
225 IF( xnrm.GT.bignum / tnrm ) THEN
226 rmul = one / max( xnrm, tnrm )
227 END IF
228 END IF
229 CALL sgemm( trana, 'N', m, n, m, rmul,
230 $ a, maxm, x, maxm, -scale*rmul,
231 $ c, maxm )
232 CALL sgemm( 'N', tranb, m, n, n,
233 $ real( isgn )*rmul, x, maxm, b,
234 $ maxn, one, c, maxm )
235 res1 = slange( 'M', m, n, c, maxm, dum )
236 res = res1 / max( smlnum, smlnum*xnrm,
237 $ ( ( rmul*tnrm )*eps )*xnrm )
238 IF( res.GT.thresh )
239 $ nfail( 1 ) = nfail( 1 ) + 1
240 IF( res.GT.rmax( 1 ) )
241 $ rmax( 1 ) = res
242*
243 CALL slacpy( 'All', m, n, c, maxm, x, maxm )
244 CALL slacpy( 'All', m, n, c, maxm, cc, maxm )
245 CALL strsyl3( trana, tranb, isgn, m, n,
246 $ a, maxm, b, maxn, x, maxm,
247 $ scale3, iwork, liwork,
248 $ swork, ldswork, info)
249 IF( info.NE.0 )
250 $ ninfo( 2 ) = ninfo( 2 ) + 1
251 xnrm = slange( 'M', m, n, x, maxm, dum )
252 rmul = one
253 IF( xnrm.GT.one .AND. tnrm.GT.one ) THEN
254 IF( xnrm.GT.bignum / tnrm ) THEN
255 rmul = one / max( xnrm, tnrm )
256 END IF
257 END IF
258 CALL sgemm( trana, 'N', m, n, m, rmul,
259 $ a, maxm, x, maxm, -scale3*rmul,
260 $ cc, maxm )
261 CALL sgemm( 'N', tranb, m, n, n,
262 $ real( isgn )*rmul, x, maxm, b,
263 $ maxn, one, cc, maxm )
264 res1 = slange( 'M', m, n, cc, maxm, dum )
265 res = res1 / max( smlnum, smlnum*xnrm,
266 $ ( ( rmul*tnrm )*eps )*xnrm )
267* Verify that TRSYL3 only flushes if TRSYL flushes (but
268* there may be cases where TRSYL3 avoid flushing).
269 IF( scale3.EQ.zero .AND. scale.GT.zero .OR.
270 $ iinfo.NE.info ) THEN
271 nfail( 3 ) = nfail( 3 ) + 1
272 END IF
273 IF( res.GT.thresh .OR. sisnan( res ) )
274 $ nfail( 2 ) = nfail( 2 ) + 1
275 IF( res.GT.rmax( 2 ) )
276 $ rmax( 2 ) = res
277 END DO
278 END DO
279 END DO
280 END DO
281 END DO
282 END DO
283*
284 RETURN
285*
286* End of SSYL01
287*
288 END
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine slatmr(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, RSIGN, GRADE, DL, MODEL, CONDL, DR, MODER, CONDR, PIVTNG, IPIVOT, KL, KU, SPARSE, ANORM, PACK, A, LDA, IWORK, INFO)
SLATMR
Definition: slatmr.f:471
subroutine strsyl(TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, LDC, SCALE, INFO)
STRSYL
Definition: strsyl.f:164
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
subroutine ssyl01(THRESH, NFAIL, RMAX, NINFO, KNT)
SSYL01
Definition: ssyl01.f:89
subroutine strsyl3(TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, LDC, SCALE, IWORK, LIWORK, SWORK, LDSWORK, INFO)
STRSYL3
Definition: strsyl3.f:181