LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
November 2011

Definition at line 105 of file sget39.f.

105 *
106 * -- LAPACK test routine (version 3.4.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 * November 2011
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 *
sdot
real function sdot(N, SX, INCX, SY, INCY)
SDOT
Definition: sdot.f:53
isamax
integer function isamax(N, SX, INCX)
ISAMAX
Definition: isamax.f:53
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
sgemv
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:158
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:116
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
sasum
real function sasum(N, SX, INCX)
SASUM
Definition: sasum.f:54
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53

Here is the call graph for this function:

Here is the caller graph for this function: