LAPACK 3.12.1
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 DUML( MAXM ), DUMR( MAXN ),
124 $ D( MAX( MAXM, MAXN ) )
125 REAL DUM( MAXN ), VM( 2 )
126 INTEGER ISEED( 4 ), IWORK( MAXM + MAXN + 2 )
127* ..
128* .. Allocatable Arrays ..
129 INTEGER AllocateStatus
130 COMPLEX, DIMENSION(:,:), ALLOCATABLE :: A, B, C, CC, X
131 REAL, DIMENSION(:,:), ALLOCATABLE :: SWORK
132* ..
133* .. External Functions ..
134 LOGICAL SISNAN
135 REAL SLAMCH, CLANGE
136 EXTERNAL sisnan, slamch, clange
137* ..
138* .. External Subroutines ..
139 EXTERNAL clatmr, clacpy, cgemm, ctrsyl, ctrsyl3
140* ..
141* .. Intrinsic Functions ..
142 INTRINSIC abs, real, max
143* ..
144* .. Allocate memory dynamically ..
145 ALLOCATE ( a( maxm, maxm ), stat = allocatestatus )
146 IF( allocatestatus /= 0 ) stop "*** Not enough memory ***"
147 ALLOCATE ( b( maxn, maxn ), stat = allocatestatus )
148 IF( allocatestatus /= 0 ) stop "*** Not enough memory ***"
149 ALLOCATE ( c( maxm, maxn ), stat = allocatestatus )
150 IF( allocatestatus /= 0 ) stop "*** Not enough memory ***"
151 ALLOCATE ( cc( maxm, maxn ), stat = allocatestatus )
152 IF( allocatestatus /= 0 ) stop "*** Not enough memory ***"
153 ALLOCATE ( x( maxm, maxn ), stat = allocatestatus )
154 IF( allocatestatus /= 0 ) stop "*** Not enough memory ***"
155 ALLOCATE ( swork( ldswork, 54 ), stat = allocatestatus )
156 IF( allocatestatus /= 0 ) stop "*** Not enough memory ***"
157* ..
158* .. Executable Statements ..
159*
160* Get machine parameters
161*
162 eps = slamch( 'P' )
163 smlnum = slamch( 'S' ) / eps
164 bignum = one / smlnum
165*
166* Expect INFO = 0
167 vm( 1 ) = one
168* Expect INFO = 1
169 vm( 2 ) = 0.5e+0
170*
171* Begin test loop
172*
173 ninfo( 1 ) = 0
174 ninfo( 2 ) = 0
175 nfail( 1 ) = 0
176 nfail( 2 ) = 0
177 nfail( 3 ) = 0
178 rmax( 1 ) = zero
179 rmax( 2 ) = zero
180 knt = 0
181 iseed( 1 ) = 1
182 iseed( 2 ) = 1
183 iseed( 3 ) = 1
184 iseed( 4 ) = 1
185 scale = one
186 scale3 = one
187 DO j = 1, 2
188 DO isgn = -1, 1, 2
189* Reset seed (overwritten by LATMR)
190 iseed( 1 ) = 1
191 iseed( 2 ) = 1
192 iseed( 3 ) = 1
193 iseed( 4 ) = 1
194 DO m = 32, maxm, 23
195 kla = 0
196 kua = m - 1
197 CALL clatmr( m, m, 'S', iseed, 'N', d,
198 $ 6, one, cone, 'T', 'N',
199 $ duml, 1, one, dumr, 1, one,
200 $ 'N', iwork, kla, kua, zero,
201 $ one, 'NO', a, maxm, iwork,
202 $ iinfo )
203 DO i = 1, m
204 a( i, i ) = a( i, i ) * vm( j )
205 END DO
206 anrm = clange( 'M', m, m, a, maxm, dum )
207 DO n = 51, maxn, 29
208 klb = 0
209 kub = n - 1
210 CALL clatmr( n, n, 'S', iseed, 'N', d,
211 $ 6, one, cone, 'T', 'N',
212 $ duml, 1, one, dumr, 1, one,
213 $ 'N', iwork, klb, kub, zero,
214 $ one, 'NO', b, maxn, iwork,
215 $ iinfo )
216 DO i = 1, n
217 b( i, i ) = b( i, i ) * vm( j )
218 END DO
219 bnrm = clange( 'M', n, n, b, maxn, dum )
220 tnrm = max( anrm, bnrm )
221 CALL clatmr( m, n, 'S', iseed, 'N', d,
222 $ 6, one, cone, 'T', 'N',
223 $ duml, 1, one, dumr, 1, one,
224 $ 'N', iwork, m, n, zero, one,
225 $ 'NO', c, maxm, iwork, iinfo )
226 DO itrana = 1, 2
227 IF( itrana.EQ.1 )
228 $ trana = 'N'
229 IF( itrana.EQ.2 )
230 $ trana = 'C'
231 DO itranb = 1, 2
232 IF( itranb.EQ.1 )
233 $ tranb = 'N'
234 IF( itranb.EQ.2 )
235 $ tranb = 'C'
236 knt = knt + 1
237*
238 CALL clacpy( 'All', m, n, c, maxm, x, maxm)
239 CALL clacpy( 'All', m, n, c, maxm, cc, maxm)
240 CALL ctrsyl( trana, tranb, isgn, m, n,
241 $ a, maxm, b, maxn, x, maxm,
242 $ scale, iinfo )
243 IF( iinfo.NE.0 )
244 $ ninfo( 1 ) = ninfo( 1 ) + 1
245 xnrm = clange( 'M', m, n, x, maxm, dum )
246 rmul = cone
247 IF( xnrm.GT.one .AND. tnrm.GT.one ) THEN
248 IF( xnrm.GT.bignum / tnrm ) THEN
249 rmul = cone / max( xnrm, tnrm )
250 END IF
251 END IF
252 CALL cgemm( trana, 'N', m, n, m, rmul,
253 $ a, maxm, x, maxm, -scale*rmul,
254 $ cc, maxm )
255 CALL cgemm( 'N', tranb, m, n, n,
256 $ real( isgn )*rmul, x, maxm, b,
257 $ maxn, cone, cc, maxm )
258 res1 = clange( 'M', m, n, cc, maxm, dum )
259 res = res1 / max( smlnum, smlnum*xnrm,
260 $ ( ( abs( rmul )*tnrm )*eps )*xnrm )
261 IF( res.GT.thresh )
262 $ nfail( 1 ) = nfail( 1 ) + 1
263 IF( res.GT.rmax( 1 ) )
264 $ rmax( 1 ) = res
265*
266 CALL clacpy( 'All', m, n, c, maxm, x, maxm )
267 CALL clacpy( 'All', m, n, c, maxm, cc, maxm )
268 CALL ctrsyl3( trana, tranb, isgn, m, n,
269 $ a, maxm, b, maxn, x, maxm,
270 $ scale3, swork, ldswork, info)
271 IF( info.NE.0 )
272 $ ninfo( 2 ) = ninfo( 2 ) + 1
273 xnrm = clange( 'M', m, n, x, maxm, dum )
274 rmul = cone
275 IF( xnrm.GT.one .AND. tnrm.GT.one ) THEN
276 IF( xnrm.GT.bignum / tnrm ) THEN
277 rmul = cone / max( xnrm, tnrm )
278 END IF
279 END IF
280 CALL cgemm( trana, 'N', m, n, m, rmul,
281 $ a, maxm, x, maxm, -scale3*rmul,
282 $ cc, maxm )
283 CALL cgemm( 'N', tranb, m, n, n,
284 $ real( isgn )*rmul, x, maxm, b,
285 $ maxn, cone, cc, maxm )
286 res1 = clange( 'M', m, n, cc, maxm, dum )
287 res = res1 / max( smlnum, smlnum*xnrm,
288 $ ( ( abs( rmul )*tnrm )*eps )*xnrm )
289* Verify that TRSYL3 only flushes if TRSYL flushes (but
290* there may be cases where TRSYL3 avoid flushing).
291 IF( scale3.EQ.zero .AND. scale.GT.zero .OR.
292 $ iinfo.NE.info ) THEN
293 nfail( 3 ) = nfail( 3 ) + 1
294 END IF
295 IF( res.GT.thresh .OR. sisnan( res ) )
296 $ nfail( 2 ) = nfail( 2 ) + 1
297 IF( res.GT.rmax( 2 ) )
298 $ rmax( 2 ) = res
299 END DO
300 END DO
301 END DO
302 END DO
303 END DO
304 END DO
305*
306 DEALLOCATE (a, stat = allocatestatus)
307 DEALLOCATE (b, stat = allocatestatus)
308 DEALLOCATE (c, stat = allocatestatus)
309 DEALLOCATE (cc, stat = allocatestatus)
310 DEALLOCATE (x, stat = allocatestatus)
311 DEALLOCATE (swork, stat = allocatestatus)
312*
313 RETURN
314*
315* End of CSYL01
316*
317 END
clatmr
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
csyl01
subroutine csyl01(thresh, nfail, rmax, ninfo, knt)
CSYL01
Definition csyl01.f:89
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:101
ctrsyl3
subroutine ctrsyl3(trana, tranb, isgn, m, n, a, lda, b, ldb, c, ldc, scale, swork, ldswork, info)
CTRSYL3
Definition ctrsyl3.f:169
ctrsyl
subroutine ctrsyl(trana, tranb, isgn, m, n, a, lda, b, ldb, c, ldc, scale, info)
CTRSYL
Definition ctrsyl.f:155