OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
sget39.f
Go to the documentation of this file.
1*> \brief \b SGET39
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 SGET39( RMAX, LMAX, NINFO, KNT )
12*
13* .. Scalar Arguments ..
14* INTEGER KNT, LMAX, NINFO
15* REAL RMAX
16* ..
17*
18*
19*> \par Purpose:
20* =============
21*>
22*> \verbatim
23*>
24*> SGET39 tests SLAQTR, a routine for solving the real or
25*> special complex quasi upper triangular system
26*>
27*> op(T)*p = scale*c,
28*> or
29*> op(T + iB)*(p+iq) = scale*(c+id),
30*>
31*> in real arithmetic. T is upper quasi-triangular.
32*> If it is complex, then the first diagonal block of T must be
33*> 1 by 1, B has the special structure
34*>
35*> B = [ b(1) b(2) ... b(n) ]
36*> [ w ]
37*> [ w ]
38*> [ . ]
39*> [ w ]
40*>
41*> op(A) = A or A', where A' denotes the conjugate transpose of
42*> the matrix A.
43*>
44*> On input, X = [ c ]. On output, X = [ p ].
45*> [ d ] [ q ]
46*>
47*> Scale is an output less than or equal to 1, chosen to avoid
48*> overflow in X.
49*> This subroutine is specially designed for the condition number
50*> estimation in the eigenproblem routine STRSNA.
51*>
52*> The test code verifies that the following residual is order 1:
53*>
54*> ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)||
55*> -----------------------------------------
56*> max(ulp*(||T||+||B||)*(||x1||+||x2||),
57*> (||T||+||B||)*smlnum/ulp,
58*> smlnum)
59*>
60*> (The (||T||+||B||)*smlnum/ulp term accounts for possible
61*> (gradual or nongradual) underflow in x1 and x2.)
62*> \endverbatim
63*
64* Arguments:
65* ==========
66*
67*> \param[out] RMAX
68*> \verbatim
69*> RMAX is REAL
70*> Value of the largest test ratio.
71*> \endverbatim
72*>
73*> \param[out] LMAX
74*> \verbatim
75*> LMAX is INTEGER
76*> Example number where largest test ratio achieved.
77*> \endverbatim
78*>
79*> \param[out] NINFO
80*> \verbatim
81*> NINFO is INTEGER
82*> Number of examples where INFO is nonzero.
83*> \endverbatim
84*>
85*> \param[out] KNT
86*> \verbatim
87*> KNT is INTEGER
88*> Total number of examples tested.
89*> \endverbatim
90*
91* Authors:
92* ========
93*
94*> \author Univ. of Tennessee
95*> \author Univ. of California Berkeley
96*> \author Univ. of Colorado Denver
97*> \author NAG Ltd.
98*
99*> \ingroup single_eig
100*
101* =====================================================================
102 SUBROUTINE sget39( RMAX, LMAX, NINFO, KNT )
103*
104* -- LAPACK test routine --
105* -- LAPACK is a software package provided by Univ. of Tennessee, --
106* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
107*
108* .. Scalar Arguments ..
109 INTEGER KNT, LMAX, NINFO
110 REAL RMAX
111* ..
112*
113* =====================================================================
114*
115* .. Parameters ..
116 INTEGER LDT, LDT2
117 parameter( ldt = 10, ldt2 = 2*ldt )
118 REAL ZERO, ONE
119 parameter( zero = 0.0, one = 1.0 )
120* ..
121* .. Local Scalars ..
122 INTEGER I, INFO, IVM1, IVM2, IVM3, IVM4, IVM5, J, K, N,
123 $ NDIM
124 REAL BIGNUM, DOMIN, DUMM, EPS, NORM, NORMTB, RESID,
125 $ SCALE, SMLNUM, W, XNORM
126* ..
127* .. External Functions ..
128 INTEGER ISAMAX
129 REAL SASUM, SDOT, SLAMCH, SLANGE
130 EXTERNAL isamax, sasum, sdot, slamch, slange
131* ..
132* .. External Subroutines ..
133 EXTERNAL scopy, sgemv, slabad, slaqtr
134* ..
135* .. Intrinsic Functions ..
136 INTRINSIC abs, cos, max, real, sin, sqrt
137* ..
138* .. Local Arrays ..
139 INTEGER IDIM( 6 ), IVAL( 5, 5, 6 )
140 REAL B( LDT ), D( LDT2 ), DUM( 1 ), T( LDT, LDT ),
141 $ VM1( 5 ), VM2( 5 ), VM3( 5 ), VM4( 5 ),
142 $ VM5( 3 ), WORK( LDT ), X( LDT2 ), Y( LDT2 )
143* ..
144* .. Data statements ..
145 DATA idim / 4, 5*5 /
146 DATA ival / 3, 4*0, 1, 1, -1, 0, 0, 3, 2, 1, 0, 0,
147 $ 4, 3, 2, 2, 0, 5*0, 1, 4*0, 2, 2, 3*0, 3, 3, 4,
148 $ 0, 0, 4, 2, 2, 3, 0, 4*1, 5, 1, 4*0, 2, 4, -2,
149 $ 0, 0, 3, 3, 4, 0, 0, 4, 2, 2, 3, 0, 5*1, 1,
150 $ 4*0, 2, 1, -1, 0, 0, 9, 8, 1, 0, 0, 4, 9, 1, 2,
151 $ -1, 5*2, 9, 4*0, 6, 4, 0, 0, 0, 3, 2, 1, 1, 0,
152 $ 5, 1, -1, 1, 0, 5*2, 4, 4*0, 2, 2, 0, 0, 0, 1,
153 $ 4, 4, 0, 0, 2, 4, 2, 2, -1, 5*2 /
154* ..
155* .. Executable Statements ..
156*
157* Get machine parameters
158*
159 eps = slamch( 'P' )
160 smlnum = slamch( 'S' )
161 bignum = one / smlnum
162 CALL slabad( smlnum, bignum )
163*
164* Set up test case parameters
165*
166 vm1( 1 ) = one
167 vm1( 2 ) = sqrt( smlnum )
168 vm1( 3 ) = sqrt( vm1( 2 ) )
169 vm1( 4 ) = sqrt( bignum )
170 vm1( 5 ) = sqrt( vm1( 4 ) )
171*
172 vm2( 1 ) = one
173 vm2( 2 ) = sqrt( smlnum )
174 vm2( 3 ) = sqrt( vm2( 2 ) )
175 vm2( 4 ) = sqrt( bignum )
176 vm2( 5 ) = sqrt( vm2( 4 ) )
177*
178 vm3( 1 ) = one
179 vm3( 2 ) = sqrt( smlnum )
180 vm3( 3 ) = sqrt( vm3( 2 ) )
181 vm3( 4 ) = sqrt( bignum )
182 vm3( 5 ) = sqrt( vm3( 4 ) )
183*
184 vm4( 1 ) = one
185 vm4( 2 ) = sqrt( smlnum )
186 vm4( 3 ) = sqrt( vm4( 2 ) )
187 vm4( 4 ) = sqrt( bignum )
188 vm4( 5 ) = sqrt( vm4( 4 ) )
189*
190 vm5( 1 ) = one
191 vm5( 2 ) = eps
192 vm5( 3 ) = sqrt( smlnum )
193*
194* Initialization
195*
196 knt = 0
197 rmax = zero
198 ninfo = 0
199 smlnum = smlnum / eps
200*
201* Begin test loop
202*
203 DO 140 ivm5 = 1, 3
204 DO 130 ivm4 = 1, 5
205 DO 120 ivm3 = 1, 5
206 DO 110 ivm2 = 1, 5
207 DO 100 ivm1 = 1, 5
208 DO 90 ndim = 1, 6
209*
210 n = idim( ndim )
211 DO 20 i = 1, n
212 DO 10 j = 1, n
213 t( i, j ) = real( ival( i, j, ndim ) )*
214 $ vm1( ivm1 )
215 IF( i.GE.j )
216 $ t( i, j ) = t( i, j )*vm5( ivm5 )
217 10 CONTINUE
218 20 CONTINUE
219*
220 w = one*vm2( ivm2 )
221*
222 DO 30 i = 1, n
223 b( i ) = cos( real( i ) )*vm3( ivm3 )
224 30 CONTINUE
225*
226 DO 40 i = 1, 2*n
227 d( i ) = sin( real( i ) )*vm4( ivm4 )
228 40 CONTINUE
229*
230 norm = slange( '1', n, n, t, ldt, work )
231 k = isamax( n, b, 1 )
232 normtb = norm + abs( b( k ) ) + abs( w )
233*
234 CALL scopy( n, d, 1, x, 1 )
235 knt = knt + 1
236 CALL slaqtr( .false., .true., n, t, ldt, dum,
237 $ dumm, scale, x, work, info )
238 IF( info.NE.0 )
239 $ ninfo = ninfo + 1
240*
241* || T*x - scale*d || /
242* max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
243*
244 CALL scopy( n, d, 1, y, 1 )
245 CALL sgemv( 'No transpose', n, n, one, t, ldt,
246 $ x, 1, -scale, y, 1 )
247 xnorm = sasum( n, x, 1 )
248 resid = sasum( n, y, 1 )
249 domin = max( smlnum, ( smlnum / eps )*norm,
250 $ ( norm*eps )*xnorm )
251 resid = resid / domin
252 IF( resid.GT.rmax ) THEN
253 rmax = resid
254 lmax = knt
255 END IF
256*
257 CALL scopy( n, d, 1, x, 1 )
258 knt = knt + 1
259 CALL slaqtr( .true., .true., n, t, ldt, dum,
260 $ dumm, scale, x, work, info )
261 IF( info.NE.0 )
262 $ ninfo = ninfo + 1
263*
264* || T*x - scale*d || /
265* max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
266*
267 CALL scopy( n, d, 1, y, 1 )
268 CALL sgemv( 'Transpose', n, n, one, t, ldt, x,
269 $ 1, -scale, y, 1 )
270 xnorm = sasum( n, x, 1 )
271 resid = sasum( n, y, 1 )
272 domin = max( smlnum, ( smlnum / eps )*norm,
273 $ ( norm*eps )*xnorm )
274 resid = resid / domin
275 IF( resid.GT.rmax ) THEN
276 rmax = resid
277 lmax = knt
278 END IF
279*
280 CALL scopy( 2*n, d, 1, x, 1 )
281 knt = knt + 1
282 CALL slaqtr( .false., .false., n, t, ldt, b, w,
283 $ scale, x, work, info )
284 IF( info.NE.0 )
285 $ ninfo = ninfo + 1
286*
287* ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
288* max(ulp*(||T||+||B||)*(||x1||+||x2||),
289* smlnum/ulp * (||T||+||B||), smlnum )
290*
291*
292 CALL scopy( 2*n, d, 1, y, 1 )
293 y( 1 ) = sdot( n, b, 1, x( 1+n ), 1 ) +
294 $ scale*y( 1 )
295 DO 50 i = 2, n
296 y( i ) = w*x( i+n ) + scale*y( i )
297 50 CONTINUE
298 CALL sgemv( 'No transpose', n, n, one, t, ldt,
299 $ x, 1, -one, y, 1 )
300*
301 y( 1+n ) = sdot( n, b, 1, x, 1 ) -
302 $ scale*y( 1+n )
303 DO 60 i = 2, n
304 y( i+n ) = w*x( i ) - scale*y( i+n )
305 60 CONTINUE
306 CALL sgemv( 'no transpose', N, N, ONE, T, LDT,
307 $ X( 1+N ), 1, ONE, Y( 1+N ), 1 )
308*
309 RESID = SASUM( 2*N, Y, 1 )
310 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB,
311 $ EPS*( NORMTB*SASUM( 2*N, X, 1 ) ) )
312 RESID = RESID / DOMIN
313.GT. IF( RESIDRMAX ) THEN
314 RMAX = RESID
315 LMAX = KNT
316 END IF
317*
318 CALL SCOPY( 2*N, D, 1, X, 1 )
319 KNT = KNT + 1
320 CALL SLAQTR( .TRUE., .FALSE., N, T, LDT, B, W,
321 $ SCALE, X, WORK, INFO )
322.NE. IF( INFO0 )
323 $ NINFO = NINFO + 1
324*
325* ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
326* max(ulp*(||T||+||B||)*(||x1||+||x2||),
327* smlnum/ulp * (||T||+||B||), smlnum )
328*
329 CALL SCOPY( 2*N, D, 1, Y, 1 )
330 Y( 1 ) = B( 1 )*X( 1+N ) - SCALE*Y( 1 )
331 DO 70 I = 2, N
332 Y( I ) = B( I )*X( 1+N ) + W*X( I+N ) -
333 $ SCALE*Y( I )
334 70 CONTINUE
335 CALL SGEMV( 'transpose', N, N, ONE, T, LDT, X,
336 $ 1, ONE, Y, 1 )
337*
338 Y( 1+N ) = B( 1 )*X( 1 ) + SCALE*Y( 1+N )
339 DO 80 I = 2, N
340 Y( I+N ) = B( I )*X( 1 ) + W*X( I ) +
341 $ SCALE*Y( I+N )
342 80 CONTINUE
343 CALL SGEMV( 'transpose', N, N, ONE, T, LDT,
344 $ X( 1+N ), 1, -ONE, Y( 1+N ), 1 )
345*
346 RESID = SASUM( 2*N, Y, 1 )
347 DOMIN = MAX( SMLNUM, ( SMLNUM / EPS )*NORMTB,
348 $ EPS*( NORMTB*SASUM( 2*N, X, 1 ) ) )
349 RESID = RESID / DOMIN
350.GT. IF( RESIDRMAX ) THEN
351 RMAX = RESID
352 LMAX = KNT
353 END IF
354*
355 90 CONTINUE
356 100 CONTINUE
357 110 CONTINUE
358 120 CONTINUE
359 130 CONTINUE
360 140 CONTINUE
361*
362 RETURN
363*
364* End of SGET39
365*
366 END
slabad
subroutine slabad(small, large)
SLABAD
Definition slabad.f:74
slaqtr
subroutine slaqtr(ltran, lreal, n, t, ldt, b, w, scale, x, work, info)
SLAQTR solves a real quasi-triangular system of equations, or a complex quasi-triangular system of sp...
Definition slaqtr.f:165
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sgemv
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:156
sget39
subroutine sget39(rmax, lmax, ninfo, knt)
SGET39
Definition sget39.f:103
max
#define max(a, b)
Definition macros.h:21