LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dsyl01.f
Go to the documentation of this file.
1*> \brief \b DSYL01
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 DSYL01( THRESH, NFAIL, RMAX, NINFO, KNT )
12*
13* .. Scalar Arguments ..
14* INTEGER KNT
15* DOUBLE PRECISION THRESH
16* ..
17* .. Array Arguments ..
18* INTEGER NFAIL( 3 ), NINFO( 2 )
19* DOUBLE PRECISION RMAX( 2 )
20* ..
21*
22*
23*> \par Purpose:
24* =============
25*>
26*> \verbatim
27*>
28*> DSYL01 tests DTRSYL and DTRSYL3, 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 DGET35 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 DOUBLE PRECISION
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 DTRSYL exceeds threshold THRESH
62*> NFAIL(2) = No. of times residual DTRSYL3 exceeds threshold THRESH
63*> NFAIL(3) = No. of times DTRSYL3 and DTRSYL deviate
64*> \endverbatim
65*>
66*> \param[out] RMAX
67*> \verbatim
68*> RMAX is DOUBLE PRECISION, dimension (2)
69*> RMAX(1) = Value of the largest test ratio of DTRSYL
70*> RMAX(2) = Value of the largest test ratio of DTRSYL3
71*> \endverbatim
72*>
73*> \param[out] NINFO
74*> \verbatim
75*> NINFO is INTEGER array, dimension (2)
76*> NINFO(1) = No. of times DTRSYL returns an expected INFO
77*> NINFO(2) = No. of times DTRSYL3 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 dsyl01( 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 DOUBLE PRECISION THRESH
98* ..
99* .. Array Arguments ..
100 INTEGER NFAIL( 3 ), NINFO( 2 )
101 DOUBLE PRECISION RMAX( 2 )
102* ..
103*
104* =====================================================================
105* ..
106* .. Parameters ..
107 DOUBLE PRECISION ZERO, ONE
108 parameter( zero = 0.0d0, one = 1.0d0 )
109 INTEGER MAXM, MAXN, LDSWORK
110 parameter( maxm = 245, maxn = 192, ldswork = 36 )
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 DOUBLE PRECISION ANRM, BNRM, BIGNUM, EPS, RES, RES1, RMUL,
117 $ SCALE, SCALE3, SMLNUM, TNRM, XNRM
118* ..
119* .. Local Arrays ..
120 DOUBLE PRECISION DUML( MAXM ), DUMR( MAXN ),
121 $ D( MAX( MAXM, MAXN ) ), DUM( MAXN ),
122 $ VM( 2 )
123 INTEGER ISEED( 4 ), IWORK( MAXM + MAXN + 2 )
124* ..
125* .. Allocatable Arrays ..
126 INTEGER AllocateStatus
127 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: A, B, C, CC, X,
128 $ SWORK
129* ..
130* .. External Functions ..
131 LOGICAL DISNAN
132 DOUBLE PRECISION DLAMCH, DLANGE
133 EXTERNAL dlamch, dlange
134* ..
135* .. External Subroutines ..
136 EXTERNAL dlatmr, dlacpy, dgemm, dtrsyl, dtrsyl3
137* ..
138* .. Intrinsic Functions ..
139 INTRINSIC abs, dble, max
140* ..
141* .. Allocate memory dynamically ..
142 ALLOCATE ( a( maxm, maxm ), stat = allocatestatus )
143 IF( allocatestatus /= 0 ) stop "*** Not enough memory ***"
144 ALLOCATE ( b( maxn, maxn ), stat = allocatestatus )
145 IF( allocatestatus /= 0 ) stop "*** Not enough memory ***"
146 ALLOCATE ( c( maxm, maxn ), stat = allocatestatus )
147 IF( allocatestatus /= 0 ) stop "*** Not enough memory ***"
148 ALLOCATE ( cc( maxm, maxn ), stat = allocatestatus )
149 IF( allocatestatus /= 0 ) stop "*** Not enough memory ***"
150 ALLOCATE ( x( maxm, maxn ), stat = allocatestatus )
151 IF( allocatestatus /= 0 ) stop "*** Not enough memory ***"
152 ALLOCATE ( swork( ldswork, 126 ), stat = allocatestatus )
153 IF( allocatestatus /= 0 ) stop "*** Not enough memory ***"
154* ..
155* .. Executable Statements ..
156*
157* Get machine parameters
158*
159 eps = dlamch( 'P' )
160 smlnum = dlamch( 'S' ) / eps
161 bignum = one / smlnum
162*
163 vm( 1 ) = one
164 vm( 2 ) = 0.000001d+0
165*
166* Begin test loop
167*
168 ninfo( 1 ) = 0
169 ninfo( 2 ) = 0
170 nfail( 1 ) = 0
171 nfail( 2 ) = 0
172 nfail( 3 ) = 0
173 rmax( 1 ) = zero
174 rmax( 2 ) = zero
175 knt = 0
176 DO i = 1, 4
177 iseed( i ) = 1
178 END DO
179 scale = one
180 scale3 = one
181 liwork = maxm + maxn + 2
182 DO j = 1, 2
183 DO isgn = -1, 1, 2
184* Reset seed (overwritten by LATMR)
185 DO i = 1, 4
186 iseed( i ) = 1
187 END DO
188 DO m = 32, maxm, 71
189 kla = 0
190 kua = m - 1
191 CALL dlatmr( m, m, 'S', iseed, 'N', d,
192 $ 6, one, one, 'T', 'N',
193 $ duml, 1, one, dumr, 1, one,
194 $ 'N', iwork, kla, kua, zero,
195 $ one, 'NO', a, maxm, iwork, iinfo )
196 DO i = 1, m
197 a( i, i ) = a( i, i ) * vm( j )
198 END DO
199 anrm = dlange( 'M', m, m, a, maxm, dum )
200 DO n = 51, maxn, 47
201 klb = 0
202 kub = n - 1
203 CALL dlatmr( n, n, 'S', iseed, 'N', d,
204 $ 6, one, one, 'T', 'N',
205 $ duml, 1, one, dumr, 1, one,
206 $ 'N', iwork, klb, kub, zero,
207 $ one, 'NO', b, maxn, iwork, iinfo )
208 bnrm = dlange( 'M', n, n, b, maxn, dum )
209 tnrm = max( anrm, bnrm )
210 CALL dlatmr( m, n, 'S', iseed, 'N', d,
211 $ 6, one, one, 'T', 'N',
212 $ duml, 1, one, dumr, 1, one,
213 $ 'N', iwork, m, n, zero, one,
214 $ 'NO', c, maxm, iwork, iinfo )
215 DO itrana = 1, 2
216 IF( itrana.EQ.1 ) THEN
217 trana = 'N'
218 END IF
219 IF( itrana.EQ.2 ) THEN
220 trana = 'T'
221 END IF
222 DO itranb = 1, 2
223 IF( itranb.EQ.1 ) THEN
224 tranb = 'N'
225 END IF
226 IF( itranb.EQ.2 ) THEN
227 tranb = 'T'
228 END IF
229 knt = knt + 1
230*
231 CALL dlacpy( 'All', m, n, c, maxm, x, maxm)
232 CALL dlacpy( 'All', m, n, c, maxm, cc, maxm)
233 CALL dtrsyl( trana, tranb, isgn, m, n,
234 $ a, maxm, b, maxn, x, maxm,
235 $ scale, iinfo )
236 IF( iinfo.NE.0 )
237 $ ninfo( 1 ) = ninfo( 1 ) + 1
238 xnrm = dlange( 'M', m, n, x, maxm, dum )
239 rmul = one
240 IF( xnrm.GT.one .AND. tnrm.GT.one ) THEN
241 IF( xnrm.GT.bignum / tnrm ) THEN
242 rmul = one / max( xnrm, tnrm )
243 END IF
244 END IF
245 CALL dgemm( trana, 'N', m, n, m, rmul,
246 $ a, maxm, x, maxm, -scale*rmul,
247 $ cc, maxm )
248 CALL dgemm( 'N', tranb, m, n, n,
249 $ dble( isgn )*rmul, x, maxm, b,
250 $ maxn, one, cc, maxm )
251 res1 = dlange( 'M', m, n, cc, maxm, dum )
252 res = res1 / max( smlnum, smlnum*xnrm,
253 $ ( ( rmul*tnrm )*eps )*xnrm )
254 IF( res.GT.thresh )
255 $ nfail( 1 ) = nfail( 1 ) + 1
256 IF( res.GT.rmax( 1 ) )
257 $ rmax( 1 ) = res
258*
259 CALL dlacpy( 'All', m, n, c, maxm, x, maxm )
260 CALL dlacpy( 'All', m, n, c, maxm, cc, maxm )
261 CALL dtrsyl3( trana, tranb, isgn, m, n,
262 $ a, maxm, b, maxn, x, maxm,
263 $ scale3, iwork, liwork,
264 $ swork, ldswork, info)
265 IF( info.NE.0 )
266 $ ninfo( 2 ) = ninfo( 2 ) + 1
267 xnrm = dlange( 'M', m, n, x, maxm, dum )
268 rmul = one
269 IF( xnrm.GT.one .AND. tnrm.GT.one ) THEN
270 IF( xnrm.GT.bignum / tnrm ) THEN
271 rmul = one / max( xnrm, tnrm )
272 END IF
273 END IF
274 CALL dgemm( trana, 'N', m, n, m, rmul,
275 $ a, maxm, x, maxm, -scale3*rmul,
276 $ cc, maxm )
277 CALL dgemm( 'N', tranb, m, n, n,
278 $ dble( isgn )*rmul, x, maxm, b,
279 $ maxn, one, cc, maxm )
280 res1 = dlange( 'M', m, n, cc, maxm, dum )
281 res = res1 / max( smlnum, smlnum*xnrm,
282 $ ( ( rmul*tnrm )*eps )*xnrm )
283* Verify that TRSYL3 only flushes if TRSYL flushes (but
284* there may be cases where TRSYL3 avoid flushing).
285 IF( scale3.EQ.zero .AND. scale.GT.zero .OR.
286 $ iinfo.NE.info ) THEN
287 nfail( 3 ) = nfail( 3 ) + 1
288 END IF
289 IF( res.GT.thresh .OR. disnan( res ) )
290 $ nfail( 2 ) = nfail( 2 ) + 1
291 IF( res.GT.rmax( 2 ) )
292 $ rmax( 2 ) = res
293 END DO
294 END DO
295 END DO
296 END DO
297 END DO
298 END DO
299*
300 DEALLOCATE (a, stat = allocatestatus)
301 DEALLOCATE (b, stat = allocatestatus)
302 DEALLOCATE (c, stat = allocatestatus)
303 DEALLOCATE (cc, stat = allocatestatus)
304 DEALLOCATE (x, stat = allocatestatus)
305 DEALLOCATE (swork, stat = allocatestatus)
306*
307 RETURN
308*
309* End of DSYL01
310*
311 END
dlatmr
subroutine dlatmr(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)
DLATMR
Definition dlatmr.f:471
dsyl01
subroutine dsyl01(thresh, nfail, rmax, ninfo, knt)
DSYL01
Definition dsyl01.f:89
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:101
dtrsyl3
subroutine dtrsyl3(trana, tranb, isgn, m, n, a, lda, b, ldb, c, ldc, scale, iwork, liwork, swork, ldswork, info)
DTRSYL3
Definition dtrsyl3.f:200
dtrsyl
subroutine dtrsyl(trana, tranb, isgn, m, n, a, lda, b, ldb, c, ldc, scale, info)
DTRSYL
Definition dtrsyl.f:162