LAPACK  3.10.1
LAPACK: Linear Algebra PACKage

◆ dget39()

subroutine dget39 ( double precision  RMAX,
integer  LMAX,
integer  NINFO,
integer  KNT 
)

DGET39

Purpose:
 DGET39 tests DLAQTR, 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 DTRSNA.

 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 DOUBLE PRECISION
          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 dget39.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  DOUBLE PRECISION RMAX
111 * ..
112 *
113 * =====================================================================
114 *
115 * .. Parameters ..
116  INTEGER LDT, LDT2
117  parameter( ldt = 10, ldt2 = 2*ldt )
118  DOUBLE PRECISION ZERO, ONE
119  parameter( zero = 0.0d0, one = 1.0d0 )
120 * ..
121 * .. Local Scalars ..
122  INTEGER I, INFO, IVM1, IVM2, IVM3, IVM4, IVM5, J, K, N,
123  $ NDIM
124  DOUBLE PRECISION BIGNUM, DOMIN, DUMM, EPS, NORM, NORMTB, RESID,
125  $ SCALE, SMLNUM, W, XNORM
126 * ..
127 * .. External Functions ..
128  INTEGER IDAMAX
129  DOUBLE PRECISION DASUM, DDOT, DLAMCH, DLANGE
130  EXTERNAL idamax, dasum, ddot, dlamch, dlange
131 * ..
132 * .. External Subroutines ..
133  EXTERNAL dcopy, dgemv, dlabad, dlaqtr
134 * ..
135 * .. Intrinsic Functions ..
136  INTRINSIC abs, cos, dble, max, sin, sqrt
137 * ..
138 * .. Local Arrays ..
139  INTEGER IDIM( 6 ), IVAL( 5, 5, 6 )
140  DOUBLE PRECISION 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 = dlamch( 'P' )
160  smlnum = dlamch( 'S' )
161  bignum = one / smlnum
162  CALL dlabad( 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 ) = dble( 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( dble( i ) )*vm3( ivm3 )
224  30 CONTINUE
225 *
226  DO 40 i = 1, 2*n
227  d( i ) = sin( dble( i ) )*vm4( ivm4 )
228  40 CONTINUE
229 *
230  norm = dlange( '1', n, n, t, ldt, work )
231  k = idamax( n, b, 1 )
232  normtb = norm + abs( b( k ) ) + abs( w )
233 *
234  CALL dcopy( n, d, 1, x, 1 )
235  knt = knt + 1
236  CALL dlaqtr( .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 dcopy( n, d, 1, y, 1 )
245  CALL dgemv( 'No transpose', n, n, one, t, ldt,
246  $ x, 1, -scale, y, 1 )
247  xnorm = dasum( n, x, 1 )
248  resid = dasum( 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 dcopy( n, d, 1, x, 1 )
258  knt = knt + 1
259  CALL dlaqtr( .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 dcopy( n, d, 1, y, 1 )
268  CALL dgemv( 'Transpose', n, n, one, t, ldt, x,
269  $ 1, -scale, y, 1 )
270  xnorm = dasum( n, x, 1 )
271  resid = dasum( 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 dcopy( 2*n, d, 1, x, 1 )
281  knt = knt + 1
282  CALL dlaqtr( .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 dcopy( 2*n, d, 1, y, 1 )
293  y( 1 ) = ddot( 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 dgemv( 'No transpose', n, n, one, t, ldt,
299  $ x, 1, -one, y, 1 )
300 *
301  y( 1+n ) = ddot( 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 dgemv( 'No transpose', n, n, one, t, ldt,
307  $ x( 1+n ), 1, one, y( 1+n ), 1 )
308 *
309  resid = dasum( 2*n, y, 1 )
310  domin = max( smlnum, ( smlnum / eps )*normtb,
311  $ eps*( normtb*dasum( 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 dcopy( 2*n, d, 1, x, 1 )
319  knt = knt + 1
320  CALL dlaqtr( .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 dcopy( 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 dgemv( '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 dgemv( 'Transpose', n, n, one, t, ldt,
344  $ x( 1+n ), 1, -one, y( 1+n ), 1 )
345 *
346  resid = dasum( 2*n, y, 1 )
347  domin = max( smlnum, ( smlnum / eps )*normtb,
348  $ eps*( normtb*dasum( 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 DGET39
365 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:71
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
double precision function dasum(N, DX, INCX)
DASUM
Definition: dasum.f:71
double precision function ddot(N, DX, INCX, DY, INCY)
DDOT
Definition: ddot.f:82
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:156
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: dlange.f:114
subroutine dlaqtr(LTRAN, LREAL, N, T, LDT, B, W, SCALE, X, WORK, INFO)
DLAQTR solves a real quasi-triangular system of equations, or a complex quasi-triangular system of sp...
Definition: dlaqtr.f:165
Here is the call graph for this function:
Here is the caller graph for this function: