LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
zchkhb2stg.f
Go to the documentation of this file.
1 *> \brief \b ZCHKHBSTG
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 ZCHKHBSTG( NSIZES, NN, NWDTHS, KK, NTYPES, DOTYPE,
12 * ISEED, THRESH, NOUNIT, A, LDA, SD, SE, D1,
13 * D2, D3, U, LDU, WORK, LWORK, RWORK RESULT,
14 * INFO )
15 *
16 * .. Scalar Arguments ..
17 * INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES,
18 * $ NWDTHS
19 * DOUBLE PRECISION THRESH
20 * ..
21 * .. Array Arguments ..
22 * LOGICAL DOTYPE( * )
23 * INTEGER ISEED( 4 ), KK( * ), NN( * )
24 * DOUBLE PRECISION RESULT( * ), RWORK( * ), SD( * ), SE( * )
25 * COMPLEX*16 A( LDA, * ), U( LDU, * ), WORK( * )
26 * ..
27 *
28 *
29 *> \par Purpose:
30 * =============
31 *>
32 *> \verbatim
33 *>
34 *> ZCHKHBSTG tests the reduction of a Hermitian band matrix to tridiagonal
35 *> from, used with the Hermitian eigenvalue problem.
36 *>
37 *> ZHBTRD factors a Hermitian band matrix A as U S U* , where * means
38 *> conjugate transpose, S is symmetric tridiagonal, and U is unitary.
39 *> ZHBTRD can use either just the lower or just the upper triangle
40 *> of A; ZCHKHBSTG checks both cases.
41 *>
42 *> ZHETRD_HB2ST factors a Hermitian band matrix A as U S U* ,
43 *> where * means conjugate transpose, S is symmetric tridiagonal, and U is
44 *> unitary. ZHETRD_HB2ST can use either just the lower or just
45 *> the upper triangle of A; ZCHKHBSTG checks both cases.
46 *>
47 *> DSTEQR factors S as Z D1 Z'.
48 *> D1 is the matrix of eigenvalues computed when Z is not computed
49 *> and from the S resulting of DSBTRD "U" (used as reference for DSYTRD_SB2ST)
50 *> D2 is the matrix of eigenvalues computed when Z is not computed
51 *> and from the S resulting of DSYTRD_SB2ST "U".
52 *> D3 is the matrix of eigenvalues computed when Z is not computed
53 *> and from the S resulting of DSYTRD_SB2ST "L".
54 *>
55 *> When ZCHKHBSTG is called, a number of matrix "sizes" ("n's"), a number
56 *> of bandwidths ("k's"), and a number of matrix "types" are
57 *> specified. For each size ("n"), each bandwidth ("k") less than or
58 *> equal to "n", and each type of matrix, one matrix will be generated
59 *> and used to test the hermitian banded reduction routine. For each
60 *> matrix, a number of tests will be performed:
61 *>
62 *> (1) | A - V S V* | / ( |A| n ulp ) computed by ZHBTRD with
63 *> UPLO='U'
64 *>
65 *> (2) | I - UU* | / ( n ulp )
66 *>
67 *> (3) | A - V S V* | / ( |A| n ulp ) computed by ZHBTRD with
68 *> UPLO='L'
69 *>
70 *> (4) | I - UU* | / ( n ulp )
71 *>
72 *> (5) | D1 - D2 | / ( |D1| ulp ) where D1 is computed by
73 *> DSBTRD with UPLO='U' and
74 *> D2 is computed by
75 *> ZHETRD_HB2ST with UPLO='U'
76 *>
77 *> (6) | D1 - D3 | / ( |D1| ulp ) where D1 is computed by
78 *> DSBTRD with UPLO='U' and
79 *> D3 is computed by
80 *> ZHETRD_HB2ST with UPLO='L'
81 *>
82 *> The "sizes" are specified by an array NN(1:NSIZES); the value of
83 *> each element NN(j) specifies one size.
84 *> The "types" are specified by a logical array DOTYPE( 1:NTYPES );
85 *> if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
86 *> Currently, the list of possible types is:
87 *>
88 *> (1) The zero matrix.
89 *> (2) The identity matrix.
90 *>
91 *> (3) A diagonal matrix with evenly spaced entries
92 *> 1, ..., ULP and random signs.
93 *> (ULP = (first number larger than 1) - 1 )
94 *> (4) A diagonal matrix with geometrically spaced entries
95 *> 1, ..., ULP and random signs.
96 *> (5) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
97 *> and random signs.
98 *>
99 *> (6) Same as (4), but multiplied by SQRT( overflow threshold )
100 *> (7) Same as (4), but multiplied by SQRT( underflow threshold )
101 *>
102 *> (8) A matrix of the form U* D U, where U is unitary and
103 *> D has evenly spaced entries 1, ..., ULP with random signs
104 *> on the diagonal.
105 *>
106 *> (9) A matrix of the form U* D U, where U is unitary and
107 *> D has geometrically spaced entries 1, ..., ULP with random
108 *> signs on the diagonal.
109 *>
110 *> (10) A matrix of the form U* D U, where U is unitary and
111 *> D has "clustered" entries 1, ULP,..., ULP with random
112 *> signs on the diagonal.
113 *>
114 *> (11) Same as (8), but multiplied by SQRT( overflow threshold )
115 *> (12) Same as (8), but multiplied by SQRT( underflow threshold )
116 *>
117 *> (13) Hermitian matrix with random entries chosen from (-1,1).
118 *> (14) Same as (13), but multiplied by SQRT( overflow threshold )
119 *> (15) Same as (13), but multiplied by SQRT( underflow threshold )
120 *> \endverbatim
121 *
122 * Arguments:
123 * ==========
124 *
125 *> \param[in] NSIZES
126 *> \verbatim
127 *> NSIZES is INTEGER
128 *> The number of sizes of matrices to use. If it is zero,
129 *> ZCHKHBSTG does nothing. It must be at least zero.
130 *> \endverbatim
131 *>
132 *> \param[in] NN
133 *> \verbatim
134 *> NN is INTEGER array, dimension (NSIZES)
135 *> An array containing the sizes to be used for the matrices.
136 *> Zero values will be skipped. The values must be at least
137 *> zero.
138 *> \endverbatim
139 *>
140 *> \param[in] NWDTHS
141 *> \verbatim
142 *> NWDTHS is INTEGER
143 *> The number of bandwidths to use. If it is zero,
144 *> ZCHKHBSTG does nothing. It must be at least zero.
145 *> \endverbatim
146 *>
147 *> \param[in] KK
148 *> \verbatim
149 *> KK is INTEGER array, dimension (NWDTHS)
150 *> An array containing the bandwidths to be used for the band
151 *> matrices. The values must be at least zero.
152 *> \endverbatim
153 *>
154 *> \param[in] NTYPES
155 *> \verbatim
156 *> NTYPES is INTEGER
157 *> The number of elements in DOTYPE. If it is zero, ZCHKHBSTG
158 *> does nothing. It must be at least zero. If it is MAXTYP+1
159 *> and NSIZES is 1, then an additional type, MAXTYP+1 is
160 *> defined, which is to use whatever matrix is in A. This
161 *> is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
162 *> DOTYPE(MAXTYP+1) is .TRUE. .
163 *> \endverbatim
164 *>
165 *> \param[in] DOTYPE
166 *> \verbatim
167 *> DOTYPE is LOGICAL array, dimension (NTYPES)
168 *> If DOTYPE(j) is .TRUE., then for each size in NN a
169 *> matrix of that size and of type j will be generated.
170 *> If NTYPES is smaller than the maximum number of types
171 *> defined (PARAMETER MAXTYP), then types NTYPES+1 through
172 *> MAXTYP will not be generated. If NTYPES is larger
173 *> than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
174 *> will be ignored.
175 *> \endverbatim
176 *>
177 *> \param[in,out] ISEED
178 *> \verbatim
179 *> ISEED is INTEGER array, dimension (4)
180 *> On entry ISEED specifies the seed of the random number
181 *> generator. The array elements should be between 0 and 4095;
182 *> if not they will be reduced mod 4096. Also, ISEED(4) must
183 *> be odd. The random number generator uses a linear
184 *> congruential sequence limited to small integers, and so
185 *> should produce machine independent random numbers. The
186 *> values of ISEED are changed on exit, and can be used in the
187 *> next call to ZCHKHBSTG to continue the same random number
188 *> sequence.
189 *> \endverbatim
190 *>
191 *> \param[in] THRESH
192 *> \verbatim
193 *> THRESH is DOUBLE PRECISION
194 *> A test will count as "failed" if the "error", computed as
195 *> described above, exceeds THRESH. Note that the error
196 *> is scaled to be O(1), so THRESH should be a reasonably
197 *> small multiple of 1, e.g., 10 or 100. In particular,
198 *> it should not depend on the precision (single vs. double)
199 *> or the size of the matrix. It must be at least zero.
200 *> \endverbatim
201 *>
202 *> \param[in] NOUNIT
203 *> \verbatim
204 *> NOUNIT is INTEGER
205 *> The FORTRAN unit number for printing out error messages
206 *> (e.g., if a routine returns IINFO not equal to 0.)
207 *> \endverbatim
208 *>
209 *> \param[in,out] A
210 *> \verbatim
211 *> A is COMPLEX*16 array, dimension
212 *> (LDA, max(NN))
213 *> Used to hold the matrix whose eigenvalues are to be
214 *> computed.
215 *> \endverbatim
216 *>
217 *> \param[in] LDA
218 *> \verbatim
219 *> LDA is INTEGER
220 *> The leading dimension of A. It must be at least 2 (not 1!)
221 *> and at least max( KK )+1.
222 *> \endverbatim
223 *>
224 *> \param[out] SD
225 *> \verbatim
226 *> SD is DOUBLE PRECISION array, dimension (max(NN))
227 *> Used to hold the diagonal of the tridiagonal matrix computed
228 *> by ZHBTRD.
229 *> \endverbatim
230 *>
231 *> \param[out] SE
232 *> \verbatim
233 *> SE is DOUBLE PRECISION array, dimension (max(NN))
234 *> Used to hold the off-diagonal of the tridiagonal matrix
235 *> computed by ZHBTRD.
236 *> \endverbatim
237 *>
238 *> \param[out] U
239 *> \verbatim
240 *> U is COMPLEX*16 array, dimension (LDU, max(NN))
241 *> Used to hold the unitary matrix computed by ZHBTRD.
242 *> \endverbatim
243 *>
244 *> \param[in] LDU
245 *> \verbatim
246 *> LDU is INTEGER
247 *> The leading dimension of U. It must be at least 1
248 *> and at least max( NN ).
249 *> \endverbatim
250 *>
251 *> \param[out] WORK
252 *> \verbatim
253 *> WORK is COMPLEX*16 array, dimension (LWORK)
254 *> \endverbatim
255 *>
256 *> \param[in] LWORK
257 *> \verbatim
258 *> LWORK is INTEGER
259 *> The number of entries in WORK. This must be at least
260 *> max( LDA+1, max(NN)+1 )*max(NN).
261 *> \endverbatim
262 *>
263 *> \param[out] RWORK
264 *> \verbatim
265 *> RWORK is DOUBLE PRECISION array
266 *> \endverbatim
267 *>
268 *> \param[out] RESULT
269 *> \verbatim
270 *> RESULT is DOUBLE PRECISION array, dimension (4)
271 *> The values computed by the tests described above.
272 *> The values are currently limited to 1/ulp, to avoid
273 *> overflow.
274 *> \endverbatim
275 *>
276 *> \param[out] INFO
277 *> \verbatim
278 *> INFO is INTEGER
279 *> If 0, then everything ran OK.
280 *>
281 *>-----------------------------------------------------------------------
282 *>
283 *> Some Local Variables and Parameters:
284 *> ---- ----- --------- --- ----------
285 *> ZERO, ONE Real 0 and 1.
286 *> MAXTYP The number of types defined.
287 *> NTEST The number of tests performed, or which can
288 *> be performed so far, for the current matrix.
289 *> NTESTT The total number of tests performed so far.
290 *> NMAX Largest value in NN.
291 *> NMATS The number of matrices generated so far.
292 *> NERRS The number of tests which have exceeded THRESH
293 *> so far.
294 *> COND, IMODE Values to be passed to the matrix generators.
295 *> ANORM Norm of A; passed to matrix generators.
296 *>
297 *> OVFL, UNFL Overflow and underflow thresholds.
298 *> ULP, ULPINV Finest relative precision and its inverse.
299 *> RTOVFL, RTUNFL Square roots of the previous 2 values.
300 *> The following four arrays decode JTYPE:
301 *> KTYPE(j) The general type (1-10) for type "j".
302 *> KMODE(j) The MODE value to be passed to the matrix
303 *> generator for type "j".
304 *> KMAGN(j) The order of magnitude ( O(1),
305 *> O(overflow^(1/2) ), O(underflow^(1/2) )
306 *> \endverbatim
307 *
308 * Authors:
309 * ========
310 *
311 *> \author Univ. of Tennessee
312 *> \author Univ. of California Berkeley
313 *> \author Univ. of Colorado Denver
314 *> \author NAG Ltd.
315 *
316 *> \date June 2017
317 *
318 *> \ingroup complex16_eig
319 *
320 * =====================================================================
321  SUBROUTINE zchkhb2stg( NSIZES, NN, NWDTHS, KK, NTYPES, DOTYPE,
322  $ ISEED, THRESH, NOUNIT, A, LDA, SD, SE, D1,
323  $ D2, D3, U, LDU, WORK, LWORK, RWORK, RESULT,
324  $ INFO )
325 *
326 * -- LAPACK test routine (version 3.7.1) --
327 * -- LAPACK is a software package provided by Univ. of Tennessee, --
328 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
329 * June 2017
330 *
331 * .. Scalar Arguments ..
332  INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES,
333  $ NWDTHS
334  DOUBLE PRECISION THRESH
335 * ..
336 * .. Array Arguments ..
337  LOGICAL DOTYPE( * )
338  INTEGER ISEED( 4 ), KK( * ), NN( * )
339  DOUBLE PRECISION RESULT( * ), RWORK( * ), SD( * ), SE( * ),
340  $ d1( * ), d2( * ), d3( * )
341  COMPLEX*16 A( LDA, * ), U( LDU, * ), WORK( * )
342 * ..
343 *
344 * =====================================================================
345 *
346 * .. Parameters ..
347  COMPLEX*16 CZERO, CONE
348  PARAMETER ( CZERO = ( 0.0d+0, 0.0d+0 ),
349  $ cone = ( 1.0d+0, 0.0d+0 ) )
350  DOUBLE PRECISION ZERO, ONE, TWO, TEN
351  parameter( zero = 0.0d+0, one = 1.0d+0, two = 2.0d+0,
352  $ ten = 10.0d+0 )
353  DOUBLE PRECISION HALF
354  parameter( half = one / two )
355  INTEGER MAXTYP
356  parameter( maxtyp = 15 )
357 * ..
358 * .. Local Scalars ..
359  LOGICAL BADNN, BADNNB
360  INTEGER I, IINFO, IMODE, ITYPE, J, JC, JCOL, JR, JSIZE,
361  $ JTYPE, JWIDTH, K, KMAX, LH, LW, MTYPES, N,
362  $ nerrs, nmats, nmax, ntest, ntestt
363  DOUBLE PRECISION ANINV, ANORM, COND, OVFL, RTOVFL, RTUNFL,
364  $ TEMP1, TEMP2, TEMP3, TEMP4, ULP, ULPINV, UNFL
365 * ..
366 * .. Local Arrays ..
367  INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KMAGN( MAXTYP ),
368  $ KMODE( MAXTYP ), KTYPE( MAXTYP )
369 * ..
370 * .. External Functions ..
371  DOUBLE PRECISION DLAMCH
372  EXTERNAL DLAMCH
373 * ..
374 * .. External Subroutines ..
375  EXTERNAL dlasum, xerbla, zhbt21, zhbtrd, zlacpy, zlaset,
376  $ zlatmr, zlatms, zhetrd_hb2st, zsteqr
377 * ..
378 * .. Intrinsic Functions ..
379  INTRINSIC abs, dble, dconjg, max, min, sqrt
380 * ..
381 * .. Data statements ..
382  DATA ktype / 1, 2, 5*4, 5*5, 3*8 /
383  DATA kmagn / 2*1, 1, 1, 1, 2, 3, 1, 1, 1, 2, 3, 1,
384  $ 2, 3 /
385  DATA kmode / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
386  $ 0, 0 /
387 * ..
388 * .. Executable Statements ..
389 *
390 * Check for errors
391 *
392  ntestt = 0
393  info = 0
394 *
395 * Important constants
396 *
397  badnn = .false.
398  nmax = 1
399  DO 10 j = 1, nsizes
400  nmax = max( nmax, nn( j ) )
401  IF( nn( j ).LT.0 )
402  $ badnn = .true.
403  10 CONTINUE
404 *
405  badnnb = .false.
406  kmax = 0
407  DO 20 j = 1, nsizes
408  kmax = max( kmax, kk( j ) )
409  IF( kk( j ).LT.0 )
410  $ badnnb = .true.
411  20 CONTINUE
412  kmax = min( nmax-1, kmax )
413 *
414 * Check for errors
415 *
416  IF( nsizes.LT.0 ) THEN
417  info = -1
418  ELSE IF( badnn ) THEN
419  info = -2
420  ELSE IF( nwdths.LT.0 ) THEN
421  info = -3
422  ELSE IF( badnnb ) THEN
423  info = -4
424  ELSE IF( ntypes.LT.0 ) THEN
425  info = -5
426  ELSE IF( lda.LT.kmax+1 ) THEN
427  info = -11
428  ELSE IF( ldu.LT.nmax ) THEN
429  info = -15
430  ELSE IF( ( max( lda, nmax )+1 )*nmax.GT.lwork ) THEN
431  info = -17
432  END IF
433 *
434  IF( info.NE.0 ) THEN
435  CALL xerbla( 'ZCHKHBSTG', -info )
436  RETURN
437  END IF
438 *
439 * Quick return if possible
440 *
441  IF( nsizes.EQ.0 .OR. ntypes.EQ.0 .OR. nwdths.EQ.0 )
442  $ RETURN
443 *
444 * More Important constants
445 *
446  unfl = dlamch( 'Safe minimum' )
447  ovfl = one / unfl
448  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
449  ulpinv = one / ulp
450  rtunfl = sqrt( unfl )
451  rtovfl = sqrt( ovfl )
452 *
453 * Loop over sizes, types
454 *
455  nerrs = 0
456  nmats = 0
457 *
458  DO 190 jsize = 1, nsizes
459  n = nn( jsize )
460  aninv = one / dble( max( 1, n ) )
461 *
462  DO 180 jwidth = 1, nwdths
463  k = kk( jwidth )
464  IF( k.GT.n )
465  $ GO TO 180
466  k = max( 0, min( n-1, k ) )
467 *
468  IF( nsizes.NE.1 ) THEN
469  mtypes = min( maxtyp, ntypes )
470  ELSE
471  mtypes = min( maxtyp+1, ntypes )
472  END IF
473 *
474  DO 170 jtype = 1, mtypes
475  IF( .NOT.dotype( jtype ) )
476  $ GO TO 170
477  nmats = nmats + 1
478  ntest = 0
479 *
480  DO 30 j = 1, 4
481  ioldsd( j ) = iseed( j )
482  30 CONTINUE
483 *
484 * Compute "A".
485 * Store as "Upper"; later, we will copy to other format.
486 *
487 * Control parameters:
488 *
489 * KMAGN KMODE KTYPE
490 * =1 O(1) clustered 1 zero
491 * =2 large clustered 2 identity
492 * =3 small exponential (none)
493 * =4 arithmetic diagonal, (w/ eigenvalues)
494 * =5 random log hermitian, w/ eigenvalues
495 * =6 random (none)
496 * =7 random diagonal
497 * =8 random hermitian
498 * =9 positive definite
499 * =10 diagonally dominant tridiagonal
500 *
501  IF( mtypes.GT.maxtyp )
502  $ GO TO 100
503 *
504  itype = ktype( jtype )
505  imode = kmode( jtype )
506 *
507 * Compute norm
508 *
509  GO TO ( 40, 50, 60 )kmagn( jtype )
510 *
511  40 CONTINUE
512  anorm = one
513  GO TO 70
514 *
515  50 CONTINUE
516  anorm = ( rtovfl*ulp )*aninv
517  GO TO 70
518 *
519  60 CONTINUE
520  anorm = rtunfl*n*ulpinv
521  GO TO 70
522 *
523  70 CONTINUE
524 *
525  CALL zlaset( 'Full', lda, n, czero, czero, a, lda )
526  iinfo = 0
527  IF( jtype.LE.15 ) THEN
528  cond = ulpinv
529  ELSE
530  cond = ulpinv*aninv / ten
531  END IF
532 *
533 * Special Matrices -- Identity & Jordan block
534 *
535 * Zero
536 *
537  IF( itype.EQ.1 ) THEN
538  iinfo = 0
539 *
540  ELSE IF( itype.EQ.2 ) THEN
541 *
542 * Identity
543 *
544  DO 80 jcol = 1, n
545  a( k+1, jcol ) = anorm
546  80 CONTINUE
547 *
548  ELSE IF( itype.EQ.4 ) THEN
549 *
550 * Diagonal Matrix, [Eigen]values Specified
551 *
552  CALL zlatms( n, n, 'S', iseed, 'H', rwork, imode,
553  $ cond, anorm, 0, 0, 'Q', a( k+1, 1 ), lda,
554  $ work, iinfo )
555 *
556  ELSE IF( itype.EQ.5 ) THEN
557 *
558 * Hermitian, eigenvalues specified
559 *
560  CALL zlatms( n, n, 'S', iseed, 'H', rwork, imode,
561  $ cond, anorm, k, k, 'Q', a, lda, work,
562  $ iinfo )
563 *
564  ELSE IF( itype.EQ.7 ) THEN
565 *
566 * Diagonal, random eigenvalues
567 *
568  CALL zlatmr( n, n, 'S', iseed, 'H', work, 6, one,
569  $ cone, 'T', 'N', work( n+1 ), 1, one,
570  $ work( 2*n+1 ), 1, one, 'N', idumma, 0, 0,
571  $ zero, anorm, 'Q', a( k+1, 1 ), lda,
572  $ idumma, iinfo )
573 *
574  ELSE IF( itype.EQ.8 ) THEN
575 *
576 * Hermitian, random eigenvalues
577 *
578  CALL zlatmr( n, n, 'S', iseed, 'H', work, 6, one,
579  $ cone, 'T', 'N', work( n+1 ), 1, one,
580  $ work( 2*n+1 ), 1, one, 'N', idumma, k, k,
581  $ zero, anorm, 'Q', a, lda, idumma, iinfo )
582 *
583  ELSE IF( itype.EQ.9 ) THEN
584 *
585 * Positive definite, eigenvalues specified.
586 *
587  CALL zlatms( n, n, 'S', iseed, 'P', rwork, imode,
588  $ cond, anorm, k, k, 'Q', a, lda,
589  $ work( n+1 ), iinfo )
590 *
591  ELSE IF( itype.EQ.10 ) THEN
592 *
593 * Positive definite tridiagonal, eigenvalues specified.
594 *
595  IF( n.GT.1 )
596  $ k = max( 1, k )
597  CALL zlatms( n, n, 'S', iseed, 'P', rwork, imode,
598  $ cond, anorm, 1, 1, 'Q', a( k, 1 ), lda,
599  $ work, iinfo )
600  DO 90 i = 2, n
601  temp1 = abs( a( k, i ) ) /
602  $ sqrt( abs( a( k+1, i-1 )*a( k+1, i ) ) )
603  IF( temp1.GT.half ) THEN
604  a( k, i ) = half*sqrt( abs( a( k+1,
605  $ i-1 )*a( k+1, i ) ) )
606  END IF
607  90 CONTINUE
608 *
609  ELSE
610 *
611  iinfo = 1
612  END IF
613 *
614  IF( iinfo.NE.0 ) THEN
615  WRITE( nounit, fmt = 9999 )'Generator', iinfo, n,
616  $ jtype, ioldsd
617  info = abs( iinfo )
618  RETURN
619  END IF
620 *
621  100 CONTINUE
622 *
623 * Call ZHBTRD to compute S and U from upper triangle.
624 *
625  CALL zlacpy( ' ', k+1, n, a, lda, work, lda )
626 *
627  ntest = 1
628  CALL zhbtrd( 'V', 'U', n, k, work, lda, sd, se, u, ldu,
629  $ work( lda*n+1 ), iinfo )
630 *
631  IF( iinfo.NE.0 ) THEN
632  WRITE( nounit, fmt = 9999 )'ZHBTRD(U)', iinfo, n,
633  $ jtype, ioldsd
634  info = abs( iinfo )
635  IF( iinfo.LT.0 ) THEN
636  RETURN
637  ELSE
638  result( 1 ) = ulpinv
639  GO TO 150
640  END IF
641  END IF
642 *
643 * Do tests 1 and 2
644 *
645  CALL zhbt21( 'Upper', n, k, 1, a, lda, sd, se, u, ldu,
646  $ work, rwork, result( 1 ) )
647 *
648 * Before converting A into lower for DSBTRD, run DSYTRD_SB2ST
649 * otherwise matrix A will be converted to lower and then need
650 * to be converted back to upper in order to run the upper case
651 * ofDSYTRD_SB2ST
652 *
653 * Compute D1 the eigenvalues resulting from the tridiagonal
654 * form using the DSBTRD and used as reference to compare
655 * with the DSYTRD_SB2ST routine
656 *
657 * Compute D1 from the DSBTRD and used as reference for the
658 * DSYTRD_SB2ST
659 *
660  CALL dcopy( n, sd, 1, d1, 1 )
661  IF( n.GT.0 )
662  $ CALL dcopy( n-1, se, 1, rwork, 1 )
663 *
664  CALL zsteqr( 'N', n, d1, rwork, work, ldu,
665  $ rwork( n+1 ), iinfo )
666  IF( iinfo.NE.0 ) THEN
667  WRITE( nounit, fmt = 9999 )'ZSTEQR(N)', iinfo, n,
668  $ jtype, ioldsd
669  info = abs( iinfo )
670  IF( iinfo.LT.0 ) THEN
671  RETURN
672  ELSE
673  result( 5 ) = ulpinv
674  GO TO 150
675  END IF
676  END IF
677 *
678 * DSYTRD_SB2ST Upper case is used to compute D2.
679 * Note to set SD and SE to zero to be sure not reusing
680 * the one from above. Compare it with D1 computed
681 * using the DSBTRD.
682 *
683  CALL dlaset( 'Full', n, 1, zero, zero, sd, 1 )
684  CALL dlaset( 'Full', n, 1, zero, zero, se, 1 )
685  CALL zlacpy( ' ', k+1, n, a, lda, u, ldu )
686  lh = max(1, 4*n)
687  lw = lwork - lh
688  CALL zhetrd_hb2st( 'N', 'N', "U", n, k, u, ldu, sd, se,
689  $ work, lh, work( lh+1 ), lw, iinfo )
690 *
691 * Compute D2 from the DSYTRD_SB2ST Upper case
692 *
693  CALL dcopy( n, sd, 1, d2, 1 )
694  IF( n.GT.0 )
695  $ CALL dcopy( n-1, se, 1, rwork, 1 )
696 *
697  CALL zsteqr( 'N', n, d2, rwork, work, ldu,
698  $ rwork( n+1 ), iinfo )
699  IF( iinfo.NE.0 ) THEN
700  WRITE( nounit, fmt = 9999 )'ZSTEQR(N)', iinfo, n,
701  $ jtype, ioldsd
702  info = abs( iinfo )
703  IF( iinfo.LT.0 ) THEN
704  RETURN
705  ELSE
706  result( 5 ) = ulpinv
707  GO TO 150
708  END IF
709  END IF
710 *
711 * Convert A from Upper-Triangle-Only storage to
712 * Lower-Triangle-Only storage.
713 *
714  DO 120 jc = 1, n
715  DO 110 jr = 0, min( k, n-jc )
716  a( jr+1, jc ) = dconjg( a( k+1-jr, jc+jr ) )
717  110 CONTINUE
718  120 CONTINUE
719  DO 140 jc = n + 1 - k, n
720  DO 130 jr = min( k, n-jc ) + 1, k
721  a( jr+1, jc ) = zero
722  130 CONTINUE
723  140 CONTINUE
724 *
725 * Call ZHBTRD to compute S and U from lower triangle
726 *
727  CALL zlacpy( ' ', k+1, n, a, lda, work, lda )
728 *
729  ntest = 3
730  CALL zhbtrd( 'V', 'L', n, k, work, lda, sd, se, u, ldu,
731  $ work( lda*n+1 ), iinfo )
732 *
733  IF( iinfo.NE.0 ) THEN
734  WRITE( nounit, fmt = 9999 )'ZHBTRD(L)', iinfo, n,
735  $ jtype, ioldsd
736  info = abs( iinfo )
737  IF( iinfo.LT.0 ) THEN
738  RETURN
739  ELSE
740  result( 3 ) = ulpinv
741  GO TO 150
742  END IF
743  END IF
744  ntest = 4
745 *
746 * Do tests 3 and 4
747 *
748  CALL zhbt21( 'Lower', n, k, 1, a, lda, sd, se, u, ldu,
749  $ work, rwork, result( 3 ) )
750 *
751 * DSYTRD_SB2ST Lower case is used to compute D3.
752 * Note to set SD and SE to zero to be sure not reusing
753 * the one from above. Compare it with D1 computed
754 * using the DSBTRD.
755 *
756  CALL dlaset( 'Full', n, 1, zero, zero, sd, 1 )
757  CALL dlaset( 'Full', n, 1, zero, zero, se, 1 )
758  CALL zlacpy( ' ', k+1, n, a, lda, u, ldu )
759  lh = max(1, 4*n)
760  lw = lwork - lh
761  CALL zhetrd_hb2st( 'N', 'N', "L", n, k, u, ldu, sd, se,
762  $ work, lh, work( lh+1 ), lw, iinfo )
763 *
764 * Compute D3 from the 2-stage Upper case
765 *
766  CALL dcopy( n, sd, 1, d3, 1 )
767  IF( n.GT.0 )
768  $ CALL dcopy( n-1, se, 1, rwork, 1 )
769 *
770  CALL zsteqr( 'N', n, d3, rwork, work, ldu,
771  $ rwork( n+1 ), iinfo )
772  IF( iinfo.NE.0 ) THEN
773  WRITE( nounit, fmt = 9999 )'ZSTEQR(N)', iinfo, n,
774  $ jtype, ioldsd
775  info = abs( iinfo )
776  IF( iinfo.LT.0 ) THEN
777  RETURN
778  ELSE
779  result( 6 ) = ulpinv
780  GO TO 150
781  END IF
782  END IF
783 *
784 *
785 * Do Tests 3 and 4 which are similar to 11 and 12 but with the
786 * D1 computed using the standard 1-stage reduction as reference
787 *
788  ntest = 6
789  temp1 = zero
790  temp2 = zero
791  temp3 = zero
792  temp4 = zero
793 *
794  DO 151 j = 1, n
795  temp1 = max( temp1, abs( d1( j ) ), abs( d2( j ) ) )
796  temp2 = max( temp2, abs( d1( j )-d2( j ) ) )
797  temp3 = max( temp3, abs( d1( j ) ), abs( d3( j ) ) )
798  temp4 = max( temp4, abs( d1( j )-d3( j ) ) )
799  151 CONTINUE
800 *
801  result(5) = temp2 / max( unfl, ulp*max( temp1, temp2 ) )
802  result(6) = temp4 / max( unfl, ulp*max( temp3, temp4 ) )
803 *
804 * End of Loop -- Check for RESULT(j) > THRESH
805 *
806  150 CONTINUE
807  ntestt = ntestt + ntest
808 *
809 * Print out tests which fail.
810 *
811  DO 160 jr = 1, ntest
812  IF( result( jr ).GE.thresh ) THEN
813 *
814 * If this is the first test to fail,
815 * print a header to the data file.
816 *
817  IF( nerrs.EQ.0 ) THEN
818  WRITE( nounit, fmt = 9998 )'ZHB'
819  WRITE( nounit, fmt = 9997 )
820  WRITE( nounit, fmt = 9996 )
821  WRITE( nounit, fmt = 9995 )'Hermitian'
822  WRITE( nounit, fmt = 9994 )'unitary', '*',
823  $ 'conjugate transpose', ( '*', j = 1, 6 )
824  END IF
825  nerrs = nerrs + 1
826  WRITE( nounit, fmt = 9993 )n, k, ioldsd, jtype,
827  $ jr, result( jr )
828  END IF
829  160 CONTINUE
830 *
831  170 CONTINUE
832  180 CONTINUE
833  190 CONTINUE
834 *
835 * Summary
836 *
837  CALL dlasum( 'ZHB', nounit, nerrs, ntestt )
838  RETURN
839 *
840  9999 FORMAT( ' ZCHKHBSTG: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
841  $ i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ), i5, ')' )
842  9998 FORMAT( / 1x, a3,
843  $ ' -- Complex Hermitian Banded Tridiagonal Reduction Routines'
844  $ )
845  9997 FORMAT( ' Matrix types (see DCHK23 for details): ' )
846 *
847  9996 FORMAT( / ' Special Matrices:',
848  $ / ' 1=Zero matrix. ',
849  $ ' 5=Diagonal: clustered entries.',
850  $ / ' 2=Identity matrix. ',
851  $ ' 6=Diagonal: large, evenly spaced.',
852  $ / ' 3=Diagonal: evenly spaced entries. ',
853  $ ' 7=Diagonal: small, evenly spaced.',
854  $ / ' 4=Diagonal: geometr. spaced entries.' )
855  9995 FORMAT( ' Dense ', a, ' Banded Matrices:',
856  $ / ' 8=Evenly spaced eigenvals. ',
857  $ ' 12=Small, evenly spaced eigenvals.',
858  $ / ' 9=Geometrically spaced eigenvals. ',
859  $ ' 13=Matrix with random O(1) entries.',
860  $ / ' 10=Clustered eigenvalues. ',
861  $ ' 14=Matrix with large random entries.',
862  $ / ' 11=Large, evenly spaced eigenvals. ',
863  $ ' 15=Matrix with small random entries.' )
864 *
865  9994 FORMAT( / ' Tests performed: (S is Tridiag, U is ', a, ',',
866  $ / 20x, a, ' means ', a, '.', / ' UPLO=''U'':',
867  $ / ' 1= | A - U S U', a1, ' | / ( |A| n ulp ) ',
868  $ ' 2= | I - U U', a1, ' | / ( n ulp )', / ' UPLO=''L'':',
869  $ / ' 3= | A - U S U', a1, ' | / ( |A| n ulp ) ',
870  $ ' 4= | I - U U', a1, ' | / ( n ulp )' / ' Eig check:',
871  $ /' 5= | D1 - D2', '', ' | / ( |D1| ulp ) ',
872  $ ' 6= | D1 - D3', '', ' | / ( |D1| ulp ) ' )
873  9993 FORMAT( ' N=', i5, ', K=', i4, ', seed=', 4( i4, ',' ), ' type ',
874  $ i2, ', test(', i2, ')=', g10.3 )
875 *
876 * End of ZCHKHBSTG
877 *
878  END
dlasum
subroutine dlasum(TYPE, IOUNIT, IE, NRUN)
DLASUM
Definition: dlasum.f:45
zhbt21
subroutine zhbt21(UPLO, N, KA, KS, A, LDA, D, E, U, LDU, WORK, RWORK, RESULT)
ZHBT21
Definition: zhbt21.f:154
zsteqr
subroutine zsteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
ZSTEQR
Definition: zsteqr.f:134
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:84
zhetrd_hb2st
subroutine zhetrd_hb2st(STAGE1, VECT, UPLO, N, KD, AB, LDAB, D, E, HOUS, LHOUS, WORK, LWORK, INFO)
ZHETRD_HB2ST reduces a complex Hermitian band matrix A to real symmetric tridiagonal form T
Definition: zhetrd_hb2st.F:232
zhbtrd
subroutine zhbtrd(VECT, UPLO, N, KD, AB, LDAB, D, E, Q, LDQ, WORK, INFO)
ZHBTRD
Definition: zhbtrd.f:165
zlaset
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: zlaset.f:108
zlacpy
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
zlatmr
subroutine zlatmr(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, RSIGN, GRADE, DL, MODEL, CONDL, DR, MODER, CONDR, PIVTNG, IPIVOT, KL, KU, SPARSE, ANORM, PACK, A, LDA, IWORK, INFO)
ZLATMR
Definition: zlatmr.f:492
zchkhb2stg
subroutine zchkhb2stg(NSIZES, NN, NWDTHS, KK, NTYPES, DOTYPE, ISEED, THRESH, NOUNIT, A, LDA, SD, SE, D1, D2, D3, U, LDU, WORK, LWORK, RWORK, RESULT, INFO)
ZCHKHBSTG
Definition: zchkhb2stg.f:325
dlaset
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
zlatms
subroutine zlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
ZLATMS
Definition: zlatms.f:334