LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sdrgsx.f
Go to the documentation of this file.
1*> \brief \b SDRGSX
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 SDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B,
12* AI, 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* REAL THRESH
19* ..
20* .. Array Arguments ..
21* LOGICAL BWORK( * )
22* INTEGER IWORK( * )
23* REAL 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*> SDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
36*> problem expert driver SGGESX.
37*>
38*> SGGESX 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 SDRGSX is called with NSIZE > 0, five (5) types of built-in
57*> matrix pairs are used to test the routine SGGESX.
58*>
59*> When SDRGSX is called with NSIZE = 0, it reads in test matrix data
60*> to test SGGESX.
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 SGGESX, 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 SGESVD) 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 SLATM5).
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 SGGESX.
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 SLAKF2
199*> \endverbatim
200*>
201*> \param[in] THRESH
202*> \verbatim
203*> THRESH is REAL
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 REAL 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 REAL 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 REAL array, dimension (LDA, NSIZE)
251*> Copy of A, modified by SGGESX.
252*> \endverbatim
253*>
254*> \param[out] BI
255*> \verbatim
256*> BI is REAL array, dimension (LDA, NSIZE)
257*> Copy of B, modified by SGGESX.
258*> \endverbatim
259*>
260*> \param[out] Z
261*> \verbatim
262*> Z is REAL array, dimension (LDA, NSIZE)
263*> Z holds the left Schur vectors computed by SGGESX.
264*> \endverbatim
265*>
266*> \param[out] Q
267*> \verbatim
268*> Q is REAL array, dimension (LDA, NSIZE)
269*> Q holds the right Schur vectors computed by SGGESX.
270*> \endverbatim
271*>
272*> \param[out] ALPHAR
273*> \verbatim
274*> ALPHAR is REAL array, dimension (NSIZE)
275*> \endverbatim
276*>
277*> \param[out] ALPHAI
278*> \verbatim
279*> ALPHAI is REAL array, dimension (NSIZE)
280*> \endverbatim
281*>
282*> \param[out] BETA
283*> \verbatim
284*> BETA is REAL array, dimension (NSIZE)
285*> \verbatim
286*> On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues.
287*> \endverbatim
288*>
289*> \param[out] C
290*> \verbatim
291*> C is REAL array, dimension (LDC, LDC)
292*> Store the matrix generated by subroutine SLAKF2, 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 REAL array, dimension (LDC)
306*> Singular values of C
307*> \endverbatim
308*>
309*> \param[out] WORK
310*> \verbatim
311*> WORK is REAL 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 single_eig
354*
355* =====================================================================
356 SUBROUTINE sdrgsx( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B,
357 $ AI, 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 REAL THRESH
368* ..
369* .. Array Arguments ..
370 LOGICAL BWORK( * )
371 INTEGER IWORK( * )
372 REAL 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 REAL ZERO, ONE, TEN
382 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0, ten = 1.0e+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 REAL ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
391 $ TEMP2, THRSH2, ULP, ULPINV, WEIGHT
392* ..
393* .. Local Arrays ..
394 REAL DIFEST( 2 ), PL( 2 ), RESULT( 10 )
395* ..
396* .. External Functions ..
397 LOGICAL SLCTSX
398 INTEGER ILAENV
399 REAL SLAMCH, SLANGE
400 EXTERNAL slctsx, ilaenv, slamch, slange
401* ..
402* .. External Subroutines ..
403 EXTERNAL alasvm, sgesvd, sget51, sget53, sggesx, slabad,
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
445c MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2-2 )
446 minwrk = max( 10*( nsize+1 ), 5*nsize*nsize / 2 )
447*
448* workspace for sggesx
449*
450 maxwrk = 9*( nsize+1 ) + nsize*
451 $ ilaenv( 1, 'SGEQRF', ' ', nsize, 1, nsize, 0 )
452 maxwrk = max( maxwrk, 9*( nsize+1 )+nsize*
453 $ ilaenv( 1, 'SORGQR', ' ', nsize, 1, nsize, -1 ) )
454*
455* workspace for sgesvd
456*
457 bdspac = 5*nsize*nsize / 2
458 maxwrk = max( maxwrk, 3*nsize*nsize / 2+nsize*nsize*
459 $ ilaenv( 1, 'SGEBRD', ' ', 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 = -19
470*
471 IF( info.NE.0 ) THEN
472 CALL xerbla( 'SDRGSX', -info )
473 RETURN
474 END IF
475*
476* Important constants
477*
478 ulp = slamch( 'P' )
479 ulpinv = one / ulp
480 smlnum = slamch( 'S' ) / ulp
481 bignum = one / smlnum
482 CALL slabad( 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 SGGESX, 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 slaset( 'Full', mplusn, mplusn, zero, zero, ai,
516 $ lda )
517 CALL slaset( 'Full', mplusn, mplusn, zero, zero, bi,
518 $ lda )
519*
520 CALL slatm5( 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 SLCTSX
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 slacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
541 CALL slacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
542*
543 CALL sggesx( 'V', 'V', 'S', slctsx, sense, mplusn, ai,
544 $ lda, bi, lda, mm, alphar, alphai, beta,
545 $ q, lda, z, lda, pl, difest, work, lwork,
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 )'SGGESX', linfo, mplusn,
551 $ prtype
552 info = linfo
553 GO TO 30
554 END IF
555*
556* Compute the norm(A, B)
557*
558 CALL slacpy( 'Full', mplusn, mplusn, ai, lda, work,
559 $ mplusn )
560 CALL slacpy( 'Full', mplusn, mplusn, bi, lda,
561 $ work( mplusn*mplusn+1 ), mplusn )
562 abnrm = slange( 'Fro', mplusn, 2*mplusn, work, mplusn,
563 $ work )
564*
565* Do tests (1) to (4)
566*
567 CALL sget51( 1, mplusn, a, lda, ai, lda, q, lda, z,
568 $ lda, work, result( 1 ) )
569 CALL sget51( 1, mplusn, b, lda, bi, lda, q, lda, z,
570 $ lda, work, result( 2 ) )
571 CALL sget51( 3, mplusn, b, lda, bi, lda, q, lda, q,
572 $ lda, work, result( 3 ) )
573 CALL sget51( 3, mplusn, b, lda, bi, lda, z, lda, z,
574 $ lda, work, result( 4 ) )
575 ntest = 4
576*
577* Do tests (5) and (6): check Schur form of A and
578* compare eigenvalues with diagonals.
579*
580 temp1 = zero
581 result( 5 ) = zero
582 result( 6 ) = zero
583*
584 DO 10 j = 1, mplusn
585 ilabad = .false.
586 IF( alphai( j ).EQ.zero ) THEN
587 temp2 = ( abs( alphar( j )-ai( j, j ) ) /
588 $ max( smlnum, abs( alphar( j ) ),
589 $ abs( ai( j, j ) ) )+
590 $ abs( beta( j )-bi( j, j ) ) /
591 $ max( smlnum, abs( beta( j ) ),
592 $ abs( 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 ELSE
606 IF( alphai( j ).GT.zero ) THEN
607 i1 = j
608 ELSE
609 i1 = j - 1
610 END IF
611 IF( i1.LE.0 .OR. i1.GE.mplusn ) THEN
612 ilabad = .true.
613 ELSE IF( i1.LT.mplusn-1 ) THEN
614 IF( ai( i1+2, i1+1 ).NE.zero ) THEN
615 ilabad = .true.
616 result( 5 ) = ulpinv
617 END IF
618 ELSE IF( i1.GT.1 ) THEN
619 IF( ai( i1, i1-1 ).NE.zero ) THEN
620 ilabad = .true.
621 result( 5 ) = ulpinv
622 END IF
623 END IF
624 IF( .NOT.ilabad ) THEN
625 CALL sget53( ai( i1, i1 ), lda, bi( i1, i1 ),
626 $ lda, beta( j ), alphar( j ),
627 $ alphai( j ), temp2, iinfo )
628 IF( iinfo.GE.3 ) THEN
629 WRITE( nout, fmt = 9997 )iinfo, j,
630 $ mplusn, prtype
631 info = abs( iinfo )
632 END IF
633 ELSE
634 temp2 = ulpinv
635 END IF
636 END IF
637 temp1 = max( temp1, temp2 )
638 IF( ilabad ) THEN
639 WRITE( nout, fmt = 9996 )j, mplusn, prtype
640 END IF
641 10 CONTINUE
642 result( 6 ) = temp1
643 ntest = ntest + 2
644*
645* Test (7) (if sorting worked)
646*
647 result( 7 ) = zero
648 IF( linfo.EQ.mplusn+3 ) THEN
649 result( 7 ) = ulpinv
650 ELSE IF( mm.NE.n ) THEN
651 result( 7 ) = ulpinv
652 END IF
653 ntest = ntest + 1
654*
655* Test (8): compare the estimated value DIF and its
656* value. first, compute the exact DIF.
657*
658 result( 8 ) = zero
659 mn2 = mm*( mplusn-mm )*2
660 IF( ifunc.GE.2 .AND. mn2.LE.ncmax*ncmax ) THEN
661*
662* Note: for either following two causes, there are
663* almost same number of test cases fail the test.
664*
665 CALL slakf2( mm, mplusn-mm, ai, lda,
666 $ ai( mm+1, mm+1 ), bi,
667 $ bi( mm+1, mm+1 ), c, ldc )
668*
669 CALL sgesvd( 'N', 'N', mn2, mn2, c, ldc, s, work,
670 $ 1, work( 2 ), 1, work( 3 ), lwork-2,
671 $ info )
672 diftru = s( mn2 )
673*
674 IF( difest( 2 ).EQ.zero ) THEN
675 IF( diftru.GT.abnrm*ulp )
676 $ result( 8 ) = ulpinv
677 ELSE IF( diftru.EQ.zero ) THEN
678 IF( difest( 2 ).GT.abnrm*ulp )
679 $ result( 8 ) = ulpinv
680 ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
681 $ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
682 result( 8 ) = max( diftru / difest( 2 ),
683 $ difest( 2 ) / diftru )
684 END IF
685 ntest = ntest + 1
686 END IF
687*
688* Test (9)
689*
690 result( 9 ) = zero
691 IF( linfo.EQ.( mplusn+2 ) ) THEN
692 IF( diftru.GT.abnrm*ulp )
693 $ result( 9 ) = ulpinv
694 IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
695 $ result( 9 ) = ulpinv
696 IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
697 $ result( 9 ) = ulpinv
698 ntest = ntest + 1
699 END IF
700*
701 ntestt = ntestt + ntest
702*
703* Print out tests which fail.
704*
705 DO 20 j = 1, 9
706 IF( result( j ).GE.thresh ) THEN
707*
708* If this is the first test to fail,
709* print a header to the data file.
710*
711 IF( nerrs.EQ.0 ) THEN
712 WRITE( nout, fmt = 9995 )'SGX'
713*
714* Matrix types
715*
716 WRITE( nout, fmt = 9993 )
717*
718* Tests performed
719*
720 WRITE( nout, fmt = 9992 )'orthogonal', '''',
721 $ 'transpose', ( '''', i = 1, 4 )
722*
723 END IF
724 nerrs = nerrs + 1
725 IF( result( j ).LT.10000.0 ) THEN
726 WRITE( nout, fmt = 9991 )mplusn, prtype,
727 $ weight, m, j, result( j )
728 ELSE
729 WRITE( nout, fmt = 9990 )mplusn, prtype,
730 $ weight, m, j, result( j )
731 END IF
732 END IF
733 20 CONTINUE
734*
735 30 CONTINUE
736 40 CONTINUE
737 50 CONTINUE
738 60 CONTINUE
739*
740 GO TO 150
741*
742 70 CONTINUE
743*
744* Read in data from file to check accuracy of condition estimation
745* Read input data until N=0
746*
747 nptknt = 0
748*
749 80 CONTINUE
750 READ( nin, fmt = *, END = 140 )mplusn
751 IF( mplusn.EQ.0 )
752 $ GO TO 140
753 READ( nin, fmt = *, END = 140 )n
754 DO 90 i = 1, mplusn
755 READ( nin, fmt = * )( ai( i, j ), j = 1, mplusn )
756 90 CONTINUE
757 DO 100 i = 1, mplusn
758 READ( nin, fmt = * )( bi( i, j ), j = 1, mplusn )
759 100 CONTINUE
760 READ( nin, fmt = * )pltru, diftru
761*
762 nptknt = nptknt + 1
763 fs = .true.
764 k = 0
765 m = mplusn - n
766*
767 CALL slacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
768 CALL slacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
769*
770* Compute the Schur factorization while swapping the
771* m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
772*
773 CALL sggesx( 'V', 'V', 'S', slctsx, 'B', mplusn, ai, lda, bi, lda,
774 $ mm, alphar, alphai, beta, q, lda, z, lda, pl, difest,
775 $ work, lwork, iwork, liwork, bwork, linfo )
776*
777 IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
778 result( 1 ) = ulpinv
779 WRITE( nout, fmt = 9998 )'SGGESX', linfo, mplusn, nptknt
780 GO TO 130
781 END IF
782*
783* Compute the norm(A, B)
784* (should this be norm of (A,B) or (AI,BI)?)
785*
786 CALL slacpy( 'Full', mplusn, mplusn, ai, lda, work, mplusn )
787 CALL slacpy( 'Full', mplusn, mplusn, bi, lda,
788 $ work( mplusn*mplusn+1 ), mplusn )
789 abnrm = slange( 'Fro', mplusn, 2*mplusn, work, mplusn, work )
790*
791* Do tests (1) to (4)
792*
793 CALL sget51( 1, mplusn, a, lda, ai, lda, q, lda, z, lda, work,
794 $ result( 1 ) )
795 CALL sget51( 1, mplusn, b, lda, bi, lda, q, lda, z, lda, work,
796 $ result( 2 ) )
797 CALL sget51( 3, mplusn, b, lda, bi, lda, q, lda, q, lda, work,
798 $ result( 3 ) )
799 CALL sget51( 3, mplusn, b, lda, bi, lda, z, lda, z, lda, work,
800 $ result( 4 ) )
801*
802* Do tests (5) and (6): check Schur form of A and compare
803* eigenvalues with diagonals.
804*
805 ntest = 6
806 temp1 = zero
807 result( 5 ) = zero
808 result( 6 ) = zero
809*
810 DO 110 j = 1, mplusn
811 ilabad = .false.
812 IF( alphai( j ).EQ.zero ) THEN
813 temp2 = ( abs( alphar( j )-ai( j, j ) ) /
814 $ max( smlnum, abs( alphar( j ) ), abs( ai( j,
815 $ j ) ) )+abs( beta( j )-bi( j, j ) ) /
816 $ max( smlnum, abs( beta( j ) ), abs( bi( j, j ) ) ) )
817 $ / ulp
818 IF( j.LT.mplusn ) THEN
819 IF( ai( j+1, j ).NE.zero ) THEN
820 ilabad = .true.
821 result( 5 ) = ulpinv
822 END IF
823 END IF
824 IF( j.GT.1 ) THEN
825 IF( ai( j, j-1 ).NE.zero ) THEN
826 ilabad = .true.
827 result( 5 ) = ulpinv
828 END IF
829 END IF
830 ELSE
831 IF( alphai( j ).GT.zero ) THEN
832 i1 = j
833 ELSE
834 i1 = j - 1
835 END IF
836 IF( i1.LE.0 .OR. i1.GE.mplusn ) THEN
837 ilabad = .true.
838 ELSE IF( i1.LT.mplusn-1 ) THEN
839 IF( ai( i1+2, i1+1 ).NE.zero ) THEN
840 ilabad = .true.
841 result( 5 ) = ulpinv
842 END IF
843 ELSE IF( i1.GT.1 ) THEN
844 IF( ai( i1, i1-1 ).NE.zero ) THEN
845 ilabad = .true.
846 result( 5 ) = ulpinv
847 END IF
848 END IF
849 IF( .NOT.ilabad ) THEN
850 CALL sget53( ai( i1, i1 ), lda, bi( i1, i1 ), lda,
851 $ beta( j ), alphar( j ), alphai( j ), temp2,
852 $ iinfo )
853 IF( iinfo.GE.3 ) THEN
854 WRITE( nout, fmt = 9997 )iinfo, j, mplusn, nptknt
855 info = abs( iinfo )
856 END IF
857 ELSE
858 temp2 = ulpinv
859 END IF
860 END IF
861 temp1 = max( temp1, temp2 )
862 IF( ilabad ) THEN
863 WRITE( nout, fmt = 9996 )j, mplusn, nptknt
864 END IF
865 110 CONTINUE
866 result( 6 ) = temp1
867*
868* Test (7) (if sorting worked) <--------- need to be checked.
869*
870 ntest = 7
871 result( 7 ) = zero
872 IF( linfo.EQ.mplusn+3 )
873 $ result( 7 ) = ulpinv
874*
875* Test (8): compare the estimated value of DIF and its true value.
876*
877 ntest = 8
878 result( 8 ) = zero
879 IF( difest( 2 ).EQ.zero ) THEN
880 IF( diftru.GT.abnrm*ulp )
881 $ result( 8 ) = ulpinv
882 ELSE IF( diftru.EQ.zero ) THEN
883 IF( difest( 2 ).GT.abnrm*ulp )
884 $ result( 8 ) = ulpinv
885 ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
886 $ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
887 result( 8 ) = max( diftru / difest( 2 ), difest( 2 ) / diftru )
888 END IF
889*
890* Test (9)
891*
892 ntest = 9
893 result( 9 ) = zero
894 IF( linfo.EQ.( mplusn+2 ) ) THEN
895 IF( diftru.GT.abnrm*ulp )
896 $ result( 9 ) = ulpinv
897 IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
898 $ result( 9 ) = ulpinv
899 IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
900 $ result( 9 ) = ulpinv
901 END IF
902*
903* Test (10): compare the estimated value of PL and it true value.
904*
905 ntest = 10
906 result( 10 ) = zero
907 IF( pl( 1 ).EQ.zero ) THEN
908 IF( pltru.GT.abnrm*ulp )
909 $ result( 10 ) = ulpinv
910 ELSE IF( pltru.EQ.zero ) THEN
911 IF( pl( 1 ).GT.abnrm*ulp )
912 $ result( 10 ) = ulpinv
913 ELSE IF( ( pltru.GT.thresh*pl( 1 ) ) .OR.
914 $ ( pltru*thresh.LT.pl( 1 ) ) ) THEN
915 result( 10 ) = ulpinv
916 END IF
917*
918 ntestt = ntestt + ntest
919*
920* Print out tests which fail.
921*
922 DO 120 j = 1, ntest
923 IF( result( j ).GE.thresh ) THEN
924*
925* If this is the first test to fail,
926* print a header to the data file.
927*
928 IF( nerrs.EQ.0 ) THEN
929 WRITE( nout, fmt = 9995 )'SGX'
930*
931* Matrix types
932*
933 WRITE( nout, fmt = 9994 )
934*
935* Tests performed
936*
937 WRITE( nout, fmt = 9992 )'orthogonal', '''',
938 $ 'transpose', ( '''', i = 1, 4 )
939*
940 END IF
941 nerrs = nerrs + 1
942 IF( result( j ).LT.10000.0 ) THEN
943 WRITE( nout, fmt = 9989 )nptknt, mplusn, j, result( j )
944 ELSE
945 WRITE( nout, fmt = 9988 )nptknt, mplusn, j, result( j )
946 END IF
947 END IF
948*
949 120 CONTINUE
950*
951 130 CONTINUE
952 GO TO 80
953 140 CONTINUE
954*
955 150 CONTINUE
956*
957* Summary
958*
959 CALL alasvm( 'SGX', nout, nerrs, ntestt, 0 )
960*
961 work( 1 ) = maxwrk
962*
963 RETURN
964*
965 9999 FORMAT( ' SDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
966 $ i6, ', JTYPE=', i6, ')' )
967*
968 9998 FORMAT( ' SDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
969 $ i6, ', Input Example #', i2, ')' )
970*
971 9997 FORMAT( ' SDRGSX: SGET53 returned INFO=', i1, ' for eigenvalue ',
972 $ i6, '.', / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
973*
974 9996 FORMAT( ' SDRGSX: S not in Schur form at eigenvalue ', i6, '.',
975 $ / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
976*
977 9995 FORMAT( / 1x, a3, ' -- Real Expert Generalized Schur form',
978 $ ' problem driver' )
979*
980 9994 FORMAT( 'Input Example' )
981*
982 9993 FORMAT( ' Matrix types: ', /
983 $ ' 1: A is a block diagonal matrix of Jordan blocks ',
984 $ 'and B is the identity ', / ' matrix, ',
985 $ / ' 2: A and B are upper triangular matrices, ',
986 $ / ' 3: A and B are as type 2, but each second diagonal ',
987 $ 'block in A_11 and ', /
988 $ ' each third diaongal block in A_22 are 2x2 blocks,',
989 $ / ' 4: A and B are block diagonal matrices, ',
990 $ / ' 5: (A,B) has potentially close or common ',
991 $ 'eigenvalues.', / )
992*
993 9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
994 $ 'Q and Z are ', a, ',', / 19x,
995 $ ' a is alpha, b is beta, and ', a, ' means ', a, '.)',
996 $ / ' 1 = | A - Q S Z', a,
997 $ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
998 $ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
999 $ ' | / ( n ulp ) 4 = | I - ZZ', a,
1000 $ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ',
1001 $ 'Schur form S', / ' 6 = difference between (alpha,beta)',
1002 $ ' and diagonals of (S,T)', /
1003 $ ' 7 = 1/ULP if SDIM is not the correct number of ',
1004 $ 'selected eigenvalues', /
1005 $ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
1006 $ 'DIFTRU/DIFEST > 10*THRESH',
1007 $ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
1008 $ 'when reordering fails', /
1009 $ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
1010 $ 'PLTRU/PLEST > THRESH', /
1011 $ ' ( Test 10 is only for input examples )', / )
1012 9991 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', e10.3,
1013 $ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, f8.2 )
1014 9990 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', e10.3,
1015 $ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, e10.3 )
1016 9989 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
1017 $ ' result ', i2, ' is', 0p, f8.2 )
1018 9988 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
1019 $ ' result ', i2, ' is', 1p, e10.3 )
1020*
1021* End of SDRGSX
1022*
1023 END
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:73
subroutine slakf2(M, N, A, LDA, B, D, E, Z, LDZ)
SLAKF2
Definition: slakf2.f:105
subroutine slatm5(PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA, QBLCKB)
SLATM5
Definition: slatm5.f:268
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: slange.f:114
subroutine sggesx(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)
SGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
Definition: sggesx.f:365
subroutine sgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO)
SGESVD computes the singular value decomposition (SVD) for GE matrices
Definition: sgesvd.f:211
subroutine sget53(A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO)
SGET53
Definition: sget53.f:126
logical function slctsx(AR, AI, BETA)
SLCTSX
Definition: slctsx.f:65
subroutine sget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RESULT)
SGET51
Definition: sget51.f:149
subroutine sdrgsx(NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI, BI, Z, Q, ALPHAR, ALPHAI, BETA, C, LDC, S, WORK, LWORK, IWORK, LIWORK, BWORK, INFO)
SDRGSX
Definition: sdrgsx.f:359