LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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: