LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
ddrgsx.f
Go to the documentation of this file.
1*> \brief \b DDRGSX
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 DDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI,
12* BI, Z, Q, ALPHAR, ALPHAI, BETA, C, LDC, S,
13* WORK, LWORK, IWORK, LIWORK, BWORK, INFO )
14*
15* .. Scalar Arguments ..
16* INTEGER INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
17* $ NOUT, NSIZE
18* DOUBLE PRECISION THRESH
19* ..
20* .. Array Arguments ..
21* LOGICAL BWORK( * )
22* INTEGER IWORK( * )
23* DOUBLE PRECISION A( LDA, * ), AI( LDA, * ), ALPHAI( * ),
24* $ ALPHAR( * ), B( LDA, * ), BETA( * ),
25* $ BI( LDA, * ), C( LDC, * ), Q( LDA, * ), S( * ),
26* $ WORK( * ), Z( LDA, * )
27* ..
28*
29*
30*> \par Purpose:
31* =============
32*>
33*> \verbatim
34*>
35*> DDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
36*> problem expert driver DGGESX.
37*>
38*> DGGESX factors A and B as Q S Z' and Q T Z', where ' means
39*> transpose, T is upper triangular, S is in generalized Schur form
40*> (block upper triangular, with 1x1 and 2x2 blocks on the diagonal,
41*> the 2x2 blocks corresponding to complex conjugate pairs of
42*> generalized eigenvalues), and Q and Z are orthogonal. It also
43*> computes the generalized eigenvalues (alpha(1),beta(1)), ...,
44*> (alpha(n),beta(n)). Thus, w(j) = alpha(j)/beta(j) is a root of the
45*> characteristic equation
46*>
47*> det( A - w(j) B ) = 0
48*>
49*> Optionally it also reorders the eigenvalues so that a selected
50*> cluster of eigenvalues appears in the leading diagonal block of the
51*> Schur forms; computes a reciprocal condition number for the average
52*> of the selected eigenvalues; and computes a reciprocal condition
53*> number for the right and left deflating subspaces corresponding to
54*> the selected eigenvalues.
55*>
56*> When DDRGSX is called with NSIZE > 0, five (5) types of built-in
57*> matrix pairs are used to test the routine DGGESX.
58*>
59*> When DDRGSX is called with NSIZE = 0, it reads in test matrix data
60*> to test DGGESX.
61*>
62*> For each matrix pair, the following tests will be performed and
63*> compared with the threshold THRESH except for the tests (7) and (9):
64*>
65*> (1) | A - Q S Z' | / ( |A| n ulp )
66*>
67*> (2) | B - Q T Z' | / ( |B| n ulp )
68*>
69*> (3) | I - QQ' | / ( n ulp )
70*>
71*> (4) | I - ZZ' | / ( n ulp )
72*>
73*> (5) if A is in Schur form (i.e. quasi-triangular form)
74*>
75*> (6) maximum over j of D(j) where:
76*>
77*> if alpha(j) is real:
78*> |alpha(j) - S(j,j)| |beta(j) - T(j,j)|
79*> D(j) = ------------------------ + -----------------------
80*> max(|alpha(j)|,|S(j,j)|) max(|beta(j)|,|T(j,j)|)
81*>
82*> if alpha(j) is complex:
83*> | det( s S - w T ) |
84*> D(j) = ---------------------------------------------------
85*> ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )
86*>
87*> and S and T are here the 2 x 2 diagonal blocks of S and T
88*> corresponding to the j-th and j+1-th eigenvalues.
89*>
90*> (7) if sorting worked and SDIM is the number of eigenvalues
91*> which were selected.
92*>
93*> (8) the estimated value DIF does not differ from the true values of
94*> Difu and Difl more than a factor 10*THRESH. If the estimate DIF
95*> equals zero the corresponding true values of Difu and Difl
96*> should be less than EPS*norm(A, B). If the true value of Difu
97*> and Difl equal zero, the estimate DIF should be less than
98*> EPS*norm(A, B).
99*>
100*> (9) If INFO = N+3 is returned by DGGESX, the reordering "failed"
101*> and we check that DIF = PL = PR = 0 and that the true value of
102*> Difu and Difl is < EPS*norm(A, B). We count the events when
103*> INFO=N+3.
104*>
105*> For read-in test matrices, the above tests are run except that the
106*> exact value for DIF (and PL) is input data. Additionally, there is
107*> one more test run for read-in test matrices:
108*>
109*> (10) the estimated value PL does not differ from the true value of
110*> PLTRU more than a factor THRESH. If the estimate PL equals
111*> zero the corresponding true value of PLTRU should be less than
112*> EPS*norm(A, B). If the true value of PLTRU equal zero, the
113*> estimate PL should be less than EPS*norm(A, B).
114*>
115*> Note that for the built-in tests, a total of 10*NSIZE*(NSIZE-1)
116*> matrix pairs are generated and tested. NSIZE should be kept small.
117*>
118*> SVD (routine DGESVD) is used for computing the true value of DIF_u
119*> and DIF_l when testing the built-in test problems.
120*>
121*> Built-in Test Matrices
122*> ======================
123*>
124*> All built-in test matrices are the 2 by 2 block of triangular
125*> matrices
126*>
127*> A = [ A11 A12 ] and B = [ B11 B12 ]
128*> [ A22 ] [ B22 ]
129*>
130*> where for different type of A11 and A22 are given as the following.
131*> A12 and B12 are chosen so that the generalized Sylvester equation
132*>
133*> A11*R - L*A22 = -A12
134*> B11*R - L*B22 = -B12
135*>
136*> have prescribed solution R and L.
137*>
138*> Type 1: A11 = J_m(1,-1) and A_22 = J_k(1-a,1).
139*> B11 = I_m, B22 = I_k
140*> where J_k(a,b) is the k-by-k Jordan block with ``a'' on
141*> diagonal and ``b'' on superdiagonal.
142*>
143*> Type 2: A11 = (a_ij) = ( 2(.5-sin(i)) ) and
144*> B11 = (b_ij) = ( 2(.5-sin(ij)) ) for i=1,...,m, j=i,...,m
145*> A22 = (a_ij) = ( 2(.5-sin(i+j)) ) and
146*> B22 = (b_ij) = ( 2(.5-sin(ij)) ) for i=m+1,...,k, j=i,...,k
147*>
148*> Type 3: A11, A22 and B11, B22 are chosen as for Type 2, but each
149*> second diagonal block in A_11 and each third diagonal block
150*> in A_22 are made as 2 by 2 blocks.
151*>
152*> Type 4: A11 = ( 20(.5 - sin(ij)) ) and B22 = ( 2(.5 - sin(i+j)) )
153*> for i=1,...,m, j=1,...,m and
154*> A22 = ( 20(.5 - sin(i+j)) ) and B22 = ( 2(.5 - sin(ij)) )
155*> for i=m+1,...,k, j=m+1,...,k
156*>
157*> Type 5: (A,B) and have potentially close or common eigenvalues and
158*> very large departure from block diagonality A_11 is chosen
159*> as the m x m leading submatrix of A_1:
160*> | 1 b |
161*> | -b 1 |
162*> | 1+d b |
163*> | -b 1+d |
164*> A_1 = | d 1 |
165*> | -1 d |
166*> | -d 1 |
167*> | -1 -d |
168*> | 1 |
169*> and A_22 is chosen as the k x k leading submatrix of A_2:
170*> | -1 b |
171*> | -b -1 |
172*> | 1-d b |
173*> | -b 1-d |
174*> A_2 = | d 1+b |
175*> | -1-b d |
176*> | -d 1+b |
177*> | -1+b -d |
178*> | 1-d |
179*> and matrix B are chosen as identity matrices (see DLATM5).
180*>
181*> \endverbatim
182*
183* Arguments:
184* ==========
185*
186*> \param[in] NSIZE
187*> \verbatim
188*> NSIZE is INTEGER
189*> The maximum size of the matrices to use. NSIZE >= 0.
190*> If NSIZE = 0, no built-in tests matrices are used, but
191*> read-in test matrices are used to test DGGESX.
192*> \endverbatim
193*>
194*> \param[in] NCMAX
195*> \verbatim
196*> NCMAX is INTEGER
197*> Maximum allowable NMAX for generating Kroneker matrix
198*> in call to DLAKF2
199*> \endverbatim
200*>
201*> \param[in] THRESH
202*> \verbatim
203*> THRESH is DOUBLE PRECISION
204*> A test will count as "failed" if the "error", computed as
205*> described above, exceeds THRESH. Note that the error
206*> is scaled to be O(1), so THRESH should be a reasonably
207*> small multiple of 1, e.g., 10 or 100. In particular,
208*> it should not depend on the precision (single vs. double)
209*> or the size of the matrix. THRESH >= 0.
210*> \endverbatim
211*>
212*> \param[in] NIN
213*> \verbatim
214*> NIN is INTEGER
215*> The FORTRAN unit number for reading in the data file of
216*> problems to solve.
217*> \endverbatim
218*>
219*> \param[in] NOUT
220*> \verbatim
221*> NOUT is INTEGER
222*> The FORTRAN unit number for printing out error messages
223*> (e.g., if a routine returns IINFO not equal to 0.)
224*> \endverbatim
225*>
226*> \param[out] A
227*> \verbatim
228*> A is DOUBLE PRECISION array, dimension (LDA, NSIZE)
229*> Used to store the matrix whose eigenvalues are to be
230*> computed. On exit, A contains the last matrix actually used.
231*> \endverbatim
232*>
233*> \param[in] LDA
234*> \verbatim
235*> LDA is INTEGER
236*> The leading dimension of A, B, AI, BI, Z and Q,
237*> LDA >= max( 1, NSIZE ). For the read-in test,
238*> LDA >= max( 1, N ), N is the size of the test matrices.
239*> \endverbatim
240*>
241*> \param[out] B
242*> \verbatim
243*> B is DOUBLE PRECISION array, dimension (LDA, NSIZE)
244*> Used to store the matrix whose eigenvalues are to be
245*> computed. On exit, B contains the last matrix actually used.
246*> \endverbatim
247*>
248*> \param[out] AI
249*> \verbatim
250*> AI is DOUBLE PRECISION array, dimension (LDA, NSIZE)
251*> Copy of A, modified by DGGESX.
252*> \endverbatim
253*>
254*> \param[out] BI
255*> \verbatim
256*> BI is DOUBLE PRECISION array, dimension (LDA, NSIZE)
257*> Copy of B, modified by DGGESX.
258*> \endverbatim
259*>
260*> \param[out] Z
261*> \verbatim
262*> Z is DOUBLE PRECISION array, dimension (LDA, NSIZE)
263*> Z holds the left Schur vectors computed by DGGESX.
264*> \endverbatim
265*>
266*> \param[out] Q
267*> \verbatim
268*> Q is DOUBLE PRECISION array, dimension (LDA, NSIZE)
269*> Q holds the right Schur vectors computed by DGGESX.
270*> \endverbatim
271*>
272*> \param[out] ALPHAR
273*> \verbatim
274*> ALPHAR is DOUBLE PRECISION array, dimension (NSIZE)
275*> \endverbatim
276*>
277*> \param[out] ALPHAI
278*> \verbatim
279*> ALPHAI is DOUBLE PRECISION array, dimension (NSIZE)
280*> \endverbatim
281*>
282*> \param[out] BETA
283*> \verbatim
284*> BETA is DOUBLE PRECISION array, dimension (NSIZE)
285*>
286*> On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues.
287*> \endverbatim
288*>
289*> \param[out] C
290*> \verbatim
291*> C is DOUBLE PRECISION array, dimension (LDC, LDC)
292*> Store the matrix generated by subroutine DLAKF2, this is the
293*> matrix formed by Kronecker products used for estimating
294*> DIF.
295*> \endverbatim
296*>
297*> \param[in] LDC
298*> \verbatim
299*> LDC is INTEGER
300*> The leading dimension of C. LDC >= max(1, LDA*LDA/2 ).
301*> \endverbatim
302*>
303*> \param[out] S
304*> \verbatim
305*> S is DOUBLE PRECISION array, dimension (LDC)
306*> Singular values of C
307*> \endverbatim
308*>
309*> \param[out] WORK
310*> \verbatim
311*> WORK is DOUBLE PRECISION array, dimension (LWORK)
312*> \endverbatim
313*>
314*> \param[in] LWORK
315*> \verbatim
316*> LWORK is INTEGER
317*> The dimension of the array WORK.
318*> LWORK >= MAX( 5*NSIZE*NSIZE/2 - 2, 10*(NSIZE+1) )
319*> \endverbatim
320*>
321*> \param[out] IWORK
322*> \verbatim
323*> IWORK is INTEGER array, dimension (LIWORK)
324*> \endverbatim
325*>
326*> \param[in] LIWORK
327*> \verbatim
328*> LIWORK is INTEGER
329*> The dimension of the array IWORK. LIWORK >= NSIZE + 6.
330*> \endverbatim
331*>
332*> \param[out] BWORK
333*> \verbatim
334*> BWORK is LOGICAL array, dimension (LDA)
335*> \endverbatim
336*>
337*> \param[out] INFO
338*> \verbatim
339*> INFO is INTEGER
340*> = 0: successful exit
341*> < 0: if INFO = -i, the i-th argument had an illegal value.
342*> > 0: A routine returned an error code.
343*> \endverbatim
344*
345* Authors:
346* ========
347*
348*> \author Univ. of Tennessee
349*> \author Univ. of California Berkeley
350*> \author Univ. of Colorado Denver
351*> \author NAG Ltd.
352*
353*> \ingroup double_eig
354*
355* =====================================================================
356 SUBROUTINE ddrgsx( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI,
357 $ BI, Z, Q, ALPHAR, ALPHAI, BETA, C, LDC, S,
358 $ WORK, LWORK, IWORK, LIWORK, BWORK, INFO )
359*
360* -- LAPACK test routine --
361* -- LAPACK is a software package provided by Univ. of Tennessee, --
362* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
363*
364* .. Scalar Arguments ..
365 INTEGER INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
366 $ NOUT, NSIZE
367 DOUBLE PRECISION THRESH
368* ..
369* .. Array Arguments ..
370 LOGICAL BWORK( * )
371 INTEGER IWORK( * )
372 DOUBLE PRECISION A( LDA, * ), AI( LDA, * ), ALPHAI( * ),
373 $ alphar( * ), b( lda, * ), beta( * ),
374 $ bi( lda, * ), c( ldc, * ), q( lda, * ), s( * ),
375 $ work( * ), z( lda, * )
376* ..
377*
378* =====================================================================
379*
380* .. Parameters ..
381 DOUBLE PRECISION ZERO, ONE, TEN
382 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0, ten = 1.0d+1 )
383* ..
384* .. Local Scalars ..
385 LOGICAL ILABAD
386 CHARACTER SENSE
387 INTEGER BDSPAC, I, I1, IFUNC, IINFO, J, LINFO, MAXWRK,
388 $ minwrk, mm, mn2, nerrs, nptknt, ntest, ntestt,
389 $ prtype, qba, qbb
390 DOUBLE PRECISION ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
391 $ TEMP2, THRSH2, ULP, ULPINV, WEIGHT
392* ..
393* .. Local Arrays ..
394 DOUBLE PRECISION DIFEST( 2 ), PL( 2 ), RESULT( 10 )
395* ..
396* .. External Functions ..
397 LOGICAL DLCTSX
398 INTEGER ILAENV
399 DOUBLE PRECISION DLAMCH, DLANGE
400 EXTERNAL dlctsx, ilaenv, dlamch, dlange
401* ..
402* .. External Subroutines ..
403 EXTERNAL alasvm, dgesvd, dget51, dget53, dggesx, dlabad,
405* ..
406* .. Intrinsic Functions ..
407 INTRINSIC abs, max, sqrt
408* ..
409* .. Scalars in Common ..
410 LOGICAL FS
411 INTEGER K, M, MPLUSN, N
412* ..
413* .. Common blocks ..
414 COMMON / mn / m, n, mplusn, k, fs
415* ..
416* .. Executable Statements ..
417*
418* Check for errors
419*
420 IF( nsize.LT.0 ) THEN
421 info = -1
422 ELSE IF( thresh.LT.zero ) THEN
423 info = -2
424 ELSE IF( nin.LE.0 ) THEN
425 info = -3
426 ELSE IF( nout.LE.0 ) THEN
427 info = -4
428 ELSE IF( lda.LT.1 .OR. lda.LT.nsize ) THEN
429 info = -6
430 ELSE IF( ldc.LT.1 .OR. ldc.LT.nsize*nsize / 2 ) THEN
431 info = -17
432 ELSE IF( liwork.LT.nsize+6 ) THEN
433 info = -21
434 END IF
435*
436* Compute workspace
437* (Note: Comments in the code beginning "Workspace:" describe the
438* minimal amount of workspace needed at that point in the code,
439* as well as the preferred amount for good performance.
440* NB refers to the optimal block size for the immediately
441* following subroutine, as returned by ILAENV.)
442*
443 minwrk = 1
444 IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
445 minwrk = max( 10*( nsize+1 ), 5*nsize*nsize / 2 )
446*
447* workspace for sggesx
448*
449 maxwrk = 9*( nsize+1 ) + nsize*
450 $ ilaenv( 1, 'DGEQRF', ' ', nsize, 1, nsize, 0 )
451 maxwrk = max( maxwrk, 9*( nsize+1 )+nsize*
452 $ ilaenv( 1, 'DORGQR', ' ', nsize, 1, nsize, -1 ) )
453*
454* workspace for dgesvd
455*
456 bdspac = 5*nsize*nsize / 2
457 maxwrk = max( maxwrk, 3*nsize*nsize / 2+nsize*nsize*
458 $ ilaenv( 1, 'DGEBRD', ' ', nsize*nsize / 2,
459 $ nsize*nsize / 2, -1, -1 ) )
460 maxwrk = max( maxwrk, bdspac )
461*
462 maxwrk = max( maxwrk, minwrk )
463*
464 work( 1 ) = maxwrk
465 END IF
466*
467 IF( lwork.LT.minwrk )
468 $ info = -19
469*
470 IF( info.NE.0 ) THEN
471 CALL xerbla( 'DDRGSX', -info )
472 RETURN
473 END IF
474*
475* Important constants
476*
477 ulp = dlamch( 'P' )
478 ulpinv = one / ulp
479 smlnum = dlamch( 'S' ) / ulp
480 bignum = one / smlnum
481 CALL dlabad( smlnum, bignum )
482 thrsh2 = ten*thresh
483 ntestt = 0
484 nerrs = 0
485*
486* Go to the tests for read-in matrix pairs
487*
488 ifunc = 0
489 IF( nsize.EQ.0 )
490 $ GO TO 70
491*
492* Test the built-in matrix pairs.
493* Loop over different functions (IFUNC) of DGGESX, types (PRTYPE)
494* of test matrices, different size (M+N)
495*
496 prtype = 0
497 qba = 3
498 qbb = 4
499 weight = sqrt( ulp )
500*
501 DO 60 ifunc = 0, 3
502 DO 50 prtype = 1, 5
503 DO 40 m = 1, nsize - 1
504 DO 30 n = 1, nsize - m
505*
506 weight = one / weight
507 mplusn = m + n
508*
509* Generate test matrices
510*
511 fs = .true.
512 k = 0
513*
514 CALL dlaset( 'Full', mplusn, mplusn, zero, zero, ai,
515 $ lda )
516 CALL dlaset( 'Full', mplusn, mplusn, zero, zero, bi,
517 $ lda )
518*
519 CALL dlatm5( prtype, m, n, ai, lda, ai( m+1, m+1 ),
520 $ lda, ai( 1, m+1 ), lda, bi, lda,
521 $ bi( m+1, m+1 ), lda, bi( 1, m+1 ), lda,
522 $ q, lda, z, lda, weight, qba, qbb )
523*
524* Compute the Schur factorization and swapping the
525* m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
526* Swapping is accomplished via the function DLCTSX
527* which is supplied below.
528*
529 IF( ifunc.EQ.0 ) THEN
530 sense = 'N'
531 ELSE IF( ifunc.EQ.1 ) THEN
532 sense = 'E'
533 ELSE IF( ifunc.EQ.2 ) THEN
534 sense = 'V'
535 ELSE IF( ifunc.EQ.3 ) THEN
536 sense = 'B'
537 END IF
538*
539 CALL dlacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
540 CALL dlacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
541*
542 CALL dggesx( 'V', 'V', 'S', dlctsx, sense, mplusn, ai,
543 $ lda, bi, lda, mm, alphar, alphai, beta,
544 $ q, lda, z, lda, pl, difest, work, lwork,
545 $ iwork, liwork, bwork, linfo )
546*
547 IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
548 result( 1 ) = ulpinv
549 WRITE( nout, fmt = 9999 )'DGGESX', linfo, mplusn,
550 $ prtype
551 info = linfo
552 GO TO 30
553 END IF
554*
555* Compute the norm(A, B)
556*
557 CALL dlacpy( 'Full', mplusn, mplusn, ai, lda, work,
558 $ mplusn )
559 CALL dlacpy( 'Full', mplusn, mplusn, bi, lda,
560 $ work( mplusn*mplusn+1 ), mplusn )
561 abnrm = dlange( 'Fro', mplusn, 2*mplusn, work, mplusn,
562 $ work )
563*
564* Do tests (1) to (4)
565*
566 CALL dget51( 1, mplusn, a, lda, ai, lda, q, lda, z,
567 $ lda, work, result( 1 ) )
568 CALL dget51( 1, mplusn, b, lda, bi, lda, q, lda, z,
569 $ lda, work, result( 2 ) )
570 CALL dget51( 3, mplusn, b, lda, bi, lda, q, lda, q,
571 $ lda, work, result( 3 ) )
572 CALL dget51( 3, mplusn, b, lda, bi, lda, z, lda, z,
573 $ lda, work, result( 4 ) )
574 ntest = 4
575*
576* Do tests (5) and (6): check Schur form of A and
577* compare eigenvalues with diagonals.
578*
579 temp1 = zero
580 result( 5 ) = zero
581 result( 6 ) = zero
582*
583 DO 10 j = 1, mplusn
584 ilabad = .false.
585 IF( alphai( j ).EQ.zero ) THEN
586 temp2 = ( abs( alphar( j )-ai( j, j ) ) /
587 $ max( smlnum, abs( alphar( j ) ),
588 $ abs( ai( j, j ) ) )+
589 $ abs( beta( j )-bi( j, j ) ) /
590 $ max( smlnum, abs( beta( j ) ),
591 $ abs( bi( j, j ) ) ) ) / ulp
592 IF( j.LT.mplusn ) THEN
593 IF( ai( j+1, j ).NE.zero ) THEN
594 ilabad = .true.
595 result( 5 ) = ulpinv
596 END IF
597 END IF
598 IF( j.GT.1 ) THEN
599 IF( ai( j, j-1 ).NE.zero ) THEN
600 ilabad = .true.
601 result( 5 ) = ulpinv
602 END IF
603 END IF
604 ELSE
605 IF( alphai( j ).GT.zero ) THEN
606 i1 = j
607 ELSE
608 i1 = j - 1
609 END IF
610 IF( i1.LE.0 .OR. i1.GE.mplusn ) THEN
611 ilabad = .true.
612 ELSE IF( i1.LT.mplusn-1 ) THEN
613 IF( ai( i1+2, i1+1 ).NE.zero ) THEN
614 ilabad = .true.
615 result( 5 ) = ulpinv
616 END IF
617 ELSE IF( i1.GT.1 ) THEN
618 IF( ai( i1, i1-1 ).NE.zero ) THEN
619 ilabad = .true.
620 result( 5 ) = ulpinv
621 END IF
622 END IF
623 IF( .NOT.ilabad ) THEN
624 CALL dget53( ai( i1, i1 ), lda, bi( i1, i1 ),
625 $ lda, beta( j ), alphar( j ),
626 $ alphai( j ), temp2, iinfo )
627 IF( iinfo.GE.3 ) THEN
628 WRITE( nout, fmt = 9997 )iinfo, j,
629 $ mplusn, prtype
630 info = abs( iinfo )
631 END IF
632 ELSE
633 temp2 = ulpinv
634 END IF
635 END IF
636 temp1 = max( temp1, temp2 )
637 IF( ilabad ) THEN
638 WRITE( nout, fmt = 9996 )j, mplusn, prtype
639 END IF
640 10 CONTINUE
641 result( 6 ) = temp1
642 ntest = ntest + 2
643*
644* Test (7) (if sorting worked)
645*
646 result( 7 ) = zero
647 IF( linfo.EQ.mplusn+3 ) THEN
648 result( 7 ) = ulpinv
649 ELSE IF( mm.NE.n ) THEN
650 result( 7 ) = ulpinv
651 END IF
652 ntest = ntest + 1
653*
654* Test (8): compare the estimated value DIF and its
655* value. first, compute the exact DIF.
656*
657 result( 8 ) = zero
658 mn2 = mm*( mplusn-mm )*2
659 IF( ifunc.GE.2 .AND. mn2.LE.ncmax*ncmax ) THEN
660*
661* Note: for either following two causes, there are
662* almost same number of test cases fail the test.
663*
664 CALL dlakf2( mm, mplusn-mm, ai, lda,
665 $ ai( mm+1, mm+1 ), bi,
666 $ bi( mm+1, mm+1 ), c, ldc )
667*
668 CALL dgesvd( 'N', 'N', mn2, mn2, c, ldc, s, work,
669 $ 1, work( 2 ), 1, work( 3 ), lwork-2,
670 $ info )
671 diftru = s( mn2 )
672*
673 IF( difest( 2 ).EQ.zero ) THEN
674 IF( diftru.GT.abnrm*ulp )
675 $ result( 8 ) = ulpinv
676 ELSE IF( diftru.EQ.zero ) THEN
677 IF( difest( 2 ).GT.abnrm*ulp )
678 $ result( 8 ) = ulpinv
679 ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
680 $ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
681 result( 8 ) = max( diftru / difest( 2 ),
682 $ difest( 2 ) / diftru )
683 END IF
684 ntest = ntest + 1
685 END IF
686*
687* Test (9)
688*
689 result( 9 ) = zero
690 IF( linfo.EQ.( mplusn+2 ) ) THEN
691 IF( diftru.GT.abnrm*ulp )
692 $ result( 9 ) = ulpinv
693 IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
694 $ result( 9 ) = ulpinv
695 IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
696 $ result( 9 ) = ulpinv
697 ntest = ntest + 1
698 END IF
699*
700 ntestt = ntestt + ntest
701*
702* Print out tests which fail.
703*
704 DO 20 j = 1, 9
705 IF( result( j ).GE.thresh ) THEN
706*
707* If this is the first test to fail,
708* print a header to the data file.
709*
710 IF( nerrs.EQ.0 ) THEN
711 WRITE( nout, fmt = 9995 )'DGX'
712*
713* Matrix types
714*
715 WRITE( nout, fmt = 9993 )
716*
717* Tests performed
718*
719 WRITE( nout, fmt = 9992 )'orthogonal', '''',
720 $ 'transpose', ( '''', i = 1, 4 )
721*
722 END IF
723 nerrs = nerrs + 1
724 IF( result( j ).LT.10000.0d0 ) THEN
725 WRITE( nout, fmt = 9991 )mplusn, prtype,
726 $ weight, m, j, result( j )
727 ELSE
728 WRITE( nout, fmt = 9990 )mplusn, prtype,
729 $ weight, m, j, result( j )
730 END IF
731 END IF
732 20 CONTINUE
733*
734 30 CONTINUE
735 40 CONTINUE
736 50 CONTINUE
737 60 CONTINUE
738*
739 GO TO 150
740*
741 70 CONTINUE
742*
743* Read in data from file to check accuracy of condition estimation
744* Read input data until N=0
745*
746 nptknt = 0
747*
748 80 CONTINUE
749 READ( nin, fmt = *, END = 140 )mplusn
750 IF( mplusn.EQ.0 )
751 $ GO TO 140
752 READ( nin, fmt = *, END = 140 )n
753 DO 90 i = 1, mplusn
754 READ( nin, fmt = * )( ai( i, j ), j = 1, mplusn )
755 90 CONTINUE
756 DO 100 i = 1, mplusn
757 READ( nin, fmt = * )( bi( i, j ), j = 1, mplusn )
758 100 CONTINUE
759 READ( nin, fmt = * )pltru, diftru
760*
761 nptknt = nptknt + 1
762 fs = .true.
763 k = 0
764 m = mplusn - n
765*
766 CALL dlacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
767 CALL dlacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
768*
769* Compute the Schur factorization while swapping the
770* m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
771*
772 CALL dggesx( 'V', 'V', 'S', dlctsx, 'B', mplusn, ai, lda, bi, lda,
773 $ mm, alphar, alphai, beta, q, lda, z, lda, pl, difest,
774 $ work, lwork, iwork, liwork, bwork, linfo )
775*
776 IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
777 result( 1 ) = ulpinv
778 WRITE( nout, fmt = 9998 )'DGGESX', linfo, mplusn, nptknt
779 GO TO 130
780 END IF
781*
782* Compute the norm(A, B)
783* (should this be norm of (A,B) or (AI,BI)?)
784*
785 CALL dlacpy( 'Full', mplusn, mplusn, ai, lda, work, mplusn )
786 CALL dlacpy( 'Full', mplusn, mplusn, bi, lda,
787 $ work( mplusn*mplusn+1 ), mplusn )
788 abnrm = dlange( 'Fro', mplusn, 2*mplusn, work, mplusn, work )
789*
790* Do tests (1) to (4)
791*
792 CALL dget51( 1, mplusn, a, lda, ai, lda, q, lda, z, lda, work,
793 $ result( 1 ) )
794 CALL dget51( 1, mplusn, b, lda, bi, lda, q, lda, z, lda, work,
795 $ result( 2 ) )
796 CALL dget51( 3, mplusn, b, lda, bi, lda, q, lda, q, lda, work,
797 $ result( 3 ) )
798 CALL dget51( 3, mplusn, b, lda, bi, lda, z, lda, z, lda, work,
799 $ result( 4 ) )
800*
801* Do tests (5) and (6): check Schur form of A and compare
802* eigenvalues with diagonals.
803*
804 ntest = 6
805 temp1 = zero
806 result( 5 ) = zero
807 result( 6 ) = zero
808*
809 DO 110 j = 1, mplusn
810 ilabad = .false.
811 IF( alphai( j ).EQ.zero ) THEN
812 temp2 = ( abs( alphar( j )-ai( j, j ) ) /
813 $ max( smlnum, abs( alphar( j ) ), abs( ai( j,
814 $ j ) ) )+abs( beta( j )-bi( j, j ) ) /
815 $ max( smlnum, abs( beta( j ) ), abs( bi( j, j ) ) ) )
816 $ / ulp
817 IF( j.LT.mplusn ) THEN
818 IF( ai( j+1, j ).NE.zero ) THEN
819 ilabad = .true.
820 result( 5 ) = ulpinv
821 END IF
822 END IF
823 IF( j.GT.1 ) THEN
824 IF( ai( j, j-1 ).NE.zero ) THEN
825 ilabad = .true.
826 result( 5 ) = ulpinv
827 END IF
828 END IF
829 ELSE
830 IF( alphai( j ).GT.zero ) THEN
831 i1 = j
832 ELSE
833 i1 = j - 1
834 END IF
835 IF( i1.LE.0 .OR. i1.GE.mplusn ) THEN
836 ilabad = .true.
837 ELSE IF( i1.LT.mplusn-1 ) THEN
838 IF( ai( i1+2, i1+1 ).NE.zero ) THEN
839 ilabad = .true.
840 result( 5 ) = ulpinv
841 END IF
842 ELSE IF( i1.GT.1 ) THEN
843 IF( ai( i1, i1-1 ).NE.zero ) THEN
844 ilabad = .true.
845 result( 5 ) = ulpinv
846 END IF
847 END IF
848 IF( .NOT.ilabad ) THEN
849 CALL dget53( ai( i1, i1 ), lda, bi( i1, i1 ), lda,
850 $ beta( j ), alphar( j ), alphai( j ), temp2,
851 $ iinfo )
852 IF( iinfo.GE.3 ) THEN
853 WRITE( nout, fmt = 9997 )iinfo, j, mplusn, nptknt
854 info = abs( iinfo )
855 END IF
856 ELSE
857 temp2 = ulpinv
858 END IF
859 END IF
860 temp1 = max( temp1, temp2 )
861 IF( ilabad ) THEN
862 WRITE( nout, fmt = 9996 )j, mplusn, nptknt
863 END IF
864 110 CONTINUE
865 result( 6 ) = temp1
866*
867* Test (7) (if sorting worked) <--------- need to be checked.
868*
869 ntest = 7
870 result( 7 ) = zero
871 IF( linfo.EQ.mplusn+3 )
872 $ result( 7 ) = ulpinv
873*
874* Test (8): compare the estimated value of DIF and its true value.
875*
876 ntest = 8
877 result( 8 ) = zero
878 IF( difest( 2 ).EQ.zero ) THEN
879 IF( diftru.GT.abnrm*ulp )
880 $ result( 8 ) = ulpinv
881 ELSE IF( diftru.EQ.zero ) THEN
882 IF( difest( 2 ).GT.abnrm*ulp )
883 $ result( 8 ) = ulpinv
884 ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
885 $ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
886 result( 8 ) = max( diftru / difest( 2 ), difest( 2 ) / diftru )
887 END IF
888*
889* Test (9)
890*
891 ntest = 9
892 result( 9 ) = zero
893 IF( linfo.EQ.( mplusn+2 ) ) THEN
894 IF( diftru.GT.abnrm*ulp )
895 $ result( 9 ) = ulpinv
896 IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
897 $ result( 9 ) = ulpinv
898 IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
899 $ result( 9 ) = ulpinv
900 END IF
901*
902* Test (10): compare the estimated value of PL and it true value.
903*
904 ntest = 10
905 result( 10 ) = zero
906 IF( pl( 1 ).EQ.zero ) THEN
907 IF( pltru.GT.abnrm*ulp )
908 $ result( 10 ) = ulpinv
909 ELSE IF( pltru.EQ.zero ) THEN
910 IF( pl( 1 ).GT.abnrm*ulp )
911 $ result( 10 ) = ulpinv
912 ELSE IF( ( pltru.GT.thresh*pl( 1 ) ) .OR.
913 $ ( pltru*thresh.LT.pl( 1 ) ) ) THEN
914 result( 10 ) = ulpinv
915 END IF
916*
917 ntestt = ntestt + ntest
918*
919* Print out tests which fail.
920*
921 DO 120 j = 1, ntest
922 IF( result( j ).GE.thresh ) THEN
923*
924* If this is the first test to fail,
925* print a header to the data file.
926*
927 IF( nerrs.EQ.0 ) THEN
928 WRITE( nout, fmt = 9995 )'DGX'
929*
930* Matrix types
931*
932 WRITE( nout, fmt = 9994 )
933*
934* Tests performed
935*
936 WRITE( nout, fmt = 9992 )'orthogonal', '''',
937 $ 'transpose', ( '''', i = 1, 4 )
938*
939 END IF
940 nerrs = nerrs + 1
941 IF( result( j ).LT.10000.0d0 ) THEN
942 WRITE( nout, fmt = 9989 )nptknt, mplusn, j, result( j )
943 ELSE
944 WRITE( nout, fmt = 9988 )nptknt, mplusn, j, result( j )
945 END IF
946 END IF
947*
948 120 CONTINUE
949*
950 130 CONTINUE
951 GO TO 80
952 140 CONTINUE
953*
954 150 CONTINUE
955*
956* Summary
957*
958 CALL alasvm( 'DGX', nout, nerrs, ntestt, 0 )
959*
960 work( 1 ) = maxwrk
961*
962 RETURN
963*
964 9999 FORMAT( ' DDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
965 $ i6, ', JTYPE=', i6, ')' )
966*
967 9998 FORMAT( ' DDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
968 $ i6, ', Input Example #', i2, ')' )
969*
970 9997 FORMAT( ' DDRGSX: DGET53 returned INFO=', i1, ' for eigenvalue ',
971 $ i6, '.', / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
972*
973 9996 FORMAT( ' DDRGSX: S not in Schur form at eigenvalue ', i6, '.',
974 $ / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
975*
976 9995 FORMAT( / 1x, a3, ' -- Real Expert Generalized Schur form',
977 $ ' problem driver' )
978*
979 9994 FORMAT( 'Input Example' )
980*
981 9993 FORMAT( ' Matrix types: ', /
982 $ ' 1: A is a block diagonal matrix of Jordan blocks ',
983 $ 'and B is the identity ', / ' matrix, ',
984 $ / ' 2: A and B are upper triangular matrices, ',
985 $ / ' 3: A and B are as type 2, but each second diagonal ',
986 $ 'block in A_11 and ', /
987 $ ' each third diaongal block in A_22 are 2x2 blocks,',
988 $ / ' 4: A and B are block diagonal matrices, ',
989 $ / ' 5: (A,B) has potentially close or common ',
990 $ 'eigenvalues.', / )
991*
992 9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
993 $ 'Q and Z are ', a, ',', / 19x,
994 $ ' a is alpha, b is beta, and ', a, ' means ', a, '.)',
995 $ / ' 1 = | A - Q S Z', a,
996 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
997 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
998 $ ' | / ( n ulp ) 4 = | I - ZZ', a,
999 $ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ',
1000 $ 'Schur form S', / ' 6 = difference between (alpha,beta)',
1001 $ ' and diagonals of (S,T)', /
1002 $ ' 7 = 1/ULP if SDIM is not the correct number of ',
1003 $ 'selected eigenvalues', /
1004 $ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
1005 $ 'DIFTRU/DIFEST > 10*THRESH',
1006 $ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
1007 $ 'when reordering fails', /
1008 $ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
1009 $ 'PLTRU/PLEST > THRESH', /
1010 $ ' ( Test 10 is only for input examples )', / )
1011 9991 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', d10.3,
1012 $ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, f8.2 )
1013 9990 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', d10.3,
1014 $ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, d10.3 )
1015 9989 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
1016 $ ' result ', i2, ' is', 0p, f8.2 )
1017 9988 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
1018 $ ' result ', i2, ' is', 1p, d10.3 )
1019*
1020* End of DDRGSX
1021*
1022 END
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
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
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:73
logical function dlctsx(AR, AI, BETA)
DLCTSX
Definition: dlctsx.f:65
subroutine dget53(A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO)
DGET53
Definition: dget53.f:126
subroutine ddrgsx(NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI, BI, Z, Q, ALPHAR, ALPHAI, BETA, C, LDC, S, WORK, LWORK, IWORK, LIWORK, BWORK, INFO)
DDRGSX
Definition: ddrgsx.f:359
subroutine dget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RESULT)
DGET51
Definition: dget51.f:149
subroutine dlatm5(PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA, QBLCKB)
DLATM5
Definition: dlatm5.f:268
subroutine dlakf2(M, N, A, LDA, B, D, E, Z, LDZ)
DLAKF2
Definition: dlakf2.f:105
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 dggesx(JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA, B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, IWORK, LIWORK, BWORK, INFO)
DGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
Definition: dggesx.f:365
subroutine dgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO)
DGESVD computes the singular value decomposition (SVD) for GE matrices
Definition: dgesvd.f:211