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