LAPACK  3.7.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 *> \date December 2016
100 *
101 *> \ingroup single_eig
102 *
103 * =====================================================================
104  SUBROUTINE sget39( RMAX, LMAX, NINFO, KNT )
105 *
106 * -- LAPACK test routine (version 3.7.0) --
107 * -- LAPACK is a software package provided by Univ. of Tennessee, --
108 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
109 * December 2016
110 *
111 * .. Scalar Arguments ..
112  INTEGER KNT, LMAX, NINFO
113  REAL RMAX
114 * ..
115 *
116 * =====================================================================
117 *
118 * .. Parameters ..
119  INTEGER LDT, LDT2
120  parameter( ldt = 10, ldt2 = 2*ldt )
121  REAL ZERO, ONE
122  parameter( zero = 0.0, one = 1.0 )
123 * ..
124 * .. Local Scalars ..
125  INTEGER I, INFO, IVM1, IVM2, IVM3, IVM4, IVM5, J, K, N,
126  $ ndim
127  REAL BIGNUM, DOMIN, DUMM, EPS, NORM, NORMTB, RESID,
128  $ scale, smlnum, w, xnorm
129 * ..
130 * .. External Functions ..
131  INTEGER ISAMAX
132  REAL SASUM, SDOT, SLAMCH, SLANGE
133  EXTERNAL isamax, sasum, sdot, slamch, slange
134 * ..
135 * .. External Subroutines ..
136  EXTERNAL scopy, sgemv, slabad, slaqtr
137 * ..
138 * .. Intrinsic Functions ..
139  INTRINSIC abs, cos, max, REAL, SIN, SQRT
140 * ..
141 * .. Local Arrays ..
142  INTEGER IDIM( 6 ), IVAL( 5, 5, 6 )
143  REAL B( ldt ), D( ldt2 ), DUM( 1 ), T( ldt, ldt ),
144  $ vm1( 5 ), vm2( 5 ), vm3( 5 ), vm4( 5 ),
145  $ vm5( 3 ), work( ldt ), x( ldt2 ), y( ldt2 )
146 * ..
147 * .. Data statements ..
148  DATA idim / 4, 5*5 /
149  DATA ival / 3, 4*0, 1, 1, -1, 0, 0, 3, 2, 1, 0, 0,
150  $ 4, 3, 2, 2, 0, 5*0, 1, 4*0, 2, 2, 3*0, 3, 3, 4,
151  $ 0, 0, 4, 2, 2, 3, 0, 4*1, 5, 1, 4*0, 2, 4, -2,
152  $ 0, 0, 3, 3, 4, 0, 0, 4, 2, 2, 3, 0, 5*1, 1,
153  $ 4*0, 2, 1, -1, 0, 0, 9, 8, 1, 0, 0, 4, 9, 1, 2,
154  $ -1, 5*2, 9, 4*0, 6, 4, 0, 0, 0, 3, 2, 1, 1, 0,
155  $ 5, 1, -1, 1, 0, 5*2, 4, 4*0, 2, 2, 0, 0, 0, 1,
156  $ 4, 4, 0, 0, 2, 4, 2, 2, -1, 5*2 /
157 * ..
158 * .. Executable Statements ..
159 *
160 * Get machine parameters
161 *
162  eps = slamch( 'P' )
163  smlnum = slamch( 'S' )
164  bignum = one / smlnum
165  CALL slabad( smlnum, bignum )
166 *
167 * Set up test case parameters
168 *
169  vm1( 1 ) = one
170  vm1( 2 ) = sqrt( smlnum )
171  vm1( 3 ) = sqrt( vm1( 2 ) )
172  vm1( 4 ) = sqrt( bignum )
173  vm1( 5 ) = sqrt( vm1( 4 ) )
174 *
175  vm2( 1 ) = one
176  vm2( 2 ) = sqrt( smlnum )
177  vm2( 3 ) = sqrt( vm2( 2 ) )
178  vm2( 4 ) = sqrt( bignum )
179  vm2( 5 ) = sqrt( vm2( 4 ) )
180 *
181  vm3( 1 ) = one
182  vm3( 2 ) = sqrt( smlnum )
183  vm3( 3 ) = sqrt( vm3( 2 ) )
184  vm3( 4 ) = sqrt( bignum )
185  vm3( 5 ) = sqrt( vm3( 4 ) )
186 *
187  vm4( 1 ) = one
188  vm4( 2 ) = sqrt( smlnum )
189  vm4( 3 ) = sqrt( vm4( 2 ) )
190  vm4( 4 ) = sqrt( bignum )
191  vm4( 5 ) = sqrt( vm4( 4 ) )
192 *
193  vm5( 1 ) = one
194  vm5( 2 ) = eps
195  vm5( 3 ) = sqrt( smlnum )
196 *
197 * Initalization
198 *
199  knt = 0
200  rmax = zero
201  ninfo = 0
202  smlnum = smlnum / eps
203 *
204 * Begin test loop
205 *
206  DO 140 ivm5 = 1, 3
207  DO 130 ivm4 = 1, 5
208  DO 120 ivm3 = 1, 5
209  DO 110 ivm2 = 1, 5
210  DO 100 ivm1 = 1, 5
211  DO 90 ndim = 1, 6
212 *
213  n = idim( ndim )
214  DO 20 i = 1, n
215  DO 10 j = 1, n
216  t( i, j ) = REAL( IVAL( I, J, NDIM ) )*
217  $ vm1( ivm1 )
218  IF( i.GE.j )
219  $ t( i, j ) = t( i, j )*vm5( ivm5 )
220  10 CONTINUE
221  20 CONTINUE
222 *
223  w = one*vm2( ivm2 )
224 *
225  DO 30 i = 1, n
226  b( i ) = cos( REAL( I ) )*vm3( IVM3 )
227  30 CONTINUE
228 *
229  DO 40 i = 1, 2*n
230  d( i ) = sin( REAL( I ) )*vm4( IVM4 )
231  40 CONTINUE
232 *
233  norm = slange( '1', n, n, t, ldt, work )
234  k = isamax( n, b, 1 )
235  normtb = norm + abs( b( k ) ) + abs( w )
236 *
237  CALL scopy( n, d, 1, x, 1 )
238  knt = knt + 1
239  CALL slaqtr( .false., .true., n, t, ldt, dum,
240  $ dumm, scale, x, work, info )
241  IF( info.NE.0 )
242  $ ninfo = ninfo + 1
243 *
244 * || T*x - scale*d || /
245 * max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
246 *
247  CALL scopy( n, d, 1, y, 1 )
248  CALL sgemv( 'No transpose', n, n, one, t, ldt,
249  $ x, 1, -scale, y, 1 )
250  xnorm = sasum( n, x, 1 )
251  resid = sasum( n, y, 1 )
252  domin = max( smlnum, ( smlnum / eps )*norm,
253  $ ( norm*eps )*xnorm )
254  resid = resid / domin
255  IF( resid.GT.rmax ) THEN
256  rmax = resid
257  lmax = knt
258  END IF
259 *
260  CALL scopy( n, d, 1, x, 1 )
261  knt = knt + 1
262  CALL slaqtr( .true., .true., n, t, ldt, dum,
263  $ dumm, scale, x, work, info )
264  IF( info.NE.0 )
265  $ ninfo = ninfo + 1
266 *
267 * || T*x - scale*d || /
268 * max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
269 *
270  CALL scopy( n, d, 1, y, 1 )
271  CALL sgemv( 'Transpose', n, n, one, t, ldt, x,
272  $ 1, -scale, y, 1 )
273  xnorm = sasum( n, x, 1 )
274  resid = sasum( n, y, 1 )
275  domin = max( smlnum, ( smlnum / eps )*norm,
276  $ ( norm*eps )*xnorm )
277  resid = resid / domin
278  IF( resid.GT.rmax ) THEN
279  rmax = resid
280  lmax = knt
281  END IF
282 *
283  CALL scopy( 2*n, d, 1, x, 1 )
284  knt = knt + 1
285  CALL slaqtr( .false., .false., n, t, ldt, b, w,
286  $ scale, x, work, info )
287  IF( info.NE.0 )
288  $ ninfo = ninfo + 1
289 *
290 * ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
291 * max(ulp*(||T||+||B||)*(||x1||+||x2||),
292 * smlnum/ulp * (||T||+||B||), smlnum )
293 *
294 *
295  CALL scopy( 2*n, d, 1, y, 1 )
296  y( 1 ) = sdot( n, b, 1, x( 1+n ), 1 ) +
297  $ scale*y( 1 )
298  DO 50 i = 2, n
299  y( i ) = w*x( i+n ) + scale*y( i )
300  50 CONTINUE
301  CALL sgemv( 'No transpose', n, n, one, t, ldt,
302  $ x, 1, -one, y, 1 )
303 *
304  y( 1+n ) = sdot( n, b, 1, x, 1 ) -
305  $ scale*y( 1+n )
306  DO 60 i = 2, n
307  y( i+n ) = w*x( i ) - scale*y( i+n )
308  60 CONTINUE
309  CALL sgemv( 'No transpose', n, n, one, t, ldt,
310  $ x( 1+n ), 1, one, y( 1+n ), 1 )
311 *
312  resid = sasum( 2*n, y, 1 )
313  domin = max( smlnum, ( smlnum / eps )*normtb,
314  $ eps*( normtb*sasum( 2*n, x, 1 ) ) )
315  resid = resid / domin
316  IF( resid.GT.rmax ) THEN
317  rmax = resid
318  lmax = knt
319  END IF
320 *
321  CALL scopy( 2*n, d, 1, x, 1 )
322  knt = knt + 1
323  CALL slaqtr( .true., .false., n, t, ldt, b, w,
324  $ scale, x, work, info )
325  IF( info.NE.0 )
326  $ ninfo = ninfo + 1
327 *
328 * ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
329 * max(ulp*(||T||+||B||)*(||x1||+||x2||),
330 * smlnum/ulp * (||T||+||B||), smlnum )
331 *
332  CALL scopy( 2*n, d, 1, y, 1 )
333  y( 1 ) = b( 1 )*x( 1+n ) - scale*y( 1 )
334  DO 70 i = 2, n
335  y( i ) = b( i )*x( 1+n ) + w*x( i+n ) -
336  $ scale*y( i )
337  70 CONTINUE
338  CALL sgemv( 'Transpose', n, n, one, t, ldt, x,
339  $ 1, one, y, 1 )
340 *
341  y( 1+n ) = b( 1 )*x( 1 ) + scale*y( 1+n )
342  DO 80 i = 2, n
343  y( i+n ) = b( i )*x( 1 ) + w*x( i ) +
344  $ scale*y( i+n )
345  80 CONTINUE
346  CALL sgemv( 'Transpose', n, n, one, t, ldt,
347  $ x( 1+n ), 1, -one, y( 1+n ), 1 )
348 *
349  resid = sasum( 2*n, y, 1 )
350  domin = max( smlnum, ( smlnum / eps )*normtb,
351  $ eps*( normtb*sasum( 2*n, x, 1 ) ) )
352  resid = resid / domin
353  IF( resid.GT.rmax ) THEN
354  rmax = resid
355  lmax = knt
356  END IF
357 *
358  90 CONTINUE
359  100 CONTINUE
360  110 CONTINUE
361  120 CONTINUE
362  130 CONTINUE
363  140 CONTINUE
364 *
365  RETURN
366 *
367 * End of SGET39
368 *
369  END
sgemv
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:158
sget39
subroutine sget39(RMAX, LMAX, NINFO, KNT)
SGET39
Definition: sget39.f:105
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
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:167
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:84