LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ 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, 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*
163* Set up test case parameters
164*
165 vm1( 1 ) = one
166 vm1( 2 ) = sqrt( smlnum )
167 vm1( 3 ) = sqrt( vm1( 2 ) )
168 vm1( 4 ) = sqrt( bignum )
169 vm1( 5 ) = sqrt( vm1( 4 ) )
170*
171 vm2( 1 ) = one
172 vm2( 2 ) = sqrt( smlnum )
173 vm2( 3 ) = sqrt( vm2( 2 ) )
174 vm2( 4 ) = sqrt( bignum )
175 vm2( 5 ) = sqrt( vm2( 4 ) )
176*
177 vm3( 1 ) = one
178 vm3( 2 ) = sqrt( smlnum )
179 vm3( 3 ) = sqrt( vm3( 2 ) )
180 vm3( 4 ) = sqrt( bignum )
181 vm3( 5 ) = sqrt( vm3( 4 ) )
182*
183 vm4( 1 ) = one
184 vm4( 2 ) = sqrt( smlnum )
185 vm4( 3 ) = sqrt( vm4( 2 ) )
186 vm4( 4 ) = sqrt( bignum )
187 vm4( 5 ) = sqrt( vm4( 4 ) )
188*
189 vm5( 1 ) = one
190 vm5( 2 ) = eps
191 vm5( 3 ) = sqrt( smlnum )
192*
193* Initialization
194*
195 knt = 0
196 rmax = zero
197 ninfo = 0
198 smlnum = smlnum / eps
199*
200* Begin test loop
201*
202 DO 140 ivm5 = 1, 3
203 DO 130 ivm4 = 1, 5
204 DO 120 ivm3 = 1, 5
205 DO 110 ivm2 = 1, 5
206 DO 100 ivm1 = 1, 5
207 DO 90 ndim = 1, 6
208*
209 n = idim( ndim )
210 DO 20 i = 1, n
211 DO 10 j = 1, n
212 t( i, j ) = dble( ival( i, j, ndim ) )*
213 $ vm1( ivm1 )
214 IF( i.GE.j )
215 $ t( i, j ) = t( i, j )*vm5( ivm5 )
216 10 CONTINUE
217 20 CONTINUE
218*
219 w = one*vm2( ivm2 )
220*
221 DO 30 i = 1, n
222 b( i ) = cos( dble( i ) )*vm3( ivm3 )
223 30 CONTINUE
224*
225 DO 40 i = 1, 2*n
226 d( i ) = sin( dble( i ) )*vm4( ivm4 )
227 40 CONTINUE
228*
229 norm = dlange( '1', n, n, t, ldt, work )
230 k = idamax( n, b, 1 )
231 normtb = norm + abs( b( k ) ) + abs( w )
232*
233 CALL dcopy( n, d, 1, x, 1 )
234 knt = knt + 1
235 CALL dlaqtr( .false., .true., n, t, ldt, dum,
236 $ dumm, scale, x, work, info )
237 IF( info.NE.0 )
238 $ ninfo = ninfo + 1
239*
240* || T*x - scale*d || /
241* max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
242*
243 CALL dcopy( n, d, 1, y, 1 )
244 CALL dgemv( 'No transpose', n, n, one, t, ldt,
245 $ x, 1, -scale, y, 1 )
246 xnorm = dasum( n, x, 1 )
247 resid = dasum( n, y, 1 )
248 domin = max( smlnum, ( smlnum / eps )*norm,
249 $ ( norm*eps )*xnorm )
250 resid = resid / domin
251 IF( resid.GT.rmax ) THEN
252 rmax = resid
253 lmax = knt
254 END IF
255*
256 CALL dcopy( n, d, 1, x, 1 )
257 knt = knt + 1
258 CALL dlaqtr( .true., .true., n, t, ldt, dum,
259 $ dumm, scale, x, work, info )
260 IF( info.NE.0 )
261 $ ninfo = ninfo + 1
262*
263* || T*x - scale*d || /
264* max(ulp*||T||*||x||,smlnum/ulp*||T||,smlnum)
265*
266 CALL dcopy( n, d, 1, y, 1 )
267 CALL dgemv( 'Transpose', n, n, one, t, ldt, x,
268 $ 1, -scale, y, 1 )
269 xnorm = dasum( n, x, 1 )
270 resid = dasum( n, y, 1 )
271 domin = max( smlnum, ( smlnum / eps )*norm,
272 $ ( norm*eps )*xnorm )
273 resid = resid / domin
274 IF( resid.GT.rmax ) THEN
275 rmax = resid
276 lmax = knt
277 END IF
278*
279 CALL dcopy( 2*n, d, 1, x, 1 )
280 knt = knt + 1
281 CALL dlaqtr( .false., .false., n, t, ldt, b, w,
282 $ scale, x, work, info )
283 IF( info.NE.0 )
284 $ ninfo = ninfo + 1
285*
286* ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
287* max(ulp*(||T||+||B||)*(||x1||+||x2||),
288* smlnum/ulp * (||T||+||B||), smlnum )
289*
290*
291 CALL dcopy( 2*n, d, 1, y, 1 )
292 y( 1 ) = ddot( n, b, 1, x( 1+n ), 1 ) +
293 $ scale*y( 1 )
294 DO 50 i = 2, n
295 y( i ) = w*x( i+n ) + scale*y( i )
296 50 CONTINUE
297 CALL dgemv( 'No transpose', n, n, one, t, ldt,
298 $ x, 1, -one, y, 1 )
299*
300 y( 1+n ) = ddot( n, b, 1, x, 1 ) -
301 $ scale*y( 1+n )
302 DO 60 i = 2, n
303 y( i+n ) = w*x( i ) - scale*y( i+n )
304 60 CONTINUE
305 CALL dgemv( 'No transpose', n, n, one, t, ldt,
306 $ x( 1+n ), 1, one, y( 1+n ), 1 )
307*
308 resid = dasum( 2*n, y, 1 )
309 domin = max( smlnum, ( smlnum / eps )*normtb,
310 $ eps*( normtb*dasum( 2*n, x, 1 ) ) )
311 resid = resid / domin
312 IF( resid.GT.rmax ) THEN
313 rmax = resid
314 lmax = knt
315 END IF
316*
317 CALL dcopy( 2*n, d, 1, x, 1 )
318 knt = knt + 1
319 CALL dlaqtr( .true., .false., n, t, ldt, b, w,
320 $ scale, x, work, info )
321 IF( info.NE.0 )
322 $ ninfo = ninfo + 1
323*
324* ||(T+i*B)*(x1+i*x2) - scale*(d1+i*d2)|| /
325* max(ulp*(||T||+||B||)*(||x1||+||x2||),
326* smlnum/ulp * (||T||+||B||), smlnum )
327*
328 CALL dcopy( 2*n, d, 1, y, 1 )
329 y( 1 ) = b( 1 )*x( 1+n ) - scale*y( 1 )
330 DO 70 i = 2, n
331 y( i ) = b( i )*x( 1+n ) + w*x( i+n ) -
332 $ scale*y( i )
333 70 CONTINUE
334 CALL dgemv( 'Transpose', n, n, one, t, ldt, x,
335 $ 1, one, y, 1 )
336*
337 y( 1+n ) = b( 1 )*x( 1 ) + scale*y( 1+n )
338 DO 80 i = 2, n
339 y( i+n ) = b( i )*x( 1 ) + w*x( i ) +
340 $ scale*y( i+n )
341 80 CONTINUE
342 CALL dgemv( 'Transpose', n, n, one, t, ldt,
343 $ x( 1+n ), 1, -one, y( 1+n ), 1 )
344*
345 resid = dasum( 2*n, y, 1 )
346 domin = max( smlnum, ( smlnum / eps )*normtb,
347 $ eps*( normtb*dasum( 2*n, x, 1 ) ) )
348 resid = resid / domin
349 IF( resid.GT.rmax ) THEN
350 rmax = resid
351 lmax = knt
352 END IF
353*
354 90 CONTINUE
355 100 CONTINUE
356 110 CONTINUE
357 120 CONTINUE
358 130 CONTINUE
359 140 CONTINUE
360*
361 RETURN
362*
363* End of DGET39
364*
dasum
double precision function dasum(n, dx, incx)
DASUM
Definition dasum.f:71
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
ddot
double precision function ddot(n, dx, incx, dy, incy)
DDOT
Definition ddot.f:82
dgemv
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlange
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
dlaqtr
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: