LAPACK  3.6.1 LAPACK: Linear Algebra PACKage
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 *> \date June 2016
354 *
355 *> \ingroup double_eig
356 *
357 * =====================================================================
358  SUBROUTINE ddrgsx( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI,
359  \$ bi, z, q, alphar, alphai, beta, c, ldc, s,
360  \$ work, lwork, iwork, liwork, bwork, info )
361 *
362 * -- LAPACK test routine (version 3.6.1) --
363 * -- LAPACK is a software package provided by Univ. of Tennessee, --
364 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
365 * June 2016
366 *
367 * .. Scalar Arguments ..
368  INTEGER INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
369  \$ nout, nsize
370  DOUBLE PRECISION THRESH
371 * ..
372 * .. Array Arguments ..
373  LOGICAL BWORK( * )
374  INTEGER IWORK( * )
375  DOUBLE PRECISION A( lda, * ), AI( lda, * ), ALPHAI( * ),
376  \$ alphar( * ), b( lda, * ), beta( * ),
377  \$ bi( lda, * ), c( ldc, * ), q( lda, * ), s( * ),
378  \$ work( * ), z( lda, * )
379 * ..
380 *
381 * =====================================================================
382 *
383 * .. Parameters ..
384  DOUBLE PRECISION ZERO, ONE, TEN
385  parameter ( zero = 0.0d+0, one = 1.0d+0, ten = 1.0d+1 )
386 * ..
387 * .. Local Scalars ..
388  LOGICAL ILABAD
389  CHARACTER SENSE
390  INTEGER BDSPAC, I, I1, IFUNC, IINFO, J, LINFO, MAXWRK,
391  \$ minwrk, mm, mn2, nerrs, nptknt, ntest, ntestt,
392  \$ prtype, qba, qbb
393  DOUBLE PRECISION ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
394  \$ temp2, thrsh2, ulp, ulpinv, weight
395 * ..
396 * .. Local Arrays ..
397  DOUBLE PRECISION DIFEST( 2 ), PL( 2 ), RESULT( 10 )
398 * ..
399 * .. External Functions ..
400  LOGICAL DLCTSX
401  INTEGER ILAENV
402  DOUBLE PRECISION DLAMCH, DLANGE
403  EXTERNAL dlctsx, ilaenv, dlamch, dlange
404 * ..
405 * .. External Subroutines ..
406  EXTERNAL alasvm, dgesvd, dget51, dget53, dggesx, dlabad,
408 * ..
409 * .. Intrinsic Functions ..
410  INTRINSIC abs, max, sqrt
411 * ..
412 * .. Scalars in Common ..
413  LOGICAL FS
414  INTEGER K, M, MPLUSN, N
415 * ..
416 * .. Common blocks ..
417  COMMON / mn / m, n, mplusn, k, fs
418 * ..
419 * .. Executable Statements ..
420 *
421 * Check for errors
422 *
423  IF( nsize.LT.0 ) THEN
424  info = -1
425  ELSE IF( thresh.LT.zero ) THEN
426  info = -2
427  ELSE IF( nin.LE.0 ) THEN
428  info = -3
429  ELSE IF( nout.LE.0 ) THEN
430  info = -4
431  ELSE IF( lda.LT.1 .OR. lda.LT.nsize ) THEN
432  info = -6
433  ELSE IF( ldc.LT.1 .OR. ldc.LT.nsize*nsize / 2 ) THEN
434  info = -17
435  ELSE IF( liwork.LT.nsize+6 ) THEN
436  info = -21
437  END IF
438 *
439 * Compute workspace
440 * (Note: Comments in the code beginning "Workspace:" describe the
441 * minimal amount of workspace needed at that point in the code,
442 * as well as the preferred amount for good performance.
443 * NB refers to the optimal block size for the immediately
444 * following subroutine, as returned by ILAENV.)
445 *
446  minwrk = 1
447  IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
448  minwrk = max( 10*( nsize+1 ), 5*nsize*nsize / 2 )
449 *
450 * workspace for sggesx
451 *
452  maxwrk = 9*( nsize+1 ) + nsize*
453  \$ ilaenv( 1, 'DGEQRF', ' ', nsize, 1, nsize, 0 )
454  maxwrk = max( maxwrk, 9*( nsize+1 )+nsize*
455  \$ ilaenv( 1, 'DORGQR', ' ', nsize, 1, nsize, -1 ) )
456 *
457 * workspace for dgesvd
458 *
459  bdspac = 5*nsize*nsize / 2
460  maxwrk = max( maxwrk, 3*nsize*nsize / 2+nsize*nsize*
461  \$ ilaenv( 1, 'DGEBRD', ' ', nsize*nsize / 2,
462  \$ nsize*nsize / 2, -1, -1 ) )
463  maxwrk = max( maxwrk, bdspac )
464 *
465  maxwrk = max( maxwrk, minwrk )
466 *
467  work( 1 ) = maxwrk
468  END IF
469 *
470  IF( lwork.LT.minwrk )
471  \$ info = -19
472 *
473  IF( info.NE.0 ) THEN
474  CALL xerbla( 'DDRGSX', -info )
475  RETURN
476  END IF
477 *
478 * Important constants
479 *
480  ulp = dlamch( 'P' )
481  ulpinv = one / ulp
482  smlnum = dlamch( 'S' ) / ulp
483  bignum = one / smlnum
484  CALL dlabad( smlnum, bignum )
485  thrsh2 = ten*thresh
486  ntestt = 0
487  nerrs = 0
488 *
489 * Go to the tests for read-in matrix pairs
490 *
491  ifunc = 0
492  IF( nsize.EQ.0 )
493  \$ GO TO 70
494 *
495 * Test the built-in matrix pairs.
496 * Loop over different functions (IFUNC) of DGGESX, types (PRTYPE)
497 * of test matrices, different size (M+N)
498 *
499  prtype = 0
500  qba = 3
501  qbb = 4
502  weight = sqrt( ulp )
503 *
504  DO 60 ifunc = 0, 3
505  DO 50 prtype = 1, 5
506  DO 40 m = 1, nsize - 1
507  DO 30 n = 1, nsize - m
508 *
509  weight = one / weight
510  mplusn = m + n
511 *
512 * Generate test matrices
513 *
514  fs = .true.
515  k = 0
516 *
517  CALL dlaset( 'Full', mplusn, mplusn, zero, zero, ai,
518  \$ lda )
519  CALL dlaset( 'Full', mplusn, mplusn, zero, zero, bi,
520  \$ lda )
521 *
522  CALL dlatm5( prtype, m, n, ai, lda, ai( m+1, m+1 ),
523  \$ lda, ai( 1, m+1 ), lda, bi, lda,
524  \$ bi( m+1, m+1 ), lda, bi( 1, m+1 ), lda,
525  \$ q, lda, z, lda, weight, qba, qbb )
526 *
527 * Compute the Schur factorization and swapping the
528 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
529 * Swapping is accomplished via the function DLCTSX
530 * which is supplied below.
531 *
532  IF( ifunc.EQ.0 ) THEN
533  sense = 'N'
534  ELSE IF( ifunc.EQ.1 ) THEN
535  sense = 'E'
536  ELSE IF( ifunc.EQ.2 ) THEN
537  sense = 'V'
538  ELSE IF( ifunc.EQ.3 ) THEN
539  sense = 'B'
540  END IF
541 *
542  CALL dlacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
543  CALL dlacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
544 *
545  CALL dggesx( 'V', 'V', 'S', dlctsx, sense, mplusn, ai,
546  \$ lda, bi, lda, mm, alphar, alphai, beta,
547  \$ q, lda, z, lda, pl, difest, work, lwork,
548  \$ iwork, liwork, bwork, linfo )
549 *
550  IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
551  result( 1 ) = ulpinv
552  WRITE( nout, fmt = 9999 )'DGGESX', linfo, mplusn,
553  \$ prtype
554  info = linfo
555  GO TO 30
556  END IF
557 *
558 * Compute the norm(A, B)
559 *
560  CALL dlacpy( 'Full', mplusn, mplusn, ai, lda, work,
561  \$ mplusn )
562  CALL dlacpy( 'Full', mplusn, mplusn, bi, lda,
563  \$ work( mplusn*mplusn+1 ), mplusn )
564  abnrm = dlange( 'Fro', mplusn, 2*mplusn, work, mplusn,
565  \$ work )
566 *
567 * Do tests (1) to (4)
568 *
569  CALL dget51( 1, mplusn, a, lda, ai, lda, q, lda, z,
570  \$ lda, work, result( 1 ) )
571  CALL dget51( 1, mplusn, b, lda, bi, lda, q, lda, z,
572  \$ lda, work, result( 2 ) )
573  CALL dget51( 3, mplusn, b, lda, bi, lda, q, lda, q,
574  \$ lda, work, result( 3 ) )
575  CALL dget51( 3, mplusn, b, lda, bi, lda, z, lda, z,
576  \$ lda, work, result( 4 ) )
577  ntest = 4
578 *
579 * Do tests (5) and (6): check Schur form of A and
580 * compare eigenvalues with diagonals.
581 *
582  temp1 = zero
583  result( 5 ) = zero
584  result( 6 ) = zero
585 *
586  DO 10 j = 1, mplusn
587  ilabad = .false.
588  IF( alphai( j ).EQ.zero ) THEN
589  temp2 = ( abs( alphar( j )-ai( j, j ) ) /
590  \$ max( smlnum, abs( alphar( j ) ),
591  \$ abs( ai( j, j ) ) )+
592  \$ abs( beta( j )-bi( j, j ) ) /
593  \$ max( smlnum, abs( beta( j ) ),
594  \$ abs( bi( j, j ) ) ) ) / ulp
595  IF( j.LT.mplusn ) THEN
596  IF( ai( j+1, j ).NE.zero ) THEN
597  ilabad = .true.
598  result( 5 ) = ulpinv
599  END IF
600  END IF
601  IF( j.GT.1 ) THEN
602  IF( ai( j, j-1 ).NE.zero ) THEN
603  ilabad = .true.
604  result( 5 ) = ulpinv
605  END IF
606  END IF
607  ELSE
608  IF( alphai( j ).GT.zero ) THEN
609  i1 = j
610  ELSE
611  i1 = j - 1
612  END IF
613  IF( i1.LE.0 .OR. i1.GE.mplusn ) THEN
614  ilabad = .true.
615  ELSE IF( i1.LT.mplusn-1 ) THEN
616  IF( ai( i1+2, i1+1 ).NE.zero ) THEN
617  ilabad = .true.
618  result( 5 ) = ulpinv
619  END IF
620  ELSE IF( i1.GT.1 ) THEN
621  IF( ai( i1, i1-1 ).NE.zero ) THEN
622  ilabad = .true.
623  result( 5 ) = ulpinv
624  END IF
625  END IF
626  IF( .NOT.ilabad ) THEN
627  CALL dget53( ai( i1, i1 ), lda, bi( i1, i1 ),
628  \$ lda, beta( j ), alphar( j ),
629  \$ alphai( j ), temp2, iinfo )
630  IF( iinfo.GE.3 ) THEN
631  WRITE( nout, fmt = 9997 )iinfo, j,
632  \$ mplusn, prtype
633  info = abs( iinfo )
634  END IF
635  ELSE
636  temp2 = ulpinv
637  END IF
638  END IF
639  temp1 = max( temp1, temp2 )
640  IF( ilabad ) THEN
641  WRITE( nout, fmt = 9996 )j, mplusn, prtype
642  END IF
643  10 CONTINUE
644  result( 6 ) = temp1
645  ntest = ntest + 2
646 *
647 * Test (7) (if sorting worked)
648 *
649  result( 7 ) = zero
650  IF( linfo.EQ.mplusn+3 ) THEN
651  result( 7 ) = ulpinv
652  ELSE IF( mm.NE.n ) THEN
653  result( 7 ) = ulpinv
654  END IF
655  ntest = ntest + 1
656 *
657 * Test (8): compare the estimated value DIF and its
658 * value. first, compute the exact DIF.
659 *
660  result( 8 ) = zero
661  mn2 = mm*( mplusn-mm )*2
662  IF( ifunc.GE.2 .AND. mn2.LE.ncmax*ncmax ) THEN
663 *
664 * Note: for either following two causes, there are
665 * almost same number of test cases fail the test.
666 *
667  CALL dlakf2( mm, mplusn-mm, ai, lda,
668  \$ ai( mm+1, mm+1 ), bi,
669  \$ bi( mm+1, mm+1 ), c, ldc )
670 *
671  CALL dgesvd( 'N', 'N', mn2, mn2, c, ldc, s, work,
672  \$ 1, work( 2 ), 1, work( 3 ), lwork-2,
673  \$ info )
674  diftru = s( mn2 )
675 *
676  IF( difest( 2 ).EQ.zero ) THEN
677  IF( diftru.GT.abnrm*ulp )
678  \$ result( 8 ) = ulpinv
679  ELSE IF( diftru.EQ.zero ) THEN
680  IF( difest( 2 ).GT.abnrm*ulp )
681  \$ result( 8 ) = ulpinv
682  ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
683  \$ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
684  result( 8 ) = max( diftru / difest( 2 ),
685  \$ difest( 2 ) / diftru )
686  END IF
687  ntest = ntest + 1
688  END IF
689 *
690 * Test (9)
691 *
692  result( 9 ) = zero
693  IF( linfo.EQ.( mplusn+2 ) ) THEN
694  IF( diftru.GT.abnrm*ulp )
695  \$ result( 9 ) = ulpinv
696  IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
697  \$ result( 9 ) = ulpinv
698  IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
699  \$ result( 9 ) = ulpinv
700  ntest = ntest + 1
701  END IF
702 *
703  ntestt = ntestt + ntest
704 *
705 * Print out tests which fail.
706 *
707  DO 20 j = 1, 9
708  IF( result( j ).GE.thresh ) THEN
709 *
710 * If this is the first test to fail,
711 * print a header to the data file.
712 *
713  IF( nerrs.EQ.0 ) THEN
714  WRITE( nout, fmt = 9995 )'DGX'
715 *
716 * Matrix types
717 *
718  WRITE( nout, fmt = 9993 )
719 *
720 * Tests performed
721 *
722  WRITE( nout, fmt = 9992 )'orthogonal', '''',
723  \$ 'transpose', ( '''', i = 1, 4 )
724 *
725  END IF
726  nerrs = nerrs + 1
727  IF( result( j ).LT.10000.0d0 ) THEN
728  WRITE( nout, fmt = 9991 )mplusn, prtype,
729  \$ weight, m, j, result( j )
730  ELSE
731  WRITE( nout, fmt = 9990 )mplusn, prtype,
732  \$ weight, m, j, result( j )
733  END IF
734  END IF
735  20 CONTINUE
736 *
737  30 CONTINUE
738  40 CONTINUE
739  50 CONTINUE
740  60 CONTINUE
741 *
742  GO TO 150
743 *
744  70 CONTINUE
745 *
746 * Read in data from file to check accuracy of condition estimation
747 * Read input data until N=0
748 *
749  nptknt = 0
750 *
751  80 CONTINUE
752  READ( nin, fmt = *, end = 140 )mplusn
753  IF( mplusn.EQ.0 )
754  \$ GO TO 140
755  READ( nin, fmt = *, end = 140 )n
756  DO 90 i = 1, mplusn
757  READ( nin, fmt = * )( ai( i, j ), j = 1, mplusn )
758  90 CONTINUE
759  DO 100 i = 1, mplusn
760  READ( nin, fmt = * )( bi( i, j ), j = 1, mplusn )
761  100 CONTINUE
762  READ( nin, fmt = * )pltru, diftru
763 *
764  nptknt = nptknt + 1
765  fs = .true.
766  k = 0
767  m = mplusn - n
768 *
769  CALL dlacpy( 'Full', mplusn, mplusn, ai, lda, a, lda )
770  CALL dlacpy( 'Full', mplusn, mplusn, bi, lda, b, lda )
771 *
772 * Compute the Schur factorization while swaping the
773 * m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
774 *
775  CALL dggesx( 'V', 'V', 'S', dlctsx, 'B', mplusn, ai, lda, bi, lda,
776  \$ mm, alphar, alphai, beta, q, lda, z, lda, pl, difest,
777  \$ work, lwork, iwork, liwork, bwork, linfo )
778 *
779  IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 ) THEN
780  result( 1 ) = ulpinv
781  WRITE( nout, fmt = 9998 )'DGGESX', linfo, mplusn, nptknt
782  GO TO 130
783  END IF
784 *
785 * Compute the norm(A, B)
786 * (should this be norm of (A,B) or (AI,BI)?)
787 *
788  CALL dlacpy( 'Full', mplusn, mplusn, ai, lda, work, mplusn )
789  CALL dlacpy( 'Full', mplusn, mplusn, bi, lda,
790  \$ work( mplusn*mplusn+1 ), mplusn )
791  abnrm = dlange( 'Fro', mplusn, 2*mplusn, work, mplusn, work )
792 *
793 * Do tests (1) to (4)
794 *
795  CALL dget51( 1, mplusn, a, lda, ai, lda, q, lda, z, lda, work,
796  \$ result( 1 ) )
797  CALL dget51( 1, mplusn, b, lda, bi, lda, q, lda, z, lda, work,
798  \$ result( 2 ) )
799  CALL dget51( 3, mplusn, b, lda, bi, lda, q, lda, q, lda, work,
800  \$ result( 3 ) )
801  CALL dget51( 3, mplusn, b, lda, bi, lda, z, lda, z, lda, work,
802  \$ result( 4 ) )
803 *
804 * Do tests (5) and (6): check Schur form of A and compare
805 * eigenvalues with diagonals.
806 *
807  ntest = 6
808  temp1 = zero
809  result( 5 ) = zero
810  result( 6 ) = zero
811 *
812  DO 110 j = 1, mplusn
813  ilabad = .false.
814  IF( alphai( j ).EQ.zero ) THEN
815  temp2 = ( abs( alphar( j )-ai( j, j ) ) /
816  \$ max( smlnum, abs( alphar( j ) ), abs( ai( j,
817  \$ j ) ) )+abs( beta( j )-bi( j, j ) ) /
818  \$ max( smlnum, abs( beta( j ) ), abs( bi( j, j ) ) ) )
819  \$ / ulp
820  IF( j.LT.mplusn ) THEN
821  IF( ai( j+1, j ).NE.zero ) THEN
822  ilabad = .true.
823  result( 5 ) = ulpinv
824  END IF
825  END IF
826  IF( j.GT.1 ) THEN
827  IF( ai( j, j-1 ).NE.zero ) THEN
828  ilabad = .true.
829  result( 5 ) = ulpinv
830  END IF
831  END IF
832  ELSE
833  IF( alphai( j ).GT.zero ) THEN
834  i1 = j
835  ELSE
836  i1 = j - 1
837  END IF
838  IF( i1.LE.0 .OR. i1.GE.mplusn ) THEN
839  ilabad = .true.
840  ELSE IF( i1.LT.mplusn-1 ) THEN
841  IF( ai( i1+2, i1+1 ).NE.zero ) THEN
842  ilabad = .true.
843  result( 5 ) = ulpinv
844  END IF
845  ELSE IF( i1.GT.1 ) THEN
846  IF( ai( i1, i1-1 ).NE.zero ) THEN
847  ilabad = .true.
848  result( 5 ) = ulpinv
849  END IF
850  END IF
851  IF( .NOT.ilabad ) THEN
852  CALL dget53( ai( i1, i1 ), lda, bi( i1, i1 ), lda,
853  \$ beta( j ), alphar( j ), alphai( j ), temp2,
854  \$ iinfo )
855  IF( iinfo.GE.3 ) THEN
856  WRITE( nout, fmt = 9997 )iinfo, j, mplusn, nptknt
857  info = abs( iinfo )
858  END IF
859  ELSE
860  temp2 = ulpinv
861  END IF
862  END IF
863  temp1 = max( temp1, temp2 )
864  IF( ilabad ) THEN
865  WRITE( nout, fmt = 9996 )j, mplusn, nptknt
866  END IF
867  110 CONTINUE
868  result( 6 ) = temp1
869 *
870 * Test (7) (if sorting worked) <--------- need to be checked.
871 *
872  ntest = 7
873  result( 7 ) = zero
874  IF( linfo.EQ.mplusn+3 )
875  \$ result( 7 ) = ulpinv
876 *
877 * Test (8): compare the estimated value of DIF and its true value.
878 *
879  ntest = 8
880  result( 8 ) = zero
881  IF( difest( 2 ).EQ.zero ) THEN
882  IF( diftru.GT.abnrm*ulp )
883  \$ result( 8 ) = ulpinv
884  ELSE IF( diftru.EQ.zero ) THEN
885  IF( difest( 2 ).GT.abnrm*ulp )
886  \$ result( 8 ) = ulpinv
887  ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
888  \$ ( diftru*thrsh2.LT.difest( 2 ) ) ) THEN
889  result( 8 ) = max( diftru / difest( 2 ), difest( 2 ) / diftru )
890  END IF
891 *
892 * Test (9)
893 *
894  ntest = 9
895  result( 9 ) = zero
896  IF( linfo.EQ.( mplusn+2 ) ) THEN
897  IF( diftru.GT.abnrm*ulp )
898  \$ result( 9 ) = ulpinv
899  IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
900  \$ result( 9 ) = ulpinv
901  IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
902  \$ result( 9 ) = ulpinv
903  END IF
904 *
905 * Test (10): compare the estimated value of PL and it true value.
906 *
907  ntest = 10
908  result( 10 ) = zero
909  IF( pl( 1 ).EQ.zero ) THEN
910  IF( pltru.GT.abnrm*ulp )
911  \$ result( 10 ) = ulpinv
912  ELSE IF( pltru.EQ.zero ) THEN
913  IF( pl( 1 ).GT.abnrm*ulp )
914  \$ result( 10 ) = ulpinv
915  ELSE IF( ( pltru.GT.thresh*pl( 1 ) ) .OR.
916  \$ ( pltru*thresh.LT.pl( 1 ) ) ) THEN
917  result( 10 ) = ulpinv
918  END IF
919 *
920  ntestt = ntestt + ntest
921 *
922 * Print out tests which fail.
923 *
924  DO 120 j = 1, ntest
925  IF( result( j ).GE.thresh ) THEN
926 *
927 * If this is the first test to fail,
928 * print a header to the data file.
929 *
930  IF( nerrs.EQ.0 ) THEN
931  WRITE( nout, fmt = 9995 )'DGX'
932 *
933 * Matrix types
934 *
935  WRITE( nout, fmt = 9994 )
936 *
937 * Tests performed
938 *
939  WRITE( nout, fmt = 9992 )'orthogonal', '''',
940  \$ 'transpose', ( '''', i = 1, 4 )
941 *
942  END IF
943  nerrs = nerrs + 1
944  IF( result( j ).LT.10000.0d0 ) THEN
945  WRITE( nout, fmt = 9989 )nptknt, mplusn, j, result( j )
946  ELSE
947  WRITE( nout, fmt = 9988 )nptknt, mplusn, j, result( j )
948  END IF
949  END IF
950 *
951  120 CONTINUE
952 *
953  130 CONTINUE
954  GO TO 80
955  140 CONTINUE
956 *
957  150 CONTINUE
958 *
959 * Summary
960 *
961  CALL alasvm( 'DGX', nout, nerrs, ntestt, 0 )
962 *
963  work( 1 ) = maxwrk
964 *
965  RETURN
966 *
967  9999 FORMAT( ' DDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
968  \$ i6, ', JTYPE=', i6, ')' )
969 *
970  9998 FORMAT( ' DDRGSX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
971  \$ i6, ', Input Example #', i2, ')' )
972 *
973  9997 FORMAT( ' DDRGSX: DGET53 returned INFO=', i1, ' for eigenvalue ',
974  \$ i6, '.', / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
975 *
976  9996 FORMAT( ' DDRGSX: S not in Schur form at eigenvalue ', i6, '.',
977  \$ / 9x, 'N=', i6, ', JTYPE=', i6, ')' )
978 *
979  9995 FORMAT( / 1x, a3, ' -- Real Expert Generalized Schur form',
980  \$ ' problem driver' )
981 *
982  9994 FORMAT( 'Input Example' )
983 *
984  9993 FORMAT( ' Matrix types: ', /
985  \$ ' 1: A is a block diagonal matrix of Jordan blocks ',
986  \$ 'and B is the identity ', / ' matrix, ',
987  \$ / ' 2: A and B are upper triangular matrices, ',
988  \$ / ' 3: A and B are as type 2, but each second diagonal ',
989  \$ 'block in A_11 and ', /
990  \$ ' each third diaongal block in A_22 are 2x2 blocks,',
991  \$ / ' 4: A and B are block diagonal matrices, ',
992  \$ / ' 5: (A,B) has potentially close or common ',
993  \$ 'eigenvalues.', / )
994 *
995  9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
996  \$ 'Q and Z are ', a, ',', / 19x,
997  \$ ' a is alpha, b is beta, and ', a, ' means ', a, '.)',
998  \$ / ' 1 = | A - Q S Z', a,
999  \$ ' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
1000  \$ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', a,
1001  \$ ' | / ( n ulp ) 4 = | I - ZZ', a,
1002  \$ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ',
1003  \$ 'Schur form S', / ' 6 = difference between (alpha,beta)',
1004  \$ ' and diagonals of (S,T)', /
1005  \$ ' 7 = 1/ULP if SDIM is not the correct number of ',
1006  \$ 'selected eigenvalues', /
1007  \$ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
1008  \$ 'DIFTRU/DIFEST > 10*THRESH',
1009  \$ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
1010  \$ 'when reordering fails', /
1011  \$ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
1012  \$ 'PLTRU/PLEST > THRESH', /
1013  \$ ' ( Test 10 is only for input examples )', / )
1014  9991 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', d10.3,
1015  \$ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, f8.2 )
1016  9990 FORMAT( ' Matrix order=', i2, ', type=', i2, ', a=', d10.3,
1017  \$ ', order(A_11)=', i2, ', result ', i2, ' is ', 0p, d10.3 )
1018  9989 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
1019  \$ ' result ', i2, ' is', 0p, f8.2 )
1020  9988 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
1021  \$ ' result ', i2, ' is', 1p, d10.3 )
1022 *
1023 * End of DDRGSX
1024 *
1025  END
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:75
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:112
subroutine dlakf2(M, N, A, LDA, B, D, E, Z, LDZ)
DLAKF2
Definition: dlakf2.f:107
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
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:270
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine dget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RESULT)
DGET51
Definition: dget51.f:151
subroutine dget53(A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO)
DGET53
Definition: dget53.f:128
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:367
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:213
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:361