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