LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zsyl01.f
Go to the documentation of this file.
1*> \brief \b ZSYL01
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 ZSYL01( 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*> ZSYL01 tests ZTRSYL and ZTRSYL3, 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 ZGET35 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 ZTRSYL exceeds threshold THRESH
62*> NFAIL(2) = No. of times residual ZTRSYL3 exceeds threshold THRESH
63*> NFAIL(3) = No. of times ZTRSYL3 and ZTRSYL 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 ZTRSYL
70*> RMAX(2) = Value of the largest test ratio of ZTRSYL3
71*> \endverbatim
72*>
73*> \param[out] NINFO
74*> \verbatim
75*> NINFO is INTEGER array, dimension (2)
76*> NINFO(1) = No. of times ZTRSYL returns an expected INFO
77*> NINFO(2) = No. of times ZTRSYL3 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 zsyl01( 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 COMPLEX*16 CONE
108 parameter( cone = ( 1.0d0, 0.0d+0 ) )
109 DOUBLE PRECISION ONE, ZERO
110 parameter( zero = 0.0d+0, one = 1.0d+0 )
111 INTEGER MAXM, MAXN, LDSWORK
112 parameter( maxm = 185, maxn = 192, ldswork = 36 )
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 DOUBLE PRECISION ANRM, BNRM, BIGNUM, EPS, RES, RES1,
119 $ SCALE, SCALE3, SMLNUM, TNRM, XNRM
120 COMPLEX*16 RMUL
121* ..
122* .. Local Arrays ..
123 COMPLEX*16 DUML( MAXM ), DUMR( MAXN ),
124 $ D( MAX( MAXM, MAXN ) )
125 DOUBLE PRECISION DUM( MAXN ), VM( 2 )
126 INTEGER ISEED( 4 ), IWORK( MAXM + MAXN + 2 )
127* ..
128* .. Allocatable Arrays ..
129 INTEGER AllocateStatus
130 COMPLEX*16, DIMENSION(:,:), ALLOCATABLE :: A, B, C, CC, X
131 DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: SWORK
132* ..
133* .. External Functions ..
134 LOGICAL DISNAN
135 DOUBLE PRECISION DLAMCH, ZLANGE
136 EXTERNAL disnan, dlamch, zlange
137* ..
138* .. External Subroutines ..
139 EXTERNAL zlatmr, zlacpy, zgemm, ztrsyl, ztrsyl3
140* ..
141* .. Intrinsic Functions ..
142 INTRINSIC abs, dble, max, sqrt
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, 103 ), stat = allocatestatus )
156 IF( allocatestatus /= 0 ) stop "*** Not enough memory ***"
157* ..
158* .. Executable Statements ..
159*
160* Get machine parameters
161*
162 eps = dlamch( 'P' )
163 smlnum = dlamch( 'S' ) / eps
164 bignum = one / smlnum
165*
166* Expect INFO = 0
167 vm( 1 ) = one
168* Expect INFO = 1
169 vm( 2 ) = 0.05d+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, 51
195 kla = 0
196 kua = m - 1
197 CALL zlatmr( 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 = zlange( 'M', m, m, a, maxm, dum )
207 DO n = 51, maxn, 47
208 klb = 0
209 kub = n - 1
210 CALL zlatmr( 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 = zlange( 'M', n, n, b, maxn, dum )
220 tnrm = max( anrm, bnrm )
221 CALL zlatmr( 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 zlacpy( 'All', m, n, c, maxm, x, maxm)
239 CALL zlacpy( 'All', m, n, c, maxm, cc, maxm)
240 CALL ztrsyl( 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 = zlange( '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 zgemm( trana, 'N', m, n, m, rmul,
253 $ a, maxm, x, maxm, -scale*rmul,
254 $ cc, maxm )
255 CALL zgemm( 'N', tranb, m, n, n,
256 $ dble( isgn )*rmul, x, maxm, b,
257 $ maxn, cone, cc, maxm )
258 res1 = zlange( '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 zlacpy( 'All', m, n, c, maxm, x, maxm )
267 CALL zlacpy( 'All', m, n, c, maxm, cc, maxm )
268 CALL ztrsyl3( 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 = zlange( '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 zgemm( trana, 'N', m, n, m, rmul,
281 $ a, maxm, x, maxm, -scale3*rmul,
282 $ cc, maxm )
283 CALL zgemm( 'N', tranb, m, n, n,
284 $ dble( isgn )*rmul, x, maxm, b,
285 $ maxn, cone, cc, maxm )
286 res1 = zlange( '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. disnan( 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 ZSYL01
316*
317 END
zgemm
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:101
ztrsyl3
subroutine ztrsyl3(trana, tranb, isgn, m, n, a, lda, b, ldb, c, ldc, scale, swork, ldswork, info)
ZTRSYL3
Definition ztrsyl3.f:170
ztrsyl
subroutine ztrsyl(trana, tranb, isgn, m, n, a, lda, b, ldb, c, ldc, scale, info)
ZTRSYL
Definition ztrsyl.f:155
zlatmr
subroutine zlatmr(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)
ZLATMR
Definition zlatmr.f:490
zsyl01
subroutine zsyl01(thresh, nfail, rmax, ninfo, knt)
ZSYL01
Definition zsyl01.f:89