LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
csyl01.f
Go to the documentation of this file.
1*> \brief \b CSYL01
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 CSYL01( 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*> CSYL01 tests CTRSYL and CTRSYL3, routines for solving the Sylvester matrix
29*> equation
30*>
31*> op(A)*X + ISGN*X*op(B) = scale*C,
32*>
33*> where op(A) and op(B) are both upper triangular form, op() represents an
34*> optional conjugate 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 CGET35 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 CTRSYL exceeds threshold THRESH
62*> NFAIL(2) = No. of times residual CTRSYL3 exceeds threshold THRESH
63*> NFAIL(3) = No. of times CTRSYL3 and CTRSYL deviate
64*> \endverbatim
65*>
66*> \param[out] RMAX
67*> \verbatim
68*> RMAX is DOUBLE PRECISION array, dimension (2)
69*> RMAX(1) = Value of the largest test ratio of CTRSYL
70*> RMAX(2) = Value of the largest test ratio of CTRSYL3
71*> \endverbatim
72*>
73*> \param[out] NINFO
74*> \verbatim
75*> NINFO is INTEGER array, dimension (2)
76*> NINFO(1) = No. of times CTRSYL where INFO is nonzero
77*> NINFO(2) = No. of times CTRSYL3 where INFO is nonzero
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 csyl01( 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 COMPLEX CONE
108 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
109 REAL ONE, ZERO
110 parameter( zero = 0.0e+0, one = 1.0e+0 )
111 INTEGER MAXM, MAXN, LDSWORK
112 parameter( maxm = 101, maxn = 138, ldswork = 18 )
113* ..
114* .. Local Scalars ..
115 CHARACTER TRANA, TRANB
116 INTEGER I, INFO, IINFO, ISGN, ITRANA, ITRANB, J, KLA,
117 $ KUA, KLB, KUB, M, N
118 REAL ANRM, BNRM, BIGNUM, EPS, RES, RES1,
119 $ SCALE, SCALE3, SMLNUM, TNRM, XNRM
120 COMPLEX RMUL
121* ..
122* .. Local Arrays ..
123 COMPLEX A( MAXM, MAXM ), B( MAXN, MAXN ),
124 $ C( MAXM, MAXN ), CC( MAXM, MAXN ),
125 $ X( MAXM, MAXN ),
126 $ DUML( MAXM ), DUMR( MAXN ),
127 $ D( MIN( MAXM, MAXN ) )
128 REAL SWORK( LDSWORK, 54 ), DUM( MAXN ), VM( 2 )
129 INTEGER ISEED( 4 ), IWORK( MAXM + MAXN + 2 )
130* ..
131* .. External Functions ..
132 LOGICAL SISNAN
133 REAL SLAMCH, CLANGE
134 EXTERNAL sisnan, slamch, clange
135* ..
136* .. External Subroutines ..
137 EXTERNAL clatmr, clacpy, cgemm, ctrsyl, ctrsyl3
138* ..
139* .. Intrinsic Functions ..
140 INTRINSIC abs, real, max
141* ..
142* .. Executable Statements ..
143*
144* Get machine parameters
145*
146 eps = slamch( 'P' )
147 smlnum = slamch( 'S' ) / eps
148 bignum = one / smlnum
149*
150* Expect INFO = 0
151 vm( 1 ) = one
152* Expect INFO = 1
153 vm( 2 ) = 0.5e+0
154*
155* Begin test loop
156*
157 ninfo( 1 ) = 0
158 ninfo( 2 ) = 0
159 nfail( 1 ) = 0
160 nfail( 2 ) = 0
161 nfail( 3 ) = 0
162 rmax( 1 ) = zero
163 rmax( 2 ) = zero
164 knt = 0
165 iseed( 1 ) = 1
166 iseed( 2 ) = 1
167 iseed( 3 ) = 1
168 iseed( 4 ) = 1
169 scale = one
170 scale3 = one
171 DO j = 1, 2
172 DO isgn = -1, 1, 2
173* Reset seed (overwritten by LATMR)
174 iseed( 1 ) = 1
175 iseed( 2 ) = 1
176 iseed( 3 ) = 1
177 iseed( 4 ) = 1
178 DO m = 32, maxm, 23
179 kla = 0
180 kua = m - 1
181 CALL clatmr( m, m, 'S', iseed, 'N', d,
182 $ 6, one, cone, 'T', 'N',
183 $ duml, 1, one, dumr, 1, one,
184 $ 'N', iwork, kla, kua, zero,
185 $ one, 'NO', a, maxm, iwork,
186 $ iinfo )
187 DO i = 1, m
188 a( i, i ) = a( i, i ) * vm( j )
189 END DO
190 anrm = clange( 'M', m, m, a, maxm, dum )
191 DO n = 51, maxn, 29
192 klb = 0
193 kub = n - 1
194 CALL clatmr( n, n, 'S', iseed, 'N', d,
195 $ 6, one, cone, 'T', 'N',
196 $ duml, 1, one, dumr, 1, one,
197 $ 'N', iwork, klb, kub, zero,
198 $ one, 'NO', b, maxn, iwork,
199 $ iinfo )
200 DO i = 1, n
201 b( i, i ) = b( i, i ) * vm( j )
202 END DO
203 bnrm = clange( 'M', n, n, b, maxn, dum )
204 tnrm = max( anrm, bnrm )
205 CALL clatmr( m, n, 'S', iseed, 'N', d,
206 $ 6, one, cone, 'T', 'N',
207 $ duml, 1, one, dumr, 1, one,
208 $ 'N', iwork, m, n, zero, one,
209 $ 'NO', c, maxm, iwork, iinfo )
210 DO itrana = 1, 2
211 IF( itrana.EQ.1 )
212 $ trana = 'N'
213 IF( itrana.EQ.2 )
214 $ trana = 'C'
215 DO itranb = 1, 2
216 IF( itranb.EQ.1 )
217 $ tranb = 'N'
218 IF( itranb.EQ.2 )
219 $ tranb = 'C'
220 knt = knt + 1
221*
222 CALL clacpy( 'All', m, n, c, maxm, x, maxm)
223 CALL clacpy( 'All', m, n, c, maxm, cc, maxm)
224 CALL ctrsyl( trana, tranb, isgn, m, n,
225 $ a, maxm, b, maxn, x, maxm,
226 $ scale, iinfo )
227 IF( iinfo.NE.0 )
228 $ ninfo( 1 ) = ninfo( 1 ) + 1
229 xnrm = clange( 'M', m, n, x, maxm, dum )
230 rmul = cone
231 IF( xnrm.GT.one .AND. tnrm.GT.one ) THEN
232 IF( xnrm.GT.bignum / tnrm ) THEN
233 rmul = cone / max( xnrm, tnrm )
234 END IF
235 END IF
236 CALL cgemm( trana, 'N', m, n, m, rmul,
237 $ a, maxm, x, maxm, -scale*rmul,
238 $ cc, maxm )
239 CALL cgemm( 'N', tranb, m, n, n,
240 $ real( isgn )*rmul, x, maxm, b,
241 $ maxn, cone, cc, maxm )
242 res1 = clange( 'M', m, n, cc, maxm, dum )
243 res = res1 / max( smlnum, smlnum*xnrm,
244 $ ( ( abs( rmul )*tnrm )*eps )*xnrm )
245 IF( res.GT.thresh )
246 $ nfail( 1 ) = nfail( 1 ) + 1
247 IF( res.GT.rmax( 1 ) )
248 $ rmax( 1 ) = res
249*
250 CALL clacpy( 'All', m, n, c, maxm, x, maxm )
251 CALL clacpy( 'All', m, n, c, maxm, cc, maxm )
252 CALL ctrsyl3( trana, tranb, isgn, m, n,
253 $ a, maxm, b, maxn, x, maxm,
254 $ scale3, swork, ldswork, info)
255 IF( info.NE.0 )
256 $ ninfo( 2 ) = ninfo( 2 ) + 1
257 xnrm = clange( 'M', m, n, x, maxm, dum )
258 rmul = cone
259 IF( xnrm.GT.one .AND. tnrm.GT.one ) THEN
260 IF( xnrm.GT.bignum / tnrm ) THEN
261 rmul = cone / max( xnrm, tnrm )
262 END IF
263 END IF
264 CALL cgemm( trana, 'N', m, n, m, rmul,
265 $ a, maxm, x, maxm, -scale3*rmul,
266 $ cc, maxm )
267 CALL cgemm( 'N', tranb, m, n, n,
268 $ real( isgn )*rmul, x, maxm, b,
269 $ maxn, cone, cc, maxm )
270 res1 = clange( 'M', m, n, cc, maxm, dum )
271 res = res1 / max( smlnum, smlnum*xnrm,
272 $ ( ( abs( rmul )*tnrm )*eps )*xnrm )
273* Verify that TRSYL3 only flushes if TRSYL flushes (but
274* there may be cases where TRSYL3 avoid flushing).
275 IF( scale3.EQ.zero .AND. scale.GT.zero .OR.
276 $ iinfo.NE.info ) THEN
277 nfail( 3 ) = nfail( 3 ) + 1
278 END IF
279 IF( res.GT.thresh .OR. sisnan( res ) )
280 $ nfail( 2 ) = nfail( 2 ) + 1
281 IF( res.GT.rmax( 2 ) )
282 $ rmax( 2 ) = res
283 END DO
284 END DO
285 END DO
286 END DO
287 END DO
288 END DO
289*
290 RETURN
291*
292* End of CSYL01
293*
294 END
subroutine csyl01(THRESH, NFAIL, RMAX, NINFO, KNT)
CSYL01
Definition: csyl01.f:89
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
subroutine clatmr(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)
CLATMR
Definition: clatmr.f:490
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
subroutine ctrsyl(TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, LDC, SCALE, INFO)
CTRSYL
Definition: ctrsyl.f:157
subroutine ctrsyl3(TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, LDC, SCALE, SWORK, LDSWORK, INFO)
CTRSYL3
Definition: ctrsyl3.f:156