OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
ddrgvx.f
Go to the documentation of this file.
1*> \brief \b DDRGVX
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8* Definition:
9* ===========
10*
11* SUBROUTINE DDRGVX( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI,
12* ALPHAR, ALPHAI, BETA, VL, VR, ILO, IHI, LSCALE,
13* RSCALE, S, DTRU, DIF, DIFTRU, WORK, LWORK,
14* IWORK, LIWORK, RESULT, BWORK, INFO )
15*
16* .. Scalar Arguments ..
17* INTEGER IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT,
18* $ NSIZE
19* DOUBLE PRECISION THRESH
20* ..
21* .. Array Arguments ..
22* LOGICAL BWORK( * )
23* INTEGER IWORK( * )
24* DOUBLE PRECISION A( LDA, * ), AI( LDA, * ), ALPHAI( * ),
25* $ ALPHAR( * ), B( LDA, * ), BETA( * ),
26* $ BI( LDA, * ), DIF( * ), DIFTRU( * ), DTRU( * ),
27* $ LSCALE( * ), RESULT( 4 ), RSCALE( * ), S( * ),
28* $ VL( LDA, * ), VR( LDA, * ), WORK( * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> DDRGVX checks the nonsymmetric generalized eigenvalue problem
38*> expert driver DGGEVX.
39*>
40*> DGGEVX computes the generalized eigenvalues, (optionally) the left
41*> and/or right eigenvectors, (optionally) computes a balancing
42*> transformation to improve the conditioning, and (optionally)
43*> reciprocal condition numbers for the eigenvalues and eigenvectors.
44*>
45*> When DDRGVX is called with NSIZE > 0, two types of test matrix pairs
46*> are generated by the subroutine DLATM6 and test the driver DGGEVX.
47*> The test matrices have the known exact condition numbers for
48*> eigenvalues. For the condition numbers of the eigenvectors
49*> corresponding the first and last eigenvalues are also know
50*> ``exactly'' (see DLATM6).
51*>
52*> For each matrix pair, the following tests will be performed and
53*> compared with the threshold THRESH.
54*>
55*> (1) max over all left eigenvalue/-vector pairs (beta/alpha,l) of
56*>
57*> | l**H * (beta A - alpha B) | / ( ulp max( |beta A|, |alpha B| ) )
58*>
59*> where l**H is the conjugate tranpose of l.
60*>
61*> (2) max over all right eigenvalue/-vector pairs (beta/alpha,r) of
62*>
63*> | (beta A - alpha B) r | / ( ulp max( |beta A|, |alpha B| ) )
64*>
65*> (3) The condition number S(i) of eigenvalues computed by DGGEVX
66*> differs less than a factor THRESH from the exact S(i) (see
67*> DLATM6).
68*>
69*> (4) DIF(i) computed by DTGSNA differs less than a factor 10*THRESH
70*> from the exact value (for the 1st and 5th vectors only).
71*>
72*> Test Matrices
73*> =============
74*>
75*> Two kinds of test matrix pairs
76*>
77*> (A, B) = inverse(YH) * (Da, Db) * inverse(X)
78*>
79*> are used in the tests:
80*>
81*> 1: Da = 1+a 0 0 0 0 Db = 1 0 0 0 0
82*> 0 2+a 0 0 0 0 1 0 0 0
83*> 0 0 3+a 0 0 0 0 1 0 0
84*> 0 0 0 4+a 0 0 0 0 1 0
85*> 0 0 0 0 5+a , 0 0 0 0 1 , and
86*>
87*> 2: Da = 1 -1 0 0 0 Db = 1 0 0 0 0
88*> 1 1 0 0 0 0 1 0 0 0
89*> 0 0 1 0 0 0 0 1 0 0
90*> 0 0 0 1+a 1+b 0 0 0 1 0
91*> 0 0 0 -1-b 1+a , 0 0 0 0 1 .
92*>
93*> In both cases the same inverse(YH) and inverse(X) are used to compute
94*> (A, B), giving the exact eigenvectors to (A,B) as (YH, X):
95*>
96*> YH: = 1 0 -y y -y X = 1 0 -x -x x
97*> 0 1 -y y -y 0 1 x -x -x
98*> 0 0 1 0 0 0 0 1 0 0
99*> 0 0 0 1 0 0 0 0 1 0
100*> 0 0 0 0 1, 0 0 0 0 1 , where
101*>
102*> a, b, x and y will have all values independently of each other from
103*> { sqrt(sqrt(ULP)), 0.1, 1, 10, 1/sqrt(sqrt(ULP)) }.
104*> \endverbatim
105*
106* Arguments:
107* ==========
108*
109*> \param[in] NSIZE
110*> \verbatim
111*> NSIZE is INTEGER
112*> The number of sizes of matrices to use. NSIZE must be at
113*> least zero. If it is zero, no randomly generated matrices
114*> are tested, but any test matrices read from NIN will be
115*> tested.
116*> \endverbatim
117*>
118*> \param[in] THRESH
119*> \verbatim
120*> THRESH is DOUBLE PRECISION
121*> A test will count as "failed" if the "error", computed as
122*> described above, exceeds THRESH. Note that the error
123*> is scaled to be O(1), so THRESH should be a reasonably
124*> small multiple of 1, e.g., 10 or 100. In particular,
125*> it should not depend on the precision (single vs. double)
126*> or the size of the matrix. It must be at least zero.
127*> \endverbatim
128*>
129*> \param[in] NIN
130*> \verbatim
131*> NIN is INTEGER
132*> The FORTRAN unit number for reading in the data file of
133*> problems to solve.
134*> \endverbatim
135*>
136*> \param[in] NOUT
137*> \verbatim
138*> NOUT is INTEGER
139*> The FORTRAN unit number for printing out error messages
140*> (e.g., if a routine returns IINFO not equal to 0.)
141*> \endverbatim
142*>
143*> \param[out] A
144*> \verbatim
145*> A is DOUBLE PRECISION array, dimension (LDA, NSIZE)
146*> Used to hold the matrix whose eigenvalues are to be
147*> computed. On exit, A contains the last matrix actually used.
148*> \endverbatim
149*>
150*> \param[in] LDA
151*> \verbatim
152*> LDA is INTEGER
153*> The leading dimension of A, B, AI, BI, Ao, and Bo.
154*> It must be at least 1 and at least NSIZE.
155*> \endverbatim
156*>
157*> \param[out] B
158*> \verbatim
159*> B is DOUBLE PRECISION array, dimension (LDA, NSIZE)
160*> Used to hold the matrix whose eigenvalues are to be
161*> computed. On exit, B contains the last matrix actually used.
162*> \endverbatim
163*>
164*> \param[out] AI
165*> \verbatim
166*> AI is DOUBLE PRECISION array, dimension (LDA, NSIZE)
167*> Copy of A, modified by DGGEVX.
168*> \endverbatim
169*>
170*> \param[out] BI
171*> \verbatim
172*> BI is DOUBLE PRECISION array, dimension (LDA, NSIZE)
173*> Copy of B, modified by DGGEVX.
174*> \endverbatim
175*>
176*> \param[out] ALPHAR
177*> \verbatim
178*> ALPHAR is DOUBLE PRECISION array, dimension (NSIZE)
179*> \endverbatim
180*>
181*> \param[out] ALPHAI
182*> \verbatim
183*> ALPHAI is DOUBLE PRECISION array, dimension (NSIZE)
184*> \endverbatim
185*>
186*> \param[out] BETA
187*> \verbatim
188*> BETA is DOUBLE PRECISION array, dimension (NSIZE)
189*>
190*> On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues.
191*> \endverbatim
192*>
193*> \param[out] VL
194*> \verbatim
195*> VL is DOUBLE PRECISION array, dimension (LDA, NSIZE)
196*> VL holds the left eigenvectors computed by DGGEVX.
197*> \endverbatim
198*>
199*> \param[out] VR
200*> \verbatim
201*> VR is DOUBLE PRECISION array, dimension (LDA, NSIZE)
202*> VR holds the right eigenvectors computed by DGGEVX.
203*> \endverbatim
204*>
205*> \param[out] ILO
206*> \verbatim
207*> ILO is INTEGER
208*> \endverbatim
209*>
210*> \param[out] IHI
211*> \verbatim
212*> IHI is INTEGER
213*> \endverbatim
214*>
215*> \param[out] LSCALE
216*> \verbatim
217*> LSCALE is DOUBLE PRECISION array, dimension (N)
218*> \endverbatim
219*>
220*> \param[out] RSCALE
221*> \verbatim
222*> RSCALE is DOUBLE PRECISION array, dimension (N)
223*> \endverbatim
224*>
225*> \param[out] S
226*> \verbatim
227*> S is DOUBLE PRECISION array, dimension (N)
228*> \endverbatim
229*>
230*> \param[out] DTRU
231*> \verbatim
232*> DTRU is DOUBLE PRECISION array, dimension (N)
233*> \endverbatim
234*>
235*> \param[out] DIF
236*> \verbatim
237*> DIF is DOUBLE PRECISION array, dimension (N)
238*> \endverbatim
239*>
240*> \param[out] DIFTRU
241*> \verbatim
242*> DIFTRU is DOUBLE PRECISION array, dimension (N)
243*> \endverbatim
244*>
245*> \param[out] WORK
246*> \verbatim
247*> WORK is DOUBLE PRECISION array, dimension (LWORK)
248*> \endverbatim
249*>
250*> \param[in] LWORK
251*> \verbatim
252*> LWORK is INTEGER
253*> Leading dimension of WORK. LWORK >= 2*N*N+12*N+16.
254*> \endverbatim
255*>
256*> \param[out] IWORK
257*> \verbatim
258*> IWORK is INTEGER array, dimension (LIWORK)
259*> \endverbatim
260*>
261*> \param[in] LIWORK
262*> \verbatim
263*> LIWORK is INTEGER
264*> Leading dimension of IWORK. Must be at least N+6.
265*> \endverbatim
266*>
267*> \param[out] RESULT
268*> \verbatim
269*> RESULT is DOUBLE PRECISION array, dimension (4)
270*> \endverbatim
271*>
272*> \param[out] BWORK
273*> \verbatim
274*> BWORK is LOGICAL array, dimension (N)
275*> \endverbatim
276*>
277*> \param[out] INFO
278*> \verbatim
279*> INFO is INTEGER
280*> = 0: successful exit
281*> < 0: if INFO = -i, the i-th argument had an illegal value.
282*> > 0: A routine returned an error code.
283*> \endverbatim
284*
285* Authors:
286* ========
287*
288*> \author Univ. of Tennessee
289*> \author Univ. of California Berkeley
290*> \author Univ. of Colorado Denver
291*> \author NAG Ltd.
292*
293*> \ingroup double_eig
294*
295* =====================================================================
296 SUBROUTINE ddrgvx( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI,
297 $ ALPHAR, ALPHAI, BETA, VL, VR, ILO, IHI, LSCALE,
298 $ RSCALE, S, DTRU, DIF, DIFTRU, WORK, LWORK,
299 $ IWORK, LIWORK, RESULT, BWORK, INFO )
300*
301* -- LAPACK test routine --
302* -- LAPACK is a software package provided by Univ. of Tennessee, --
303* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
304*
305* .. Scalar Arguments ..
306 INTEGER IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT,
307 $ NSIZE
308 DOUBLE PRECISION THRESH
309* ..
310* .. Array Arguments ..
311 LOGICAL BWORK( * )
312 INTEGER IWORK( * )
313 DOUBLE PRECISION A( LDA, * ), AI( LDA, * ), ALPHAI( * ),
314 $ alphar( * ), b( lda, * ), beta( * ),
315 $ bi( lda, * ), dif( * ), diftru( * ), dtru( * ),
316 $ lscale( * ), result( 4 ), rscale( * ), s( * ),
317 $ vl( lda, * ), vr( lda, * ), work( * )
318* ..
319*
320* =====================================================================
321*
322* .. Parameters ..
323 DOUBLE PRECISION ZERO, ONE, TEN, TNTH, HALF
324 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0, ten = 1.0d+1,
325 $ tnth = 1.0d-1, half = 0.5d+0 )
326* ..
327* .. Local Scalars ..
328 INTEGER I, IPTYPE, IWA, IWB, IWX, IWY, J, LINFO,
329 $ MAXWRK, MINWRK, N, NERRS, NMAX, NPTKNT, NTESTT
330 DOUBLE PRECISION ABNORM, ANORM, BNORM, RATIO1, RATIO2, THRSH2,
331 $ ulp, ulpinv
332* ..
333* .. Local Arrays ..
334 DOUBLE PRECISION WEIGHT( 5 )
335* ..
336* .. External Functions ..
337 INTEGER ILAENV
338 DOUBLE PRECISION DLAMCH, DLANGE
339 EXTERNAL ILAENV, DLAMCH, DLANGE
340* ..
341* .. External Subroutines ..
342 EXTERNAL alasvm, dget52, dggevx, dlacpy, dlatm6, xerbla
343* ..
344* .. Intrinsic Functions ..
345 INTRINSIC abs, max, sqrt
346* ..
347* .. Executable Statements ..
348*
349* Check for errors
350*
351 info = 0
352*
353 nmax = 5
354*
355 IF( nsize.LT.0 ) THEN
356 info = -1
357 ELSE IF( thresh.LT.zero ) THEN
358 info = -2
359 ELSE IF( nin.LE.0 ) THEN
360 info = -3
361 ELSE IF( nout.LE.0 ) THEN
362 info = -4
363 ELSE IF( lda.LT.1 .OR. lda.LT.nmax ) THEN
364 info = -6
365 ELSE IF( liwork.LT.nmax+6 ) THEN
366 info = -26
367 END IF
368*
369* Compute workspace
370* (Note: Comments in the code beginning "Workspace:" describe the
371* minimal amount of workspace needed at that point in the code,
372* as well as the preferred amount for good performance.
373* NB refers to the optimal block size for the immediately
374* following subroutine, as returned by ILAENV.)
375*
376 minwrk = 1
377 IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
378 minwrk = 2*nmax*nmax + 12*nmax + 16
379 maxwrk = 6*nmax + nmax*ilaenv( 1, 'DGEQRF', ' ', nmax, 1, nmax,
380 $ 0 )
381 maxwrk = max( maxwrk, 2*nmax*nmax+12*nmax+16 )
382 work( 1 ) = maxwrk
383 END IF
384*
385 IF( lwork.LT.minwrk )
386 $ info = -24
387*
388 IF( info.NE.0 ) THEN
389 CALL xerbla( 'DDRGVX', -info )
390 RETURN
391 END IF
392*
393 n = 5
394 ulp = dlamch( 'P' )
395 ulpinv = one / ulp
396 thrsh2 = ten*thresh
397 nerrs = 0
398 nptknt = 0
399 ntestt = 0
400*
401 IF( nsize.EQ.0 )
402 $ GO TO 90
403*
404* Parameters used for generating test matrices.
405*
406 weight( 1 ) = tnth
407 weight( 2 ) = half
408 weight( 3 ) = one
409 weight( 4 ) = one / weight( 2 )
410 weight( 5 ) = one / weight( 1 )
411*
412 DO 80 iptype = 1, 2
413 DO 70 iwa = 1, 5
414 DO 60 iwb = 1, 5
415 DO 50 iwx = 1, 5
416 DO 40 iwy = 1, 5
417*
418* generated a test matrix pair
419*
420 CALL dlatm6( iptype, 5, a, lda, b, vr, lda, vl,
421 $ lda, weight( iwa ), weight( iwb ),
422 $ weight( iwx ), weight( iwy ), dtru,
423 $ diftru )
424*
425* Compute eigenvalues/eigenvectors of (A, B).
426* Compute eigenvalue/eigenvector condition numbers
427* using computed eigenvectors.
428*
429 CALL dlacpy( 'F', n, n, a, lda, ai, lda )
430 CALL dlacpy( 'F', n, n, b, lda, bi, lda )
431*
432 CALL dggevx( 'N', 'V', 'V', 'B', n, ai, lda, bi,
433 $ lda, alphar, alphai, beta, vl, lda,
434 $ vr, lda, ilo, ihi, lscale, rscale,
435 $ anorm, bnorm, s, dif, work, lwork,
436 $ iwork, bwork, linfo )
437 IF( linfo.NE.0 ) THEN
438 result( 1 ) = ulpinv
439 WRITE( nout, fmt = 9999 )'DGGEVX', linfo, n,
440 $ iptype
441 GO TO 30
442 END IF
443*
444* Compute the norm(A, B)
445*
446 CALL dlacpy( 'Full', n, n, ai, lda, work, n )
447 CALL dlacpy( 'Full', n, n, bi, lda, work( n*n+1 ),
448 $ n )
449 abnorm = dlange( 'Fro', n, 2*n, work, n, work )
450*
451* Tests (1) and (2)
452*
453 result( 1 ) = zero
454 CALL dget52( .true., n, a, lda, b, lda, vl, lda,
455 $ alphar, alphai, beta, work,
456 $ result( 1 ) )
457 IF( result( 2 ).GT.thresh ) THEN
458 WRITE( nout, fmt = 9998 )'Left', 'DGGEVX',
459 $ result( 2 ), n, iptype, iwa, iwb, iwx, iwy
460 END IF
461*
462 result( 2 ) = zero
463 CALL dget52( .false., n, a, lda, b, lda, vr, lda,
464 $ alphar, alphai, beta, work,
465 $ result( 2 ) )
466 IF( result( 3 ).GT.thresh ) THEN
467 WRITE( nout, fmt = 9998 )'Right', 'DGGEVX',
468 $ result( 3 ), n, iptype, iwa, iwb, iwx, iwy
469 END IF
470*
471* Test (3)
472*
473 result( 3 ) = zero
474 DO 10 i = 1, n
475 IF( s( i ).EQ.zero ) THEN
476 IF( dtru( i ).GT.abnorm*ulp )
477 $ result( 3 ) = ulpinv
478 ELSE IF( dtru( i ).EQ.zero ) THEN
479 IF( s( i ).GT.abnorm*ulp )
480 $ result( 3 ) = ulpinv
481 ELSE
482 work( i ) = max( abs( dtru( i ) / s( i ) ),
483 $ abs( s( i ) / dtru( i ) ) )
484 result( 3 ) = max( result( 3 ), work( i ) )
485 END IF
486 10 CONTINUE
487*
488* Test (4)
489*
490 result( 4 ) = zero
491 IF( dif( 1 ).EQ.zero ) THEN
492 IF( diftru( 1 ).GT.abnorm*ulp )
493 $ result( 4 ) = ulpinv
494 ELSE IF( diftru( 1 ).EQ.zero ) THEN
495 IF( dif( 1 ).GT.abnorm*ulp )
496 $ result( 4 ) = ulpinv
497 ELSE IF( dif( 5 ).EQ.zero ) THEN
498 IF( diftru( 5 ).GT.abnorm*ulp )
499 $ result( 4 ) = ulpinv
500 ELSE IF( diftru( 5 ).EQ.zero ) THEN
501 IF( dif( 5 ).GT.abnorm*ulp )
502 $ result( 4 ) = ulpinv
503 ELSE
504 ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
505 $ abs( dif( 1 ) / diftru( 1 ) ) )
506 ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
507 $ abs( dif( 5 ) / diftru( 5 ) ) )
508 result( 4 ) = max( ratio1, ratio2 )
509 END IF
510*
511 ntestt = ntestt + 4
512*
513* Print out tests which fail.
514*
515 DO 20 j = 1, 4
516 IF( ( result( j ).GE.thrsh2 .AND. j.GE.4 ) .OR.
517 $ ( result( j ).GE.thresh .AND. j.LE.3 ) )
518 $ THEN
519*
520* If this is the first test to fail,
521* print a header to the data file.
522*
523 IF( nerrs.EQ.0 ) THEN
524 WRITE( nout, fmt = 9997 )'DXV'
525*
526* Print out messages for built-in examples
527*
528* Matrix types
529*
530 WRITE( nout, fmt = 9995 )
531 WRITE( nout, fmt = 9994 )
532 WRITE( nout, fmt = 9993 )
533*
534* Tests performed
535*
536 WRITE( nout, fmt = 9992 )'''',
537 $ 'transpose', ''''
538*
539 END IF
540 nerrs = nerrs + 1
541 IF( result( j ).LT.10000.0d0 ) THEN
542 WRITE( nout, fmt = 9991 )iptype, iwa,
543 $ iwb, iwx, iwy, j, result( j )
544 ELSE
545 WRITE( nout, fmt = 9990 )iptype, iwa,
546 $ iwb, iwx, iwy, j, result( j )
547 END IF
548 END IF
549 20 CONTINUE
550*
551 30 CONTINUE
552*
553 40 CONTINUE
554 50 CONTINUE
555 60 CONTINUE
556 70 CONTINUE
557 80 CONTINUE
558*
559 GO TO 150
560*
561 90 CONTINUE
562*
563* Read in data from file to check accuracy of condition estimation
564* Read input data until N=0
565*
566 READ( nin, fmt = *, END = 150 )n
567 IF( n.EQ.0 )
568 $ GO TO 150
569 DO 100 i = 1, n
570 READ( nin, fmt = * )( a( i, j ), j = 1, n )
571 100 CONTINUE
572 DO 110 i = 1, n
573 READ( nin, fmt = * )( b( i, j ), j = 1, n )
574 110 CONTINUE
575 READ( nin, fmt = * )( dtru( i ), i = 1, n )
576 READ( nin, fmt = * )( diftru( i ), i = 1, n )
577*
578 nptknt = nptknt + 1
579*
580* Compute eigenvalues/eigenvectors of (A, B).
581* Compute eigenvalue/eigenvector condition numbers
582* using computed eigenvectors.
583*
584 CALL dlacpy( 'F', n, n, a, lda, ai, lda )
585 CALL dlacpy( 'F', n, n, b, lda, bi, lda )
586*
587 CALL dggevx( 'n', 'v', 'v', 'b', N, AI, LDA, BI, LDA, ALPHAR,
588 $ ALPHAI, BETA, VL, LDA, VR, LDA, ILO, IHI, LSCALE,
589 $ RSCALE, ANORM, BNORM, S, DIF, WORK, LWORK, IWORK,
590 $ BWORK, LINFO )
591*
592.NE. IF( LINFO0 ) THEN
593 RESULT( 1 ) = ULPINV
594 WRITE( NOUT, FMT = 9987 )'dggevx', LINFO, N, NPTKNT
595 GO TO 140
596 END IF
597*
598* Compute the norm(A, B)
599*
600 CALL DLACPY( 'full', N, N, AI, LDA, WORK, N )
601 CALL DLACPY( 'full', N, N, BI, LDA, WORK( N*N+1 ), N )
602 ABNORM = DLANGE( 'fro', N, 2*N, WORK, N, WORK )
603*
604* Tests (1) and (2)
605*
606 RESULT( 1 ) = ZERO
607 CALL DGET52( .TRUE., N, A, LDA, B, LDA, VL, LDA, ALPHAR, ALPHAI,
608 $ BETA, WORK, RESULT( 1 ) )
609.GT. IF( RESULT( 2 )THRESH ) THEN
610 WRITE( NOUT, FMT = 9986 )'left', 'dggevx', RESULT( 2 ), N,
611 $ NPTKNT
612 END IF
613*
614 RESULT( 2 ) = ZERO
615 CALL DGET52( .FALSE., N, A, LDA, B, LDA, VR, LDA, ALPHAR, ALPHAI,
616 $ BETA, WORK, RESULT( 2 ) )
617.GT. IF( RESULT( 3 )THRESH ) THEN
618 WRITE( NOUT, FMT = 9986 )'right', 'dggevx', RESULT( 3 ), N,
619 $ NPTKNT
620 END IF
621*
622* Test (3)
623*
624 RESULT( 3 ) = ZERO
625 DO 120 I = 1, N
626.EQ. IF( S( I )ZERO ) THEN
627.GT. IF( DTRU( I )ABNORM*ULP )
628 $ RESULT( 3 ) = ULPINV
629.EQ. ELSE IF( DTRU( I )ZERO ) THEN
630.GT. IF( S( I )ABNORM*ULP )
631 $ RESULT( 3 ) = ULPINV
632 ELSE
633 WORK( I ) = MAX( ABS( DTRU( I ) / S( I ) ),
634 $ ABS( S( I ) / DTRU( I ) ) )
635 RESULT( 3 ) = MAX( RESULT( 3 ), WORK( I ) )
636 END IF
637 120 CONTINUE
638*
639* Test (4)
640*
641 RESULT( 4 ) = ZERO
642.EQ. IF( DIF( 1 )ZERO ) THEN
643.GT. IF( DIFTRU( 1 )ABNORM*ULP )
644 $ RESULT( 4 ) = ULPINV
645.EQ. ELSE IF( DIFTRU( 1 )ZERO ) THEN
646.GT. IF( DIF( 1 )ABNORM*ULP )
647 $ RESULT( 4 ) = ULPINV
648.EQ. ELSE IF( DIF( 5 )ZERO ) THEN
649.GT. IF( DIFTRU( 5 )ABNORM*ULP )
650 $ RESULT( 4 ) = ULPINV
651.EQ. ELSE IF( DIFTRU( 5 )ZERO ) THEN
652.GT. IF( DIF( 5 )ABNORM*ULP )
653 $ RESULT( 4 ) = ULPINV
654 ELSE
655 RATIO1 = MAX( ABS( DIFTRU( 1 ) / DIF( 1 ) ),
656 $ ABS( DIF( 1 ) / DIFTRU( 1 ) ) )
657 RATIO2 = MAX( ABS( DIFTRU( 5 ) / DIF( 5 ) ),
658 $ ABS( DIF( 5 ) / DIFTRU( 5 ) ) )
659 RESULT( 4 ) = MAX( RATIO1, RATIO2 )
660 END IF
661*
662 NTESTT = NTESTT + 4
663*
664* Print out tests which fail.
665*
666 DO 130 J = 1, 4
667.GE. IF( RESULT( J )THRSH2 ) THEN
668*
669* If this is the first test to fail,
670* print a header to the data file.
671*
672.EQ. IF( NERRS0 ) THEN
673 WRITE( NOUT, FMT = 9997 )'dxv'
674*
675* Print out messages for built-in examples
676*
677* Matrix types
678*
679 WRITE( NOUT, FMT = 9996 )
680*
681* Tests performed
682*
683 WRITE( NOUT, FMT = 9992 )'''', 'transpose', ''''
684*
685 END IF
686 NERRS = NERRS + 1
687.LT. IF( RESULT( J )10000.0D0 ) THEN
688 WRITE( NOUT, FMT = 9989 )NPTKNT, N, J, RESULT( J )
689 ELSE
690 WRITE( NOUT, FMT = 9988 )NPTKNT, N, J, RESULT( J )
691 END IF
692 END IF
693 130 CONTINUE
694*
695 140 CONTINUE
696*
697 GO TO 90
698 150 CONTINUE
699*
700* Summary
701*
702 CALL ALASVM( 'dxv', NOUT, NERRS, NTESTT, 0 )
703*
704 WORK( 1 ) = MAXWRK
705*
706 RETURN
707*
708 9999 FORMAT( ' ddrgvx: ', A, ' returned info=', I6, '.', / 9X, 'n=',
709 $ I6, ', jtype=', I6, ')' )
710*
711 9998 FORMAT( ' ddrgvx: ', A, ' eigenvectors from ', A, ' incorrectly ',
712 $ 'normalized.', / ' bits of error=', 0P, G10.3, ',', 9X,
713 $ 'n=', I6, ', jtype=', I6, ', iwa=', I5, ', iwb=', i5,
714 $ ', IWX=', i5, ', IWY=', i5 )
715*
716 9997 FORMAT( / 1x, a3, ' -- Real Expert Eigenvalue/vector',
717 $ ' problem driver' )
718*
719 9996 FORMAT( ' input example' )
720*
721 9995 FORMAT( ' matrix types: ', / )
722*
723 9994 FORMAT( ' TYPE 1: da is diagonal, db is identity, ',
724 $ / ' a = y^(-h) da x^(-1), b = y^(-h) db x^(-1) ',
725 $ / ' yh and x are left and right eigenvectors. ', / )
726*
727 9993 FORMAT( ' TYPE 2: da is quasi-diagonal, db is identity, ',
728 $ / ' a = y^(-h) da x^(-1), b = y^(-h) db x^(-1) ',
729 $ / ' yh and x are left and right eigenvectors. ', / )
730*
731 9992 FORMAT( / ' tests performed: ', / 4X,
732 $ ' a is alpha, b is beta, l is a left eigenvector, ', / 4X,
733 $ ' r is a right eigenvector and ', A, ' means ', A, '.',
734 $ / ' 1 = max | ( b a - a b )', A, ' l | / const.',
735 $ / ' 2 = max | ( b a - a b ) r | / const.',
736 $ / ' 3 = max( sest/stru, stru/sest ) ',
737 $ ' over all eigenvalues', /
738 $ ' 4 = max( difest/diftru, diftru/difest ) ',
739 $ ' over the 1st and 5th eigenvectors', / )
740*
741 9991 FORMAT( ' type=', I2, ',', ' iwa=', I2, ', iwb=', I2, ', iwx=',
742 $ I2, ', iwy=', I2, ', result ', I2, ' is', 0P, F8.2 )
743 9990 FORMAT( ' type=', I2, ',', ' iwa=', I2, ', iwb=', I2, ', iwx=',
744 $ I2, ', iwy=', I2, ', result ', I2, ' is', 1P, D10.3 )
745 9989 FORMAT( ' input example #', I2, ', matrix order=', I4, ',',
746 $ ' result ', i2, ' is', 0p, f8.2 )
747 9988 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
748 $ ' result ', i2, ' is', 1p, d10.3 )
749 9987 FORMAT( ' DDRGVX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
750 $ i6, ', Input example #', i2, ')' )
751*
752 9986 FORMAT( ' DDRGVX: ', a, ' Eigenvectors from ', a, ' incorrectly ',
753 $ 'normalized.', / ' Bits of error=', 0p, g10.3, ',', 9x,
754 $ 'N=', i6, ', Input Example #', i2, ')' )
755*
756*
757* End of DDRGVX
758*
759 END
the
end diagonal values have been computed in the(sparse) matrix id.SOL
alpha
#define alpha
Definition eval.h:35
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
alasvm
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
Definition alasvm.f:73
dggevx
subroutine dggevx(balanc, jobvl, jobvr, sense, n, a, lda, b, ldb, alphar, alphai, beta, vl, ldvl, vr, ldvr, ilo, ihi, lscale, rscale, abnrm, bbnrm, rconde, rcondv, work, lwork, iwork, bwork, info)
DGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition dggevx.f:391
ddrgvx
subroutine ddrgvx(nsize, thresh, nin, nout, a, lda, b, ai, bi, alphar, alphai, beta, vl, vr, ilo, ihi, lscale, rscale, s, dtru, dif, diftru, work, lwork, iwork, liwork, result, bwork, info)
DDRGVX
Definition ddrgvx.f:300
dget52
subroutine dget52(left, n, a, lda, b, ldb, e, lde, alphar, alphai, beta, work, result)
DGET52
Definition dget52.f:199
dlatm6
subroutine dlatm6(type, n, a, lda, b, x, ldx, y, ldy, alpha, beta, wx, wy, s, dif)
DLATM6
Definition dlatm6.f:176
max
#define max(a, b)
Definition macros.h:21