LAPACK 3.11.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, 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*
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
dlabad
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
idamax
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:71
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
dasum
double precision function dasum(N, DX, INCX)
DASUM
Definition: dasum.f:71
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:156
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: