LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ csyl01()

subroutine csyl01 ( real  thresh,
integer, dimension( 3 )  nfail,
real, dimension( 2 )  rmax,
integer, dimension( 2 )  ninfo,
integer  knt 
)

CSYL01

Purpose:
 CSYL01 tests CTRSYL and CTRSYL3, routines for solving the Sylvester matrix
 equation

    op(A)*X + ISGN*X*op(B) = scale*C,

 where op(A) and op(B) are both upper triangular form, op() represents an
 optional conjugate transpose, and ISGN can be -1 or +1. Scale is an output
 less than or equal to 1, chosen to avoid overflow in X.

 The test code verifies that the following residual does not exceed
 the provided threshold:

    norm(op(A)*X + ISGN*X*op(B) - scale*C) /
        (EPS*max(norm(A),norm(B))*norm(X))

 This routine complements CGET35 by testing with larger,
 random matrices, of which some require rescaling of X to avoid overflow.
Parameters
[in]THRESH
          THRESH is REAL
          A test will count as "failed" if the residual, computed as
          described above, exceeds THRESH.
[out]NFAIL
          NFAIL is INTEGER array, dimension (3)
          NFAIL(1) = No. of times residual CTRSYL exceeds threshold THRESH
          NFAIL(2) = No. of times residual CTRSYL3 exceeds threshold THRESH
          NFAIL(3) = No. of times CTRSYL3 and CTRSYL deviate
[out]RMAX
          RMAX is DOUBLE PRECISION array, dimension (2)
          RMAX(1) = Value of the largest test ratio of CTRSYL
          RMAX(2) = Value of the largest test ratio of CTRSYL3
[out]NINFO
          NINFO is INTEGER array, dimension (2)
          NINFO(1) = No. of times CTRSYL where INFO is nonzero
          NINFO(2) = No. of times CTRSYL3 where INFO is nonzero
[out]KNT
          KNT is INTEGER
          Total number of examples tested.

Definition at line 88 of file csyl01.f.

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*
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
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
sisnan
logical function sisnan(sin)
SISNAN tests input for NaN.
Definition sisnan.f:59
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:103
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
clange
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition clange.f:115
ctrsyl3
subroutine ctrsyl3(trana, tranb, isgn, m, n, a, lda, b, ldb, c, ldc, scale, swork, ldswork, info)
CTRSYL3
Definition ctrsyl3.f:156
ctrsyl
subroutine ctrsyl(trana, tranb, isgn, m, n, a, lda, b, ldb, c, ldc, scale, info)
CTRSYL
Definition ctrsyl.f:157
Here is the call graph for this function:
Here is the caller graph for this function: