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

◆ sget39()

subroutine sget39 ( real  RMAX,
integer  LMAX,
integer  NINFO,
integer  KNT 
)

SGET39

Purpose:
 SGET39 tests SLAQTR, a routine for solving the real or
 special complex quasi upper triangular system

      op(T)*p = scale*c,
 or
      op(T + iB)*(p+iq) = scale*(c+id),

 in real arithmetic. T is upper quasi-triangular.
 If it is complex, then the first diagonal block of T must be
 1 by 1, B has the special structure

                B = [ b(1) b(2) ... b(n) ]
                    [       w            ]
                    [           w        ]
                    [              .     ]
                    [                 w  ]

 op(A) = A or A', where A' denotes the conjugate transpose of
 the matrix A.

 On input, X = [ c ].  On output, X = [ p ].
               [ d ]                  [ q ]

 Scale is an output less than or equal to 1, chosen to avoid
 overflow in X.
 This subroutine is specially designed for the condition number
 estimation in the eigenproblem routine STRSNA.

 The test code verifies that the following residual is order 1:

      ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)||
    -----------------------------------------
        max(ulp*(||T||+||B||)*(||x1||+||x2||),
            (||T||+||B||)*smlnum/ulp,
            smlnum)

 (The (||T||+||B||)*smlnum/ulp term accounts for possible
  (gradual or nongradual) underflow in x1 and x2.)
Parameters
[out]RMAX
          RMAX is REAL
          Value of the largest test ratio.
[out]LMAX
          LMAX is INTEGER
          Example number where largest test ratio achieved.
[out]NINFO
          NINFO is INTEGER
          Number of examples where INFO is nonzero.
[out]KNT
          KNT is INTEGER
          Total number of examples tested.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 102 of file sget39.f.

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 IF( resid.GT.rmax ) 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 IF( info.NE.0 )
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 IF( resid.GT.rmax ) 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*
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:71
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 ...
Definition: slange.f:114
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
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:82
real function sdot(N, SX, INCX, SY, INCY)
SDOT
Definition: sdot.f:82
real function sasum(N, SX, INCX)
SASUM
Definition: sasum.f:72
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: