LAPACK  3.8.0 LAPACK: Linear Algebra PACKage
ddrvev.f
Go to the documentation of this file.
1 *> \brief \b DDRVEV
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 DDRVEV( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
12 * NOUNIT, A, LDA, H, WR, WI, WR1, WI1, VL, LDVL,
13 * VR, LDVR, LRE, LDLRE, RESULT, WORK, NWORK,
14 * IWORK, INFO )
15 *
16 * .. Scalar Arguments ..
17 * INTEGER INFO, LDA, LDLRE, LDVL, LDVR, NOUNIT, NSIZES,
18 * \$ NTYPES, NWORK
19 * DOUBLE PRECISION THRESH
20 * ..
21 * .. Array Arguments ..
22 * LOGICAL DOTYPE( * )
23 * INTEGER ISEED( 4 ), IWORK( * ), NN( * )
24 * DOUBLE PRECISION A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ),
25 * \$ RESULT( 7 ), VL( LDVL, * ), VR( LDVR, * ),
26 * \$ WI( * ), WI1( * ), WORK( * ), WR( * ), WR1( * )
27 * ..
28 *
29 *
30 *> \par Purpose:
31 * =============
32 *>
33 *> \verbatim
34 *>
35 *> DDRVEV checks the nonsymmetric eigenvalue problem driver DGEEV.
36 *>
37 *> When DDRVEV is called, a number of matrix "sizes" ("n's") and a
38 *> number of matrix "types" are specified. For each size ("n")
39 *> and each type of matrix, one matrix will be generated and used
40 *> to test the nonsymmetric eigenroutines. For each matrix, 7
41 *> tests will be performed:
42 *>
43 *> (1) | A * VR - VR * W | / ( n |A| ulp )
44 *>
45 *> Here VR is the matrix of unit right eigenvectors.
46 *> W is a block diagonal matrix, with a 1x1 block for each
47 *> real eigenvalue and a 2x2 block for each complex conjugate
48 *> pair. If eigenvalues j and j+1 are a complex conjugate pair,
49 *> so WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the
50 *> 2 x 2 block corresponding to the pair will be:
51 *>
52 *> ( wr wi )
53 *> ( -wi wr )
54 *>
55 *> Such a block multiplying an n x 2 matrix ( ur ui ) on the
56 *> right will be the same as multiplying ur + i*ui by wr + i*wi.
57 *>
58 *> (2) | A**H * VL - VL * W**H | / ( n |A| ulp )
59 *>
60 *> Here VL is the matrix of unit left eigenvectors, A**H is the
61 *> conjugate transpose of A, and W is as above.
62 *>
63 *> (3) | |VR(i)| - 1 | / ulp and whether largest component real
64 *>
65 *> VR(i) denotes the i-th column of VR.
66 *>
67 *> (4) | |VL(i)| - 1 | / ulp and whether largest component real
68 *>
69 *> VL(i) denotes the i-th column of VL.
70 *>
71 *> (5) W(full) = W(partial)
72 *>
73 *> W(full) denotes the eigenvalues computed when both VR and VL
74 *> are also computed, and W(partial) denotes the eigenvalues
75 *> computed when only W, only W and VR, or only W and VL are
76 *> computed.
77 *>
78 *> (6) VR(full) = VR(partial)
79 *>
80 *> VR(full) denotes the right eigenvectors computed when both VR
81 *> and VL are computed, and VR(partial) denotes the result
82 *> when only VR is computed.
83 *>
84 *> (7) VL(full) = VL(partial)
85 *>
86 *> VL(full) denotes the left eigenvectors computed when both VR
87 *> and VL are also computed, and VL(partial) denotes the result
88 *> when only VL is computed.
89 *>
90 *> The "sizes" are specified by an array NN(1:NSIZES); the value of
91 *> each element NN(j) specifies one size.
92 *> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
93 *> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
94 *> Currently, the list of possible types is:
95 *>
96 *> (1) The zero matrix.
97 *> (2) The identity matrix.
98 *> (3) A (transposed) Jordan block, with 1's on the diagonal.
99 *>
100 *> (4) A diagonal matrix with evenly spaced entries
101 *> 1, ..., ULP and random signs.
102 *> (ULP = (first number larger than 1) - 1 )
103 *> (5) A diagonal matrix with geometrically spaced entries
104 *> 1, ..., ULP and random signs.
105 *> (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
106 *> and random signs.
107 *>
108 *> (7) Same as (4), but multiplied by a constant near
109 *> the overflow threshold
110 *> (8) Same as (4), but multiplied by a constant near
111 *> the underflow threshold
112 *>
113 *> (9) A matrix of the form U' T U, where U is orthogonal and
114 *> T has evenly spaced entries 1, ..., ULP with random signs
115 *> on the diagonal and random O(1) entries in the upper
116 *> triangle.
117 *>
118 *> (10) A matrix of the form U' T U, where U is orthogonal and
119 *> T has geometrically spaced entries 1, ..., ULP with random
120 *> signs on the diagonal and random O(1) entries in the upper
121 *> triangle.
122 *>
123 *> (11) A matrix of the form U' T U, where U is orthogonal and
124 *> T has "clustered" entries 1, ULP,..., ULP with random
125 *> signs on the diagonal and random O(1) entries in the upper
126 *> triangle.
127 *>
128 *> (12) A matrix of the form U' T U, where U is orthogonal and
129 *> T has real or complex conjugate paired eigenvalues randomly
130 *> chosen from ( ULP, 1 ) and random O(1) entries in the upper
131 *> triangle.
132 *>
133 *> (13) A matrix of the form X' T X, where X has condition
134 *> SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP
135 *> with random signs on the diagonal and random O(1) entries
136 *> in the upper triangle.
137 *>
138 *> (14) A matrix of the form X' T X, where X has condition
139 *> SQRT( ULP ) and T has geometrically spaced entries
140 *> 1, ..., ULP with random signs on the diagonal and random
141 *> O(1) entries in the upper triangle.
142 *>
143 *> (15) A matrix of the form X' T X, where X has condition
144 *> SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP
145 *> with random signs on the diagonal and random O(1) entries
146 *> in the upper triangle.
147 *>
148 *> (16) A matrix of the form X' T X, where X has condition
149 *> SQRT( ULP ) and T has real or complex conjugate paired
150 *> eigenvalues randomly chosen from ( ULP, 1 ) and random
151 *> O(1) entries in the upper triangle.
152 *>
153 *> (17) Same as (16), but multiplied by a constant
154 *> near the overflow threshold
155 *> (18) Same as (16), but multiplied by a constant
156 *> near the underflow threshold
157 *>
158 *> (19) Nonsymmetric matrix with random entries chosen from (-1,1).
159 *> If N is at least 4, all entries in first two rows and last
160 *> row, and first column and last two columns are zero.
161 *> (20) Same as (19), but multiplied by a constant
162 *> near the overflow threshold
163 *> (21) Same as (19), but multiplied by a constant
164 *> near the underflow threshold
165 *> \endverbatim
166 *
167 * Arguments:
168 * ==========
169 *
170 *> \param[in] NSIZES
171 *> \verbatim
172 *> NSIZES is INTEGER
173 *> The number of sizes of matrices to use. If it is zero,
174 *> DDRVEV does nothing. It must be at least zero.
175 *> \endverbatim
176 *>
177 *> \param[in] NN
178 *> \verbatim
179 *> NN is INTEGER array, dimension (NSIZES)
180 *> An array containing the sizes to be used for the matrices.
181 *> Zero values will be skipped. The values must be at least
182 *> zero.
183 *> \endverbatim
184 *>
185 *> \param[in] NTYPES
186 *> \verbatim
187 *> NTYPES is INTEGER
188 *> The number of elements in DOTYPE. If it is zero, DDRVEV
189 *> does nothing. It must be at least zero. If it is MAXTYP+1
190 *> and NSIZES is 1, then an additional type, MAXTYP+1 is
191 *> defined, which is to use whatever matrix is in A. This
192 *> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
193 *> DOTYPE(MAXTYP+1) is .TRUE. .
194 *> \endverbatim
195 *>
196 *> \param[in] DOTYPE
197 *> \verbatim
198 *> DOTYPE is LOGICAL array, dimension (NTYPES)
199 *> If DOTYPE(j) is .TRUE., then for each size in NN a
200 *> matrix of that size and of type j will be generated.
201 *> If NTYPES is smaller than the maximum number of types
202 *> defined (PARAMETER MAXTYP), then types NTYPES+1 through
203 *> MAXTYP will not be generated. If NTYPES is larger
204 *> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
205 *> will be ignored.
206 *> \endverbatim
207 *>
208 *> \param[in,out] ISEED
209 *> \verbatim
210 *> ISEED is INTEGER array, dimension (4)
211 *> On entry ISEED specifies the seed of the random number
212 *> generator. The array elements should be between 0 and 4095;
213 *> if not they will be reduced mod 4096. Also, ISEED(4) must
214 *> be odd. The random number generator uses a linear
215 *> congruential sequence limited to small integers, and so
216 *> should produce machine independent random numbers. The
217 *> values of ISEED are changed on exit, and can be used in the
218 *> next call to DDRVEV to continue the same random number
219 *> sequence.
220 *> \endverbatim
221 *>
222 *> \param[in] THRESH
223 *> \verbatim
224 *> THRESH is DOUBLE PRECISION
225 *> A test will count as "failed" if the "error", computed as
226 *> described above, exceeds THRESH. Note that the error
227 *> is scaled to be O(1), so THRESH should be a reasonably
228 *> small multiple of 1, e.g., 10 or 100. In particular,
229 *> it should not depend on the precision (single vs. double)
230 *> or the size of the matrix. It must be at least zero.
231 *> \endverbatim
232 *>
233 *> \param[in] NOUNIT
234 *> \verbatim
235 *> NOUNIT is INTEGER
236 *> The FORTRAN unit number for printing out error messages
237 *> (e.g., if a routine returns INFO not equal to 0.)
238 *> \endverbatim
239 *>
240 *> \param[out] A
241 *> \verbatim
242 *> A is DOUBLE PRECISION array, dimension (LDA, max(NN))
243 *> Used to hold the matrix whose eigenvalues are to be
244 *> computed. On exit, A contains the last matrix actually used.
245 *> \endverbatim
246 *>
247 *> \param[in] LDA
248 *> \verbatim
249 *> LDA is INTEGER
250 *> The leading dimension of A, and H. LDA must be at
251 *> least 1 and at least max(NN).
252 *> \endverbatim
253 *>
254 *> \param[out] H
255 *> \verbatim
256 *> H is DOUBLE PRECISION array, dimension (LDA, max(NN))
257 *> Another copy of the test matrix A, modified by DGEEV.
258 *> \endverbatim
259 *>
260 *> \param[out] WR
261 *> \verbatim
262 *> WR is DOUBLE PRECISION array, dimension (max(NN))
263 *> \endverbatim
264 *>
265 *> \param[out] WI
266 *> \verbatim
267 *> WI is DOUBLE PRECISION array, dimension (max(NN))
268 *>
269 *> The real and imaginary parts of the eigenvalues of A.
270 *> On exit, WR + WI*i are the eigenvalues of the matrix in A.
271 *> \endverbatim
272 *>
273 *> \param[out] WR1
274 *> \verbatim
275 *> WR1 is DOUBLE PRECISION array, dimension (max(NN))
276 *> \endverbatim
277 *>
278 *> \param[out] WI1
279 *> \verbatim
280 *> WI1 is DOUBLE PRECISION array, dimension (max(NN))
281 *>
282 *> Like WR, WI, these arrays contain the eigenvalues of A,
283 *> but those computed when DGEEV only computes a partial
284 *> eigendecomposition, i.e. not the eigenvalues and left
285 *> and right eigenvectors.
286 *> \endverbatim
287 *>
288 *> \param[out] VL
289 *> \verbatim
290 *> VL is DOUBLE PRECISION array, dimension (LDVL, max(NN))
291 *> VL holds the computed left eigenvectors.
292 *> \endverbatim
293 *>
294 *> \param[in] LDVL
295 *> \verbatim
296 *> LDVL is INTEGER
297 *> Leading dimension of VL. Must be at least max(1,max(NN)).
298 *> \endverbatim
299 *>
300 *> \param[out] VR
301 *> \verbatim
302 *> VR is DOUBLE PRECISION array, dimension (LDVR, max(NN))
303 *> VR holds the computed right eigenvectors.
304 *> \endverbatim
305 *>
306 *> \param[in] LDVR
307 *> \verbatim
308 *> LDVR is INTEGER
309 *> Leading dimension of VR. Must be at least max(1,max(NN)).
310 *> \endverbatim
311 *>
312 *> \param[out] LRE
313 *> \verbatim
314 *> LRE is DOUBLE PRECISION array, dimension (LDLRE,max(NN))
315 *> LRE holds the computed right or left eigenvectors.
316 *> \endverbatim
317 *>
318 *> \param[in] LDLRE
319 *> \verbatim
320 *> LDLRE is INTEGER
321 *> Leading dimension of LRE. Must be at least max(1,max(NN)).
322 *> \endverbatim
323 *>
324 *> \param[out] RESULT
325 *> \verbatim
326 *> RESULT is DOUBLE PRECISION array, dimension (7)
327 *> The values computed by the seven tests described above.
328 *> The values are currently limited to 1/ulp, to avoid overflow.
329 *> \endverbatim
330 *>
331 *> \param[out] WORK
332 *> \verbatim
333 *> WORK is DOUBLE PRECISION array, dimension (NWORK)
334 *> \endverbatim
335 *>
336 *> \param[in] NWORK
337 *> \verbatim
338 *> NWORK is INTEGER
339 *> The number of entries in WORK. This must be at least
340 *> 5*NN(j)+2*NN(j)**2 for all j.
341 *> \endverbatim
342 *>
343 *> \param[out] IWORK
344 *> \verbatim
345 *> IWORK is INTEGER array, dimension (max(NN))
346 *> \endverbatim
347 *>
348 *> \param[out] INFO
349 *> \verbatim
350 *> INFO is INTEGER
351 *> If 0, then everything ran OK.
352 *> -1: NSIZES < 0
353 *> -2: Some NN(j) < 0
354 *> -3: NTYPES < 0
355 *> -6: THRESH < 0
356 *> -9: LDA < 1 or LDA < NMAX, where NMAX is max( NN(j) ).
357 *> -16: LDVL < 1 or LDVL < NMAX, where NMAX is max( NN(j) ).
358 *> -18: LDVR < 1 or LDVR < NMAX, where NMAX is max( NN(j) ).
359 *> -20: LDLRE < 1 or LDLRE < NMAX, where NMAX is max( NN(j) ).
360 *> -23: NWORK too small.
361 *> If DLATMR, SLATMS, SLATME or DGEEV returns an error code,
362 *> the absolute value of it is returned.
363 *>
364 *>-----------------------------------------------------------------------
365 *>
366 *> Some Local Variables and Parameters:
367 *> ---- ----- --------- --- ----------
368 *>
369 *> ZERO, ONE Real 0 and 1.
370 *> MAXTYP The number of types defined.
371 *> NMAX Largest value in NN.
372 *> NERRS The number of tests which have exceeded THRESH
373 *> COND, CONDS,
374 *> IMODE Values to be passed to the matrix generators.
375 *> ANORM Norm of A; passed to matrix generators.
376 *>
377 *> OVFL, UNFL Overflow and underflow thresholds.
378 *> ULP, ULPINV Finest relative precision and its inverse.
379 *> RTULP, RTULPI Square roots of the previous 4 values.
380 *>
381 *> The following four arrays decode JTYPE:
382 *> KTYPE(j) The general type (1-10) for type "j".
383 *> KMODE(j) The MODE value to be passed to the matrix
384 *> generator for type "j".
385 *> KMAGN(j) The order of magnitude ( O(1),
386 *> O(overflow^(1/2) ), O(underflow^(1/2) )
387 *> KCONDS(j) Selectw whether CONDS is to be 1 or
388 *> 1/sqrt(ulp). (0 means irrelevant.)
389 *> \endverbatim
390 *
391 * Authors:
392 * ========
393 *
394 *> \author Univ. of Tennessee
395 *> \author Univ. of California Berkeley
396 *> \author Univ. of Colorado Denver
397 *> \author NAG Ltd.
398 *
399 *> \date December 2016
400 *
401 *> \ingroup double_eig
402 *
403 * =====================================================================
404  SUBROUTINE ddrvev( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
405  \$ NOUNIT, A, LDA, H, WR, WI, WR1, WI1, VL, LDVL,
406  \$ VR, LDVR, LRE, LDLRE, RESULT, WORK, NWORK,
407  \$ IWORK, INFO )
408 *
409 * -- LAPACK test routine (version 3.7.0) --
410 * -- LAPACK is a software package provided by Univ. of Tennessee, --
411 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
412 * December 2016
413 *
414 * .. Scalar Arguments ..
415  INTEGER INFO, LDA, LDLRE, LDVL, LDVR, NOUNIT, NSIZES,
416  \$ ntypes, nwork
417  DOUBLE PRECISION THRESH
418 * ..
419 * .. Array Arguments ..
420  LOGICAL DOTYPE( * )
421  INTEGER ISEED( 4 ), IWORK( * ), NN( * )
422  DOUBLE PRECISION A( lda, * ), H( lda, * ), LRE( ldlre, * ),
423  \$ result( 7 ), vl( ldvl, * ), vr( ldvr, * ),
424  \$ wi( * ), wi1( * ), work( * ), wr( * ), wr1( * )
425 * ..
426 *
427 * =====================================================================
428 *
429 * .. Parameters ..
430  DOUBLE PRECISION ZERO, ONE
431  parameter( zero = 0.0d0, one = 1.0d0 )
432  DOUBLE PRECISION TWO
433  parameter( two = 2.0d0 )
434  INTEGER MAXTYP
435  parameter( maxtyp = 21 )
436 * ..
437 * .. Local Scalars ..
439  CHARACTER*3 PATH
440  INTEGER IINFO, IMODE, ITYPE, IWK, J, JCOL, JJ, JSIZE,
441  \$ jtype, mtypes, n, nerrs, nfail, nmax, nnwork,
442  \$ ntest, ntestf, ntestt
443  DOUBLE PRECISION ANORM, COND, CONDS, OVFL, RTULP, RTULPI, TNRM,
444  \$ ulp, ulpinv, unfl, vmx, vrmx, vtst
445 * ..
446 * .. Local Arrays ..
448  INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( maxtyp ),
449  \$ kmagn( maxtyp ), kmode( maxtyp ),
450  \$ ktype( maxtyp )
451  DOUBLE PRECISION DUM( 1 ), RES( 2 )
452 * ..
453 * .. External Functions ..
454  DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2
455  EXTERNAL dlamch, dlapy2, dnrm2
456 * ..
457 * .. External Subroutines ..
458  EXTERNAL dgeev, dget22, dlabad, dlacpy, dlaset, dlasum,
460 * ..
461 * .. Intrinsic Functions ..
462  INTRINSIC abs, max, min, sqrt
463 * ..
464 * .. Data statements ..
465  DATA ktype / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
466  DATA kmagn / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
467  \$ 3, 1, 2, 3 /
468  DATA kmode / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
469  \$ 1, 5, 5, 5, 4, 3, 1 /
470  DATA kconds / 3*0, 5*0, 4*1, 6*2, 3*0 /
471 * ..
472 * .. Executable Statements ..
473 *
474  path( 1: 1 ) = 'Double precision'
475  path( 2: 3 ) = 'EV'
476 *
477 * Check for errors
478 *
479  ntestt = 0
480  ntestf = 0
481  info = 0
482 *
483 * Important constants
484 *
486  nmax = 0
487  DO 10 j = 1, nsizes
488  nmax = max( nmax, nn( j ) )
489  IF( nn( j ).LT.0 )
491  10 CONTINUE
492 *
493 * Check for errors
494 *
495  IF( nsizes.LT.0 ) THEN
496  info = -1
497  ELSE IF( badnn ) THEN
498  info = -2
499  ELSE IF( ntypes.LT.0 ) THEN
500  info = -3
501  ELSE IF( thresh.LT.zero ) THEN
502  info = -6
503  ELSE IF( nounit.LE.0 ) THEN
504  info = -7
505  ELSE IF( lda.LT.1 .OR. lda.LT.nmax ) THEN
506  info = -9
507  ELSE IF( ldvl.LT.1 .OR. ldvl.LT.nmax ) THEN
508  info = -16
509  ELSE IF( ldvr.LT.1 .OR. ldvr.LT.nmax ) THEN
510  info = -18
511  ELSE IF( ldlre.LT.1 .OR. ldlre.LT.nmax ) THEN
512  info = -20
513  ELSE IF( 5*nmax+2*nmax**2.GT.nwork ) THEN
514  info = -23
515  END IF
516 *
517  IF( info.NE.0 ) THEN
518  CALL xerbla( 'DDRVEV', -info )
519  RETURN
520  END IF
521 *
522 * Quick return if nothing to do
523 *
524  IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
525  \$ RETURN
526 *
527 * More Important constants
528 *
529  unfl = dlamch( 'Safe minimum' )
530  ovfl = one / unfl
531  CALL dlabad( unfl, ovfl )
532  ulp = dlamch( 'Precision' )
533  ulpinv = one / ulp
534  rtulp = sqrt( ulp )
535  rtulpi = one / rtulp
536 *
537 * Loop over sizes, types
538 *
539  nerrs = 0
540 *
541  DO 270 jsize = 1, nsizes
542  n = nn( jsize )
543  IF( nsizes.NE.1 ) THEN
544  mtypes = min( maxtyp, ntypes )
545  ELSE
546  mtypes = min( maxtyp+1, ntypes )
547  END IF
548 *
549  DO 260 jtype = 1, mtypes
550  IF( .NOT.dotype( jtype ) )
551  \$ GO TO 260
552 *
553 * Save ISEED in case of an error.
554 *
555  DO 20 j = 1, 4
556  ioldsd( j ) = iseed( j )
557  20 CONTINUE
558 *
559 * Compute "A"
560 *
561 * Control parameters:
562 *
563 * KMAGN KCONDS KMODE KTYPE
564 * =1 O(1) 1 clustered 1 zero
565 * =2 large large clustered 2 identity
566 * =3 small exponential Jordan
567 * =4 arithmetic diagonal, (w/ eigenvalues)
568 * =5 random log symmetric, w/ eigenvalues
569 * =6 random general, w/ eigenvalues
570 * =7 random diagonal
571 * =8 random symmetric
572 * =9 random general
573 * =10 random triangular
574 *
575  IF( mtypes.GT.maxtyp )
576  \$ GO TO 90
577 *
578  itype = ktype( jtype )
579  imode = kmode( jtype )
580 *
581 * Compute norm
582 *
583  GO TO ( 30, 40, 50 )kmagn( jtype )
584 *
585  30 CONTINUE
586  anorm = one
587  GO TO 60
588 *
589  40 CONTINUE
590  anorm = ovfl*ulp
591  GO TO 60
592 *
593  50 CONTINUE
594  anorm = unfl*ulpinv
595  GO TO 60
596 *
597  60 CONTINUE
598 *
599  CALL dlaset( 'Full', lda, n, zero, zero, a, lda )
600  iinfo = 0
601  cond = ulpinv
602 *
603 * Special Matrices -- Identity & Jordan block
604 *
605 * Zero
606 *
607  IF( itype.EQ.1 ) THEN
608  iinfo = 0
609 *
610  ELSE IF( itype.EQ.2 ) THEN
611 *
612 * Identity
613 *
614  DO 70 jcol = 1, n
615  a( jcol, jcol ) = anorm
616  70 CONTINUE
617 *
618  ELSE IF( itype.EQ.3 ) THEN
619 *
620 * Jordan Block
621 *
622  DO 80 jcol = 1, n
623  a( jcol, jcol ) = anorm
624  IF( jcol.GT.1 )
625  \$ a( jcol, jcol-1 ) = one
626  80 CONTINUE
627 *
628  ELSE IF( itype.EQ.4 ) THEN
629 *
630 * Diagonal Matrix, [Eigen]values Specified
631 *
632  CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
633  \$ anorm, 0, 0, 'N', a, lda, work( n+1 ),
634  \$ iinfo )
635 *
636  ELSE IF( itype.EQ.5 ) THEN
637 *
638 * Symmetric, eigenvalues specified
639 *
640  CALL dlatms( n, n, 'S', iseed, 'S', work, imode, cond,
641  \$ anorm, n, n, 'N', a, lda, work( n+1 ),
642  \$ iinfo )
643 *
644  ELSE IF( itype.EQ.6 ) THEN
645 *
646 * General, eigenvalues specified
647 *
648  IF( kconds( jtype ).EQ.1 ) THEN
649  conds = one
650  ELSE IF( kconds( jtype ).EQ.2 ) THEN
651  conds = rtulpi
652  ELSE
653  conds = zero
654  END IF
655 *
656  adumma( 1 ) = ' '
657  CALL dlatme( n, 'S', iseed, work, imode, cond, one,
658  \$ adumma, 'T', 'T', 'T', work( n+1 ), 4,
659  \$ conds, n, n, anorm, a, lda, work( 2*n+1 ),
660  \$ iinfo )
661 *
662  ELSE IF( itype.EQ.7 ) THEN
663 *
664 * Diagonal, random eigenvalues
665 *
666  CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
667  \$ 'T', 'N', work( n+1 ), 1, one,
668  \$ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
669  \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
670 *
671  ELSE IF( itype.EQ.8 ) THEN
672 *
673 * Symmetric, random eigenvalues
674 *
675  CALL dlatmr( n, n, 'S', iseed, 'S', work, 6, one, one,
676  \$ 'T', 'N', work( n+1 ), 1, one,
677  \$ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
678  \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
679 *
680  ELSE IF( itype.EQ.9 ) THEN
681 *
682 * General, random eigenvalues
683 *
684  CALL dlatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
685  \$ 'T', 'N', work( n+1 ), 1, one,
686  \$ work( 2*n+1 ), 1, one, 'N', idumma, n, n,
687  \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
688  IF( n.GE.4 ) THEN
689  CALL dlaset( 'Full', 2, n, zero, zero, a, lda )
690  CALL dlaset( 'Full', n-3, 1, zero, zero, a( 3, 1 ),
691  \$ lda )
692  CALL dlaset( 'Full', n-3, 2, zero, zero, a( 3, n-1 ),
693  \$ lda )
694  CALL dlaset( 'Full', 1, n, zero, zero, a( n, 1 ),
695  \$ lda )
696  END IF
697 *
698  ELSE IF( itype.EQ.10 ) THEN
699 *
700 * Triangular, random eigenvalues
701 *
702  CALL dlatmr( n, n, 'S', iseed, 'N', work, 6, one, one,
703  \$ 'T', 'N', work( n+1 ), 1, one,
704  \$ work( 2*n+1 ), 1, one, 'N', idumma, n, 0,
705  \$ zero, anorm, 'NO', a, lda, iwork, iinfo )
706 *
707  ELSE
708 *
709  iinfo = 1
710  END IF
711 *
712  IF( iinfo.NE.0 ) THEN
713  WRITE( nounit, fmt = 9993 )'Generator', iinfo, n, jtype,
714  \$ ioldsd
715  info = abs( iinfo )
716  RETURN
717  END IF
718 *
719  90 CONTINUE
720 *
721 * Test for minimal and generous workspace
722 *
723  DO 250 iwk = 1, 2
724  IF( iwk.EQ.1 ) THEN
725  nnwork = 4*n
726  ELSE
727  nnwork = 5*n + 2*n**2
728  END IF
729  nnwork = max( nnwork, 1 )
730 *
731 * Initialize RESULT
732 *
733  DO 100 j = 1, 7
734  result( j ) = -one
735  100 CONTINUE
736 *
737 * Compute eigenvalues and eigenvectors, and test them
738 *
739  CALL dlacpy( 'F', n, n, a, lda, h, lda )
740  CALL dgeev( 'V', 'V', n, h, lda, wr, wi, vl, ldvl, vr,
741  \$ ldvr, work, nnwork, iinfo )
742  IF( iinfo.NE.0 ) THEN
743  result( 1 ) = ulpinv
744  WRITE( nounit, fmt = 9993 )'DGEEV1', iinfo, n, jtype,
745  \$ ioldsd
746  info = abs( iinfo )
747  GO TO 220
748  END IF
749 *
750 * Do Test (1)
751 *
752  CALL dget22( 'N', 'N', 'N', n, a, lda, vr, ldvr, wr, wi,
753  \$ work, res )
754  result( 1 ) = res( 1 )
755 *
756 * Do Test (2)
757 *
758  CALL dget22( 'T', 'N', 'T', n, a, lda, vl, ldvl, wr, wi,
759  \$ work, res )
760  result( 2 ) = res( 1 )
761 *
762 * Do Test (3)
763 *
764  DO 120 j = 1, n
765  tnrm = one
766  IF( wi( j ).EQ.zero ) THEN
767  tnrm = dnrm2( n, vr( 1, j ), 1 )
768  ELSE IF( wi( j ).GT.zero ) THEN
769  tnrm = dlapy2( dnrm2( n, vr( 1, j ), 1 ),
770  \$ dnrm2( n, vr( 1, j+1 ), 1 ) )
771  END IF
772  result( 3 ) = max( result( 3 ),
773  \$ min( ulpinv, abs( tnrm-one ) / ulp ) )
774  IF( wi( j ).GT.zero ) THEN
775  vmx = zero
776  vrmx = zero
777  DO 110 jj = 1, n
778  vtst = dlapy2( vr( jj, j ), vr( jj, j+1 ) )
779  IF( vtst.GT.vmx )
780  \$ vmx = vtst
781  IF( vr( jj, j+1 ).EQ.zero .AND.
782  \$ abs( vr( jj, j ) ).GT.vrmx )
783  \$ vrmx = abs( vr( jj, j ) )
784  110 CONTINUE
785  IF( vrmx / vmx.LT.one-two*ulp )
786  \$ result( 3 ) = ulpinv
787  END IF
788  120 CONTINUE
789 *
790 * Do Test (4)
791 *
792  DO 140 j = 1, n
793  tnrm = one
794  IF( wi( j ).EQ.zero ) THEN
795  tnrm = dnrm2( n, vl( 1, j ), 1 )
796  ELSE IF( wi( j ).GT.zero ) THEN
797  tnrm = dlapy2( dnrm2( n, vl( 1, j ), 1 ),
798  \$ dnrm2( n, vl( 1, j+1 ), 1 ) )
799  END IF
800  result( 4 ) = max( result( 4 ),
801  \$ min( ulpinv, abs( tnrm-one ) / ulp ) )
802  IF( wi( j ).GT.zero ) THEN
803  vmx = zero
804  vrmx = zero
805  DO 130 jj = 1, n
806  vtst = dlapy2( vl( jj, j ), vl( jj, j+1 ) )
807  IF( vtst.GT.vmx )
808  \$ vmx = vtst
809  IF( vl( jj, j+1 ).EQ.zero .AND.
810  \$ abs( vl( jj, j ) ).GT.vrmx )
811  \$ vrmx = abs( vl( jj, j ) )
812  130 CONTINUE
813  IF( vrmx / vmx.LT.one-two*ulp )
814  \$ result( 4 ) = ulpinv
815  END IF
816  140 CONTINUE
817 *
818 * Compute eigenvalues only, and test them
819 *
820  CALL dlacpy( 'F', n, n, a, lda, h, lda )
821  CALL dgeev( 'N', 'N', n, h, lda, wr1, wi1, dum, 1, dum,
822  \$ 1, work, nnwork, iinfo )
823  IF( iinfo.NE.0 ) THEN
824  result( 1 ) = ulpinv
825  WRITE( nounit, fmt = 9993 )'DGEEV2', iinfo, n, jtype,
826  \$ ioldsd
827  info = abs( iinfo )
828  GO TO 220
829  END IF
830 *
831 * Do Test (5)
832 *
833  DO 150 j = 1, n
834  IF( wr( j ).NE.wr1( j ) .OR. wi( j ).NE.wi1( j ) )
835  \$ result( 5 ) = ulpinv
836  150 CONTINUE
837 *
838 * Compute eigenvalues and right eigenvectors, and test them
839 *
840  CALL dlacpy( 'F', n, n, a, lda, h, lda )
841  CALL dgeev( 'N', 'V', n, h, lda, wr1, wi1, dum, 1, lre,
842  \$ ldlre, work, nnwork, iinfo )
843  IF( iinfo.NE.0 ) THEN
844  result( 1 ) = ulpinv
845  WRITE( nounit, fmt = 9993 )'DGEEV3', iinfo, n, jtype,
846  \$ ioldsd
847  info = abs( iinfo )
848  GO TO 220
849  END IF
850 *
851 * Do Test (5) again
852 *
853  DO 160 j = 1, n
854  IF( wr( j ).NE.wr1( j ) .OR. wi( j ).NE.wi1( j ) )
855  \$ result( 5 ) = ulpinv
856  160 CONTINUE
857 *
858 * Do Test (6)
859 *
860  DO 180 j = 1, n
861  DO 170 jj = 1, n
862  IF( vr( j, jj ).NE.lre( j, jj ) )
863  \$ result( 6 ) = ulpinv
864  170 CONTINUE
865  180 CONTINUE
866 *
867 * Compute eigenvalues and left eigenvectors, and test them
868 *
869  CALL dlacpy( 'F', n, n, a, lda, h, lda )
870  CALL dgeev( 'V', 'N', n, h, lda, wr1, wi1, lre, ldlre,
871  \$ dum, 1, work, nnwork, iinfo )
872  IF( iinfo.NE.0 ) THEN
873  result( 1 ) = ulpinv
874  WRITE( nounit, fmt = 9993 )'DGEEV4', iinfo, n, jtype,
875  \$ ioldsd
876  info = abs( iinfo )
877  GO TO 220
878  END IF
879 *
880 * Do Test (5) again
881 *
882  DO 190 j = 1, n
883  IF( wr( j ).NE.wr1( j ) .OR. wi( j ).NE.wi1( j ) )
884  \$ result( 5 ) = ulpinv
885  190 CONTINUE
886 *
887 * Do Test (7)
888 *
889  DO 210 j = 1, n
890  DO 200 jj = 1, n
891  IF( vl( j, jj ).NE.lre( j, jj ) )
892  \$ result( 7 ) = ulpinv
893  200 CONTINUE
894  210 CONTINUE
895 *
896 * End of Loop -- Check for RESULT(j) > THRESH
897 *
898  220 CONTINUE
899 *
900  ntest = 0
901  nfail = 0
902  DO 230 j = 1, 7
903  IF( result( j ).GE.zero )
904  \$ ntest = ntest + 1
905  IF( result( j ).GE.thresh )
906  \$ nfail = nfail + 1
907  230 CONTINUE
908 *
909  IF( nfail.GT.0 )
910  \$ ntestf = ntestf + 1
911  IF( ntestf.EQ.1 ) THEN
912  WRITE( nounit, fmt = 9999 )path
913  WRITE( nounit, fmt = 9998 )
914  WRITE( nounit, fmt = 9997 )
915  WRITE( nounit, fmt = 9996 )
916  WRITE( nounit, fmt = 9995 )thresh
917  ntestf = 2
918  END IF
919 *
920  DO 240 j = 1, 7
921  IF( result( j ).GE.thresh ) THEN
922  WRITE( nounit, fmt = 9994 )n, iwk, ioldsd, jtype,
923  \$ j, result( j )
924  END IF
925  240 CONTINUE
926 *
927  nerrs = nerrs + nfail
928  ntestt = ntestt + ntest
929 *
930  250 CONTINUE
931  260 CONTINUE
932  270 CONTINUE
933 *
934 * Summary
935 *
936  CALL dlasum( path, nounit, nerrs, ntestt )
937 *
938  9999 FORMAT( / 1x, a3, ' -- Real Eigenvalue-Eigenvector Decomposition',
939  \$ ' Driver', / ' Matrix types (see DDRVEV for details): ' )
940 *
941  9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ',
942  \$ ' ', ' 5=Diagonal: geometr. spaced entries.',
943  \$ / ' 2=Identity matrix. ', ' 6=Diagona',
944  \$ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
945  \$ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
946  \$ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
947  \$ 'mall, evenly spaced.' )
948  9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
949  \$ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
950  \$ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
951  \$ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
952  \$ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
953  \$ 'lex ', / ' 12=Well-cond., random complex ', 6x, ' ',
954  \$ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
955  \$ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
956  \$ ' complx ' )
957  9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
958  \$ 'with small random entries.', / ' 20=Matrix with large ran',
959  \$ 'dom entries. ', / )
960  9995 FORMAT( ' Tests performed with test threshold =', f8.2,
961  \$ / / ' 1 = | A VR - VR W | / ( n |A| ulp ) ',
962  \$ / ' 2 = | transpose(A) VL - VL W | / ( n |A| ulp ) ',
963  \$ / ' 3 = | |VR(i)| - 1 | / ulp ',
964  \$ / ' 4 = | |VL(i)| - 1 | / ulp ',
965  \$ / ' 5 = 0 if W same no matter if VR or VL computed,',
966  \$ ' 1/ulp otherwise', /
967  \$ ' 6 = 0 if VR same no matter if VL computed,',
968  \$ ' 1/ulp otherwise', /
969  \$ ' 7 = 0 if VL same no matter if VR computed,',
970  \$ ' 1/ulp otherwise', / )
971  9994 FORMAT( ' N=', i5, ', IWK=', i2, ', seed=', 4( i4, ',' ),
972  \$ ' type ', i2, ', test(', i2, ')=', g10.3 )
973  9993 FORMAT( ' DDRVEV: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
974  \$ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
975 *
976  RETURN
977 *
978 * End of DDRVEV
979 *
980  END
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 ddrvev(NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, NOUNIT, A, LDA, H, WR, WI, WR1, WI1, VL, LDVL, VR, LDVR, LRE, LDLRE, RESULT, WORK, NWORK, IWORK, INFO)
DDRVEV
Definition: ddrvev.f:408
subroutine dlatme(N, DIST, ISEED, D, MODE, COND, DMAX, EI, RSIGN, UPPER, SIM, DS, MODES, CONDS, KL, KU, ANORM, A, LDA, WORK, INFO)
DLATME
Definition: dlatme.f:334
subroutine dlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
DLATMS
Definition: dlatms.f:323
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 dgeev(JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR, LDVR, WORK, LWORK, INFO)
DGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices ...
Definition: dgeev.f:193
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62