LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
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  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 *
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