89 IMPLICIT NONE
90
91
92
93
94
95
96 INTEGER KNT
97 REAL THRESH
98
99
100 INTEGER NFAIL( 3 ), NINFO( 2 )
101 REAL RMAX( 2 )
102
103
104
105
106
107 REAL ZERO, ONE
108 parameter( zero = 0.0e+0, one = 1.0e+0 )
109 INTEGER MAXM, MAXN, LDSWORK
110 parameter( maxm = 101, maxn = 138, ldswork = 18 )
111
112
113 CHARACTER TRANA, TRANB
114 INTEGER I, INFO, IINFO, ISGN, ITRANA, ITRANB, J, KLA,
115 $ KUA, KLB, KUB, LIWORK, M, N
116 REAL ANRM, BNRM, BIGNUM, EPS, RES, RES1, RMUL,
117 $ SCALE, SCALE3, SMLNUM, TNRM, XNRM
118
119
120 REAL DUML( MAXM ), DUMR( MAXN ),
121 $ D( MAX( MAXM, MAXN ) ), DUM( MAXN ),
122 $ VM( 2 )
123 INTEGER ISEED( 4 ), IWORK( MAXM + MAXN + 2 )
124
125
126 INTEGER AllocateStatus
127 REAL, DIMENSION(:,:), ALLOCATABLE :: A, B, C, CC, X,
128 $ SWORK
129
130
131 LOGICAL SISNAN
132 REAL SLAMCH, SLANGE
134
135
137
138
139 INTRINSIC abs, real, max
140
141
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, 54 ), stat = allocatestatus )
153 IF( allocatestatus /= 0 ) stop "*** Not enough memory ***"
154
155
156
157
158
160 smlnum =
slamch(
'S' ) / eps
161 bignum = one / smlnum
162
163 vm( 1 ) = one
164 vm( 2 ) = 0.05e+0
165
166
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
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 slatmr( 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 =
slange(
'M', m, m, a, maxm, dum )
200 DO n = 51, maxn, 47
201 klb = 0
202 kub = n - 1
203 CALL slatmr( 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 =
slange(
'M', n, n, b, maxn, dum )
209 tnrm = max( anrm, bnrm )
210 CALL slatmr( 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 slacpy(
'All', m, n, c, maxm, x, maxm)
232 CALL slacpy(
'All', m, n, c, maxm, cc, maxm)
233 CALL strsyl( 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 =
slange(
'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 sgemm( trana,
'N', m, n, m, rmul,
246 $ a, maxm, x, maxm, -scale*rmul,
247 $ c, maxm )
248 CALL sgemm(
'N', tranb, m, n, n,
249 $ real( isgn )*rmul, x, maxm, b,
250 $ maxn, one, c, maxm )
251 res1 =
slange(
'M', m, n, c, 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 slacpy(
'All', m, n, c, maxm, x, maxm )
260 CALL slacpy(
'All', m, n, c, maxm, cc, maxm )
261 CALL strsyl3( 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 =
slange(
'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 sgemm( trana,
'N', m, n, m, rmul,
275 $ a, maxm, x, maxm, -scale3*rmul,
276 $ cc, maxm )
277 CALL sgemm(
'N', tranb, m, n, n,
278 $ real( isgn )*rmul, x, maxm, b,
279 $ maxn, one, cc, maxm )
280 res1 =
slange(
'M', m, n, cc, maxm, dum )
281 res = res1 / max( smlnum, smlnum*xnrm,
282 $ ( ( rmul*tnrm )*eps )*xnrm )
283
284
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.
sisnan( 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
310
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
logical function sisnan(sin)
SISNAN tests input for NaN.
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
real function slamch(cmach)
SLAMCH
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine strsyl(trana, tranb, isgn, m, n, a, lda, b, ldb, c, ldc, scale, info)
STRSYL
subroutine slatmr(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)
SLATMR
subroutine strsyl3(trana, tranb, isgn, m, n, a, lda, b, ldb, c, ldc, scale, iwork, liwork, swork, ldswork, info)
STRSYL3