LAPACK 3.12.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, 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*
163* Set up test case parameters
164*
165 vm1( 1 ) = one
166 vm1( 2 ) = sqrt( smlnum )
167 vm1( 3 ) = sqrt( vm1( 2 ) )
168 vm1( 4 ) = sqrt( bignum )
169 vm1( 5 ) = sqrt( vm1( 4 ) )
170*
171 vm2( 1 ) = one
172 vm2( 2 ) = sqrt( smlnum )
173 vm2( 3 ) = sqrt( vm2( 2 ) )
174 vm2( 4 ) = sqrt( bignum )
175 vm2( 5 ) = sqrt( vm2( 4 ) )
176*
177 vm3( 1 ) = one
178 vm3( 2 ) = sqrt( smlnum )
179 vm3( 3 ) = sqrt( vm3( 2 ) )
180 vm3( 4 ) = sqrt( bignum )
181 vm3( 5 ) = sqrt( vm3( 4 ) )
182*
183 vm4( 1 ) = one
184 vm4( 2 ) = sqrt( smlnum )
185 vm4( 3 ) = sqrt( vm4( 2 ) )
186 vm4( 4 ) = sqrt( bignum )
187 vm4( 5 ) = sqrt( vm4( 4 ) )
188*
189 vm5( 1 ) = one
190 vm5( 2 ) = eps
191 vm5( 3 ) = sqrt( smlnum )
192*
193* Initialization
194*
195 knt = 0
196 rmax = zero
197 ninfo = 0
198 smlnum = smlnum / eps
199*
200* Begin test loop
201*
202 DO 140 ivm5 = 1, 3
203 DO 130 ivm4 = 1, 5
204 DO 120 ivm3 = 1, 5
205 DO 110 ivm2 = 1, 5
206 DO 100 ivm1 = 1, 5
207 DO 90 ndim = 1, 6
208*
209 n = idim( ndim )
210 DO 20 i = 1, n
211 DO 10 j = 1, n
212 t( i, j ) = real( ival( i, j, ndim ) )*
213 $ vm1( ivm1 )
214 IF( i.GE.j )
215 $ t( i, j ) = t( i, j )*vm5( ivm5 )
216 10 CONTINUE
217 20 CONTINUE
218*
219 w = one*vm2( ivm2 )
220*
221 DO 30 i = 1, n
222 b( i ) = cos( real( i ) )*vm3( ivm3 )
223 30 CONTINUE
224*
225 DO 40 i = 1, 2*n
226 d( i ) = sin( real( i ) )*vm4( ivm4 )
227 40 CONTINUE
228*
229 norm = slange( '1', n, n, t, ldt, work )
230 k = isamax( n, b, 1 )
231 normtb = norm + abs( b( k ) ) + abs( w )
232*
233 CALL scopy( n, d, 1, x, 1 )
234 knt = knt + 1
235 CALL slaqtr( .false., .true., n, t, ldt, dum,
236 $ dumm, scale, x, work, info )
237 IF( info.NE.0 )
238 $ ninfo = ninfo + 1
239*
240* || T*x - scale*d || /
241* max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
242*
243 CALL scopy( n, d, 1, y, 1 )
244 CALL sgemv( 'No transpose', n, n, one, t, ldt,
245 $ x, 1, -scale, y, 1 )
246 xnorm = sasum( n, x, 1 )
247 resid = sasum( n, y, 1 )
248 domin = max( smlnum, ( smlnum / eps )*norm,
249 $ ( norm*eps )*xnorm )
250 resid = resid / domin
251 IF( resid.GT.rmax ) THEN
252 rmax = resid
253 lmax = knt
254 END IF
255*
256 CALL scopy( n, d, 1, x, 1 )
257 knt = knt + 1
258 CALL slaqtr( .true., .true., n, t, ldt, dum,
259 $ dumm, scale, x, work, info )
260 IF( info.NE.0 )
261 $ ninfo = ninfo + 1
262*
263* || T*x - scale*d || /
264* max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
265*
266 CALL scopy( n, d, 1, y, 1 )
267 CALL sgemv( 'Transpose', n, n, one, t, ldt, x,
268 $ 1, -scale, y, 1 )
269 xnorm = sasum( n, x, 1 )
270 resid = sasum( n, y, 1 )
271 domin = max( smlnum, ( smlnum / eps )*norm,
272 $ ( norm*eps )*xnorm )
273 resid = resid / domin
274 IF( resid.GT.rmax ) THEN
275 rmax = resid
276 lmax = knt
277 END IF
278*
279 CALL scopy( 2*n, d, 1, x, 1 )
280 knt = knt + 1
281 CALL slaqtr( .false., .false., n, t, ldt, b, w,
282 $ scale, x, work, info )
283 IF( info.NE.0 )
284 $ ninfo = ninfo + 1
285*
286* ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
287* max(ulp*(||T||+||B||)*(||x1||+||x2||),
288* smlnum/ulp * (||T||+||B||), smlnum )
289*
290*
291 CALL scopy( 2*n, d, 1, y, 1 )
292 y( 1 ) = sdot( n, b, 1, x( 1+n ), 1 ) +
293 $ scale*y( 1 )
294 DO 50 i = 2, n
295 y( i ) = w*x( i+n ) + scale*y( i )
296 50 CONTINUE
297 CALL sgemv( 'No transpose', n, n, one, t, ldt,
298 $ x, 1, -one, y, 1 )
299*
300 y( 1+n ) = sdot( n, b, 1, x, 1 ) -
301 $ scale*y( 1+n )
302 DO 60 i = 2, n
303 y( i+n ) = w*x( i ) - scale*y( i+n )
304 60 CONTINUE
305 CALL sgemv( 'No transpose', n, n, one, t, ldt,
306 $ x( 1+n ), 1, one, y( 1+n ), 1 )
307*
308 resid = sasum( 2*n, y, 1 )
309 domin = max( smlnum, ( smlnum / eps )*normtb,
310 $ eps*( normtb*sasum( 2*n, x, 1 ) ) )
311 resid = resid / domin
312 IF( resid.GT.rmax ) THEN
313 rmax = resid
314 lmax = knt
315 END IF
316*
317 CALL scopy( 2*n, d, 1, x, 1 )
318 knt = knt + 1
319 CALL slaqtr( .true., .false., n, t, ldt, b, w,
320 $ scale, x, work, info )
321 IF( info.NE.0 )
322 $ ninfo = ninfo + 1
323*
324* ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
325* max(ulp*(||T||+||B||)*(||x1||+||x2||),
326* smlnum/ulp * (||T||+||B||), smlnum )
327*
328 CALL scopy( 2*n, d, 1, y, 1 )
329 y( 1 ) = b( 1 )*x( 1+n ) - scale*y( 1 )
330 DO 70 i = 2, n
331 y( i ) = b( i )*x( 1+n ) + w*x( i+n ) -
332 $ scale*y( i )
333 70 CONTINUE
334 CALL sgemv( 'Transpose', n, n, one, t, ldt, x,
335 $ 1, one, y, 1 )
336*
337 y( 1+n ) = b( 1 )*x( 1 ) + scale*y( 1+n )
338 DO 80 i = 2, n
339 y( i+n ) = b( i )*x( 1 ) + w*x( i ) +
340 $ scale*y( i+n )
341 80 CONTINUE
342 CALL sgemv( 'Transpose', n, n, one, t, ldt,
343 $ x( 1+n ), 1, -one, y( 1+n ), 1 )
344*
345 resid = sasum( 2*n, y, 1 )
346 domin = max( smlnum, ( smlnum / eps )*normtb,
347 $ eps*( normtb*sasum( 2*n, x, 1 ) ) )
348 resid = resid / domin
349 IF( resid.GT.rmax ) THEN
350 rmax = resid
351 lmax = knt
352 END IF
353*
354 90 CONTINUE
355 100 CONTINUE
356 110 CONTINUE
357 120 CONTINUE
358 130 CONTINUE
359 140 CONTINUE
360*
361 RETURN
362*
363* End of SGET39
364*
sasum
real function sasum(n, sx, incx)
SASUM
Definition sasum.f:72
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
sdot
real function sdot(n, sx, incx, sy, incy)
SDOT
Definition sdot.f:82
sgemv
subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
SGEMV
Definition sgemv.f:158
isamax
integer function isamax(n, sx, incx)
ISAMAX
Definition isamax.f:71
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
slange
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
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
Here is the call graph for this function:
Here is the caller graph for this function: