LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zchkbd.f
Go to the documentation of this file.
1 *> \brief \b ZCHKBD
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 ZCHKBD( NSIZES, MVAL, NVAL, NTYPES, DOTYPE, NRHS,
12 * ISEED, THRESH, A, LDA, BD, BE, S1, S2, X, LDX,
13 * Y, Z, Q, LDQ, PT, LDPT, U, VT, WORK, LWORK,
14 * RWORK, NOUT, INFO )
15 *
16 * .. Scalar Arguments ..
17 * INTEGER INFO, LDA, LDPT, LDQ, LDX, LWORK, NOUT, NRHS,
18 * $ NSIZES, NTYPES
19 * DOUBLE PRECISION THRESH
20 * ..
21 * .. Array Arguments ..
22 * LOGICAL DOTYPE( * )
23 * INTEGER ISEED( 4 ), MVAL( * ), NVAL( * )
24 * DOUBLE PRECISION BD( * ), BE( * ), RWORK( * ), S1( * ), S2( * )
25 * COMPLEX*16 A( LDA, * ), PT( LDPT, * ), Q( LDQ, * ),
26 * $ U( LDPT, * ), VT( LDPT, * ), WORK( * ),
27 * $ X( LDX, * ), Y( LDX, * ), Z( LDX, * )
28 * ..
29 *
30 *
31 *> \par Purpose:
32 * =============
33 *>
34 *> \verbatim
35 *>
36 *> ZCHKBD checks the singular value decomposition (SVD) routines.
37 *>
38 *> ZGEBRD reduces a complex general m by n matrix A to real upper or
39 *> lower bidiagonal form by an orthogonal transformation: Q' * A * P = B
40 *> (or A = Q * B * P'). The matrix B is upper bidiagonal if m >= n
41 *> and lower bidiagonal if m < n.
42 *>
43 *> ZUNGBR generates the orthogonal matrices Q and P' from ZGEBRD.
44 *> Note that Q and P are not necessarily square.
45 *>
46 *> ZBDSQR computes the singular value decomposition of the bidiagonal
47 *> matrix B as B = U S V'. It is called three times to compute
48 *> 1) B = U S1 V', where S1 is the diagonal matrix of singular
49 *> values and the columns of the matrices U and V are the left
50 *> and right singular vectors, respectively, of B.
51 *> 2) Same as 1), but the singular values are stored in S2 and the
52 *> singular vectors are not computed.
53 *> 3) A = (UQ) S (P'V'), the SVD of the original matrix A.
54 *> In addition, ZBDSQR has an option to apply the left orthogonal matrix
55 *> U to a matrix X, useful in least squares applications.
56 *>
57 *> For each pair of matrix dimensions (M,N) and each selected matrix
58 *> type, an M by N matrix A and an M by NRHS matrix X are generated.
59 *> The problem dimensions are as follows
60 *> A: M x N
61 *> Q: M x min(M,N) (but M x M if NRHS > 0)
62 *> P: min(M,N) x N
63 *> B: min(M,N) x min(M,N)
64 *> U, V: min(M,N) x min(M,N)
65 *> S1, S2 diagonal, order min(M,N)
66 *> X: M x NRHS
67 *>
68 *> For each generated matrix, 14 tests are performed:
69 *>
70 *> Test ZGEBRD and ZUNGBR
71 *>
72 *> (1) | A - Q B PT | / ( |A| max(M,N) ulp ), PT = P'
73 *>
74 *> (2) | I - Q' Q | / ( M ulp )
75 *>
76 *> (3) | I - PT PT' | / ( N ulp )
77 *>
78 *> Test ZBDSQR on bidiagonal matrix B
79 *>
80 *> (4) | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V'
81 *>
82 *> (5) | Y - U Z | / ( |Y| max(min(M,N),k) ulp ), where Y = Q' X
83 *> and Z = U' Y.
84 *> (6) | I - U' U | / ( min(M,N) ulp )
85 *>
86 *> (7) | I - VT VT' | / ( min(M,N) ulp )
87 *>
88 *> (8) S1 contains min(M,N) nonnegative values in decreasing order.
89 *> (Return 0 if true, 1/ULP if false.)
90 *>
91 *> (9) 0 if the true singular values of B are within THRESH of
92 *> those in S1. 2*THRESH if they are not. (Tested using
93 *> DSVDCH)
94 *>
95 *> (10) | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
96 *> computing U and V.
97 *>
98 *> Test ZBDSQR on matrix A
99 *>
100 *> (11) | A - (QU) S (VT PT) | / ( |A| max(M,N) ulp )
101 *>
102 *> (12) | X - (QU) Z | / ( |X| max(M,k) ulp )
103 *>
104 *> (13) | I - (QU)'(QU) | / ( M ulp )
105 *>
106 *> (14) | I - (VT PT) (PT'VT') | / ( N ulp )
107 *>
108 *> The possible matrix types are
109 *>
110 *> (1) The zero matrix.
111 *> (2) The identity matrix.
112 *>
113 *> (3) A diagonal matrix with evenly spaced entries
114 *> 1, ..., ULP and random signs.
115 *> (ULP = (first number larger than 1) - 1 )
116 *> (4) A diagonal matrix with geometrically spaced entries
117 *> 1, ..., ULP and random signs.
118 *> (5) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
119 *> and random signs.
120 *>
121 *> (6) Same as (3), but multiplied by SQRT( overflow threshold )
122 *> (7) Same as (3), but multiplied by SQRT( underflow threshold )
123 *>
124 *> (8) A matrix of the form U D V, where U and V are orthogonal and
125 *> D has evenly spaced entries 1, ..., ULP with random signs
126 *> on the diagonal.
127 *>
128 *> (9) A matrix of the form U D V, where U and V are orthogonal and
129 *> D has geometrically spaced entries 1, ..., ULP with random
130 *> signs on the diagonal.
131 *>
132 *> (10) A matrix of the form U D V, where U and V are orthogonal and
133 *> D has "clustered" entries 1, ULP,..., ULP with random
134 *> signs on the diagonal.
135 *>
136 *> (11) Same as (8), but multiplied by SQRT( overflow threshold )
137 *> (12) Same as (8), but multiplied by SQRT( underflow threshold )
138 *>
139 *> (13) Rectangular matrix with random entries chosen from (-1,1).
140 *> (14) Same as (13), but multiplied by SQRT( overflow threshold )
141 *> (15) Same as (13), but multiplied by SQRT( underflow threshold )
142 *>
143 *> Special case:
144 *> (16) A bidiagonal matrix with random entries chosen from a
145 *> logarithmic distribution on [ulp^2,ulp^(-2)] (I.e., each
146 *> entry is e^x, where x is chosen uniformly on
147 *> [ 2 log(ulp), -2 log(ulp) ] .) For *this* type:
148 *> (a) ZGEBRD is not called to reduce it to bidiagonal form.
149 *> (b) the bidiagonal is min(M,N) x min(M,N); if M<N, the
150 *> matrix will be lower bidiagonal, otherwise upper.
151 *> (c) only tests 5--8 and 14 are performed.
152 *>
153 *> A subset of the full set of matrix types may be selected through
154 *> the logical array DOTYPE.
155 *> \endverbatim
156 *
157 * Arguments:
158 * ==========
159 *
160 *> \param[in] NSIZES
161 *> \verbatim
162 *> NSIZES is INTEGER
163 *> The number of values of M and N contained in the vectors
164 *> MVAL and NVAL. The matrix sizes are used in pairs (M,N).
165 *> \endverbatim
166 *>
167 *> \param[in] MVAL
168 *> \verbatim
169 *> MVAL is INTEGER array, dimension (NM)
170 *> The values of the matrix row dimension M.
171 *> \endverbatim
172 *>
173 *> \param[in] NVAL
174 *> \verbatim
175 *> NVAL is INTEGER array, dimension (NM)
176 *> The values of the matrix column dimension N.
177 *> \endverbatim
178 *>
179 *> \param[in] NTYPES
180 *> \verbatim
181 *> NTYPES is INTEGER
182 *> The number of elements in DOTYPE. If it is zero, ZCHKBD
183 *> does nothing. It must be at least zero. If it is MAXTYP+1
184 *> and NSIZES is 1, then an additional type, MAXTYP+1 is
185 *> defined, which is to use whatever matrices are in A and B.
186 *> This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
187 *> DOTYPE(MAXTYP+1) is .TRUE. .
188 *> \endverbatim
189 *>
190 *> \param[in] DOTYPE
191 *> \verbatim
192 *> DOTYPE is LOGICAL array, dimension (NTYPES)
193 *> If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix
194 *> of type j will be generated. If NTYPES is smaller than the
195 *> maximum number of types defined (PARAMETER MAXTYP), then
196 *> types NTYPES+1 through MAXTYP will not be generated. If
197 *> NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through
198 *> DOTYPE(NTYPES) will be ignored.
199 *> \endverbatim
200 *>
201 *> \param[in] NRHS
202 *> \verbatim
203 *> NRHS is INTEGER
204 *> The number of columns in the "right-hand side" matrices X, Y,
205 *> and Z, used in testing ZBDSQR. If NRHS = 0, then the
206 *> operations on the right-hand side will not be tested.
207 *> NRHS must be at least 0.
208 *> \endverbatim
209 *>
210 *> \param[in,out] ISEED
211 *> \verbatim
212 *> ISEED is INTEGER array, dimension (4)
213 *> On entry ISEED specifies the seed of the random number
214 *> generator. The array elements should be between 0 and 4095;
215 *> if not they will be reduced mod 4096. Also, ISEED(4) must
216 *> be odd. The values of ISEED are changed on exit, and can be
217 *> used in the next call to ZCHKBD to continue the same random
218 *> number sequence.
219 *> \endverbatim
220 *>
221 *> \param[in] THRESH
222 *> \verbatim
223 *> THRESH is DOUBLE PRECISION
224 *> The threshold value for the test ratios. A result is
225 *> included in the output file if RESULT >= THRESH. To have
226 *> every test ratio printed, use THRESH = 0. Note that the
227 *> expected value of the test ratios is O(1), so THRESH should
228 *> be a reasonably small multiple of 1, e.g., 10 or 100.
229 *> \endverbatim
230 *>
231 *> \param[out] A
232 *> \verbatim
233 *> A is COMPLEX*16 array, dimension (LDA,NMAX)
234 *> where NMAX is the maximum value of N in NVAL.
235 *> \endverbatim
236 *>
237 *> \param[in] LDA
238 *> \verbatim
239 *> LDA is INTEGER
240 *> The leading dimension of the array A. LDA >= max(1,MMAX),
241 *> where MMAX is the maximum value of M in MVAL.
242 *> \endverbatim
243 *>
244 *> \param[out] BD
245 *> \verbatim
246 *> BD is DOUBLE PRECISION array, dimension
247 *> (max(min(MVAL(j),NVAL(j))))
248 *> \endverbatim
249 *>
250 *> \param[out] BE
251 *> \verbatim
252 *> BE is DOUBLE PRECISION array, dimension
253 *> (max(min(MVAL(j),NVAL(j))))
254 *> \endverbatim
255 *>
256 *> \param[out] S1
257 *> \verbatim
258 *> S1 is DOUBLE PRECISION array, dimension
259 *> (max(min(MVAL(j),NVAL(j))))
260 *> \endverbatim
261 *>
262 *> \param[out] S2
263 *> \verbatim
264 *> S2 is DOUBLE PRECISION array, dimension
265 *> (max(min(MVAL(j),NVAL(j))))
266 *> \endverbatim
267 *>
268 *> \param[out] X
269 *> \verbatim
270 *> X is COMPLEX*16 array, dimension (LDX,NRHS)
271 *> \endverbatim
272 *>
273 *> \param[in] LDX
274 *> \verbatim
275 *> LDX is INTEGER
276 *> The leading dimension of the arrays X, Y, and Z.
277 *> LDX >= max(1,MMAX).
278 *> \endverbatim
279 *>
280 *> \param[out] Y
281 *> \verbatim
282 *> Y is COMPLEX*16 array, dimension (LDX,NRHS)
283 *> \endverbatim
284 *>
285 *> \param[out] Z
286 *> \verbatim
287 *> Z is COMPLEX*16 array, dimension (LDX,NRHS)
288 *> \endverbatim
289 *>
290 *> \param[out] Q
291 *> \verbatim
292 *> Q is COMPLEX*16 array, dimension (LDQ,MMAX)
293 *> \endverbatim
294 *>
295 *> \param[in] LDQ
296 *> \verbatim
297 *> LDQ is INTEGER
298 *> The leading dimension of the array Q. LDQ >= max(1,MMAX).
299 *> \endverbatim
300 *>
301 *> \param[out] PT
302 *> \verbatim
303 *> PT is COMPLEX*16 array, dimension (LDPT,NMAX)
304 *> \endverbatim
305 *>
306 *> \param[in] LDPT
307 *> \verbatim
308 *> LDPT is INTEGER
309 *> The leading dimension of the arrays PT, U, and V.
310 *> LDPT >= max(1, max(min(MVAL(j),NVAL(j)))).
311 *> \endverbatim
312 *>
313 *> \param[out] U
314 *> \verbatim
315 *> U is COMPLEX*16 array, dimension
316 *> (LDPT,max(min(MVAL(j),NVAL(j))))
317 *> \endverbatim
318 *>
319 *> \param[out] VT
320 *> \verbatim
321 *> VT is COMPLEX*16 array, dimension
322 *> (LDPT,max(min(MVAL(j),NVAL(j))))
323 *> \endverbatim
324 *>
325 *> \param[out] WORK
326 *> \verbatim
327 *> WORK is COMPLEX*16 array, dimension (LWORK)
328 *> \endverbatim
329 *>
330 *> \param[in] LWORK
331 *> \verbatim
332 *> LWORK is INTEGER
333 *> The number of entries in WORK. This must be at least
334 *> 3(M+N) and M(M + max(M,N,k) + 1) + N*min(M,N) for all
335 *> pairs (M,N)=(MM(j),NN(j))
336 *> \endverbatim
337 *>
338 *> \param[out] RWORK
339 *> \verbatim
340 *> RWORK is DOUBLE PRECISION array, dimension
341 *> (5*max(min(M,N)))
342 *> \endverbatim
343 *>
344 *> \param[in] NOUT
345 *> \verbatim
346 *> NOUT is INTEGER
347 *> The FORTRAN unit number for printing out error messages
348 *> (e.g., if a routine returns IINFO not equal to 0.)
349 *> \endverbatim
350 *>
351 *> \param[out] INFO
352 *> \verbatim
353 *> INFO is INTEGER
354 *> If 0, then everything ran OK.
355 *> -1: NSIZES < 0
356 *> -2: Some MM(j) < 0
357 *> -3: Some NN(j) < 0
358 *> -4: NTYPES < 0
359 *> -6: NRHS < 0
360 *> -8: THRESH < 0
361 *> -11: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ).
362 *> -17: LDB < 1 or LDB < MMAX.
363 *> -21: LDQ < 1 or LDQ < MMAX.
364 *> -23: LDP < 1 or LDP < MNMAX.
365 *> -27: LWORK too small.
366 *> If ZLATMR, CLATMS, ZGEBRD, ZUNGBR, or ZBDSQR,
367 *> returns an error code, the
368 *> absolute value of it is returned.
369 *>
370 *>-----------------------------------------------------------------------
371 *>
372 *> Some Local Variables and Parameters:
373 *> ---- ----- --------- --- ----------
374 *>
375 *> ZERO, ONE Real 0 and 1.
376 *> MAXTYP The number of types defined.
377 *> NTEST The number of tests performed, or which can
378 *> be performed so far, for the current matrix.
379 *> MMAX Largest value in NN.
380 *> NMAX Largest value in NN.
381 *> MNMIN min(MM(j), NN(j)) (the dimension of the bidiagonal
382 *> matrix.)
383 *> MNMAX The maximum value of MNMIN for j=1,...,NSIZES.
384 *> NFAIL The number of tests which have exceeded THRESH
385 *> COND, IMODE Values to be passed to the matrix generators.
386 *> ANORM Norm of A; passed to matrix generators.
387 *>
388 *> OVFL, UNFL Overflow and underflow thresholds.
389 *> RTOVFL, RTUNFL Square roots of the previous 2 values.
390 *> ULP, ULPINV Finest relative precision and its inverse.
391 *>
392 *> The following four arrays decode JTYPE:
393 *> KTYPE(j) The general type (1-10) for type "j".
394 *> KMODE(j) The MODE value to be passed to the matrix
395 *> generator for type "j".
396 *> KMAGN(j) The order of magnitude ( O(1),
397 *> O(overflow^(1/2) ), O(underflow^(1/2) )
398 *> \endverbatim
399 *
400 * Authors:
401 * ========
402 *
403 *> \author Univ. of Tennessee
404 *> \author Univ. of California Berkeley
405 *> \author Univ. of Colorado Denver
406 *> \author NAG Ltd.
407 *
408 *> \date November 2011
409 *
410 *> \ingroup complex16_eig
411 *
412 * =====================================================================
413  SUBROUTINE zchkbd( NSIZES, MVAL, NVAL, NTYPES, DOTYPE, NRHS,
414  $ iseed, thresh, a, lda, bd, be, s1, s2, x, ldx,
415  $ y, z, q, ldq, pt, ldpt, u, vt, work, lwork,
416  $ rwork, nout, info )
417 *
418 * -- LAPACK test routine (version 3.4.0) --
419 * -- LAPACK is a software package provided by Univ. of Tennessee, --
420 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
421 * November 2011
422 *
423 * .. Scalar Arguments ..
424  INTEGER info, lda, ldpt, ldq, ldx, lwork, nout, nrhs,
425  $ nsizes, ntypes
426  DOUBLE PRECISION thresh
427 * ..
428 * .. Array Arguments ..
429  LOGICAL dotype( * )
430  INTEGER iseed( 4 ), mval( * ), nval( * )
431  DOUBLE PRECISION bd( * ), be( * ), rwork( * ), s1( * ), s2( * )
432  COMPLEX*16 a( lda, * ), pt( ldpt, * ), q( ldq, * ),
433  $ u( ldpt, * ), vt( ldpt, * ), work( * ),
434  $ x( ldx, * ), y( ldx, * ), z( ldx, * )
435 * ..
436 *
437 * ======================================================================
438 *
439 * .. Parameters ..
440  DOUBLE PRECISION zero, one, two, half
441  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
442  $ half = 0.5d0 )
443  COMPLEX*16 czero, cone
444  parameter( czero = ( 0.0d+0, 0.0d+0 ),
445  $ cone = ( 1.0d+0, 0.0d+0 ) )
446  INTEGER maxtyp
447  parameter( maxtyp = 16 )
448 * ..
449 * .. Local Scalars ..
450  LOGICAL badmm, badnn, bidiag
451  CHARACTER uplo
452  CHARACTER*3 path
453  INTEGER i, iinfo, imode, itype, j, jcol, jsize, jtype,
454  $ log2ui, m, minwrk, mmax, mnmax, mnmin, mq,
455  $ mtypes, n, nfail, nmax, ntest
456  DOUBLE PRECISION amninv, anorm, cond, ovfl, rtovfl, rtunfl,
457  $ temp1, temp2, ulp, ulpinv, unfl
458 * ..
459 * .. Local Arrays ..
460  INTEGER ioldsd( 4 ), iwork( 1 ), kmagn( maxtyp ),
461  $ kmode( maxtyp ), ktype( maxtyp )
462  DOUBLE PRECISION dumma( 1 ), result( 14 )
463 * ..
464 * .. External Functions ..
465  DOUBLE PRECISION dlamch, dlarnd
466  EXTERNAL dlamch, dlarnd
467 * ..
468 * .. External Subroutines ..
469  EXTERNAL alasum, dcopy, dlabad, dlahd2, dsvdch, xerbla,
472 * ..
473 * .. Intrinsic Functions ..
474  INTRINSIC abs, exp, int, log, max, min, sqrt
475 * ..
476 * .. Scalars in Common ..
477  LOGICAL lerr, ok
478  CHARACTER*32 srnamt
479  INTEGER infot, nunit
480 * ..
481 * .. Common blocks ..
482  common / infoc / infot, nunit, ok, lerr
483  common / srnamc / srnamt
484 * ..
485 * .. Data statements ..
486  DATA ktype / 1, 2, 5*4, 5*6, 3*9, 10 /
487  DATA kmagn / 2*1, 3*1, 2, 3, 3*1, 2, 3, 1, 2, 3, 0 /
488  DATA kmode / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
489  $ 0, 0, 0 /
490 * ..
491 * .. Executable Statements ..
492 *
493 * Check for errors
494 *
495  info = 0
496 *
497  badmm = .false.
498  badnn = .false.
499  mmax = 1
500  nmax = 1
501  mnmax = 1
502  minwrk = 1
503  DO 10 j = 1, nsizes
504  mmax = max( mmax, mval( j ) )
505  IF( mval( j ).LT.0 )
506  $ badmm = .true.
507  nmax = max( nmax, nval( j ) )
508  IF( nval( j ).LT.0 )
509  $ badnn = .true.
510  mnmax = max( mnmax, min( mval( j ), nval( j ) ) )
511  minwrk = max( minwrk, 3*( mval( j )+nval( j ) ),
512  $ mval( j )*( mval( j )+max( mval( j ), nval( j ),
513  $ nrhs )+1 )+nval( j )*min( nval( j ), mval( j ) ) )
514  10 continue
515 *
516 * Check for errors
517 *
518  IF( nsizes.LT.0 ) THEN
519  info = -1
520  ELSE IF( badmm ) THEN
521  info = -2
522  ELSE IF( badnn ) THEN
523  info = -3
524  ELSE IF( ntypes.LT.0 ) THEN
525  info = -4
526  ELSE IF( nrhs.LT.0 ) THEN
527  info = -6
528  ELSE IF( lda.LT.mmax ) THEN
529  info = -11
530  ELSE IF( ldx.LT.mmax ) THEN
531  info = -17
532  ELSE IF( ldq.LT.mmax ) THEN
533  info = -21
534  ELSE IF( ldpt.LT.mnmax ) THEN
535  info = -23
536  ELSE IF( minwrk.GT.lwork ) THEN
537  info = -27
538  END IF
539 *
540  IF( info.NE.0 ) THEN
541  CALL xerbla( 'ZCHKBD', -info )
542  return
543  END IF
544 *
545 * Initialize constants
546 *
547  path( 1: 1 ) = 'Zomplex precision'
548  path( 2: 3 ) = 'BD'
549  nfail = 0
550  ntest = 0
551  unfl = dlamch( 'Safe minimum' )
552  ovfl = dlamch( 'Overflow' )
553  CALL dlabad( unfl, ovfl )
554  ulp = dlamch( 'Precision' )
555  ulpinv = one / ulp
556  log2ui = int( log( ulpinv ) / log( two ) )
557  rtunfl = sqrt( unfl )
558  rtovfl = sqrt( ovfl )
559  infot = 0
560 *
561 * Loop over sizes, types
562 *
563  DO 180 jsize = 1, nsizes
564  m = mval( jsize )
565  n = nval( jsize )
566  mnmin = min( m, n )
567  amninv = one / max( m, n, 1 )
568 *
569  IF( nsizes.NE.1 ) THEN
570  mtypes = min( maxtyp, ntypes )
571  ELSE
572  mtypes = min( maxtyp+1, ntypes )
573  END IF
574 *
575  DO 170 jtype = 1, mtypes
576  IF( .NOT.dotype( jtype ) )
577  $ go to 170
578 *
579  DO 20 j = 1, 4
580  ioldsd( j ) = iseed( j )
581  20 continue
582 *
583  DO 30 j = 1, 14
584  result( j ) = -one
585  30 continue
586 *
587  uplo = ' '
588 *
589 * Compute "A"
590 *
591 * Control parameters:
592 *
593 * KMAGN KMODE KTYPE
594 * =1 O(1) clustered 1 zero
595 * =2 large clustered 2 identity
596 * =3 small exponential (none)
597 * =4 arithmetic diagonal, (w/ eigenvalues)
598 * =5 random symmetric, w/ eigenvalues
599 * =6 nonsymmetric, w/ singular values
600 * =7 random diagonal
601 * =8 random symmetric
602 * =9 random nonsymmetric
603 * =10 random bidiagonal (log. distrib.)
604 *
605  IF( mtypes.GT.maxtyp )
606  $ go to 100
607 *
608  itype = ktype( jtype )
609  imode = kmode( jtype )
610 *
611 * Compute norm
612 *
613  go to( 40, 50, 60 )kmagn( jtype )
614 *
615  40 continue
616  anorm = one
617  go to 70
618 *
619  50 continue
620  anorm = ( rtovfl*ulp )*amninv
621  go to 70
622 *
623  60 continue
624  anorm = rtunfl*max( m, n )*ulpinv
625  go to 70
626 *
627  70 continue
628 *
629  CALL zlaset( 'Full', lda, n, czero, czero, a, lda )
630  iinfo = 0
631  cond = ulpinv
632 *
633  bidiag = .false.
634  IF( itype.EQ.1 ) THEN
635 *
636 * Zero matrix
637 *
638  iinfo = 0
639 *
640  ELSE IF( itype.EQ.2 ) THEN
641 *
642 * Identity
643 *
644  DO 80 jcol = 1, mnmin
645  a( jcol, jcol ) = anorm
646  80 continue
647 *
648  ELSE IF( itype.EQ.4 ) THEN
649 *
650 * Diagonal Matrix, [Eigen]values Specified
651 *
652  CALL zlatms( mnmin, mnmin, 'S', iseed, 'N', rwork, imode,
653  $ cond, anorm, 0, 0, 'N', a, lda, work,
654  $ iinfo )
655 *
656  ELSE IF( itype.EQ.5 ) THEN
657 *
658 * Symmetric, eigenvalues specified
659 *
660  CALL zlatms( mnmin, mnmin, 'S', iseed, 'S', rwork, imode,
661  $ cond, anorm, m, n, 'N', a, lda, work,
662  $ iinfo )
663 *
664  ELSE IF( itype.EQ.6 ) THEN
665 *
666 * Nonsymmetric, singular values specified
667 *
668  CALL zlatms( m, n, 'S', iseed, 'N', rwork, imode, cond,
669  $ anorm, m, n, 'N', a, lda, work, iinfo )
670 *
671  ELSE IF( itype.EQ.7 ) THEN
672 *
673 * Diagonal, random entries
674 *
675  CALL zlatmr( mnmin, mnmin, 'S', iseed, 'N', work, 6, one,
676  $ cone, 'T', 'N', work( mnmin+1 ), 1, one,
677  $ work( 2*mnmin+1 ), 1, one, 'N', iwork, 0, 0,
678  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
679 *
680  ELSE IF( itype.EQ.8 ) THEN
681 *
682 * Symmetric, random entries
683 *
684  CALL zlatmr( mnmin, mnmin, 'S', iseed, 'S', work, 6, one,
685  $ cone, 'T', 'N', work( mnmin+1 ), 1, one,
686  $ work( m+mnmin+1 ), 1, one, 'N', iwork, m, n,
687  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
688 *
689  ELSE IF( itype.EQ.9 ) THEN
690 *
691 * Nonsymmetric, random entries
692 *
693  CALL zlatmr( m, n, 'S', iseed, 'N', work, 6, one, cone,
694  $ 'T', 'N', work( mnmin+1 ), 1, one,
695  $ work( m+mnmin+1 ), 1, one, 'N', iwork, m, n,
696  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
697 *
698  ELSE IF( itype.EQ.10 ) THEN
699 *
700 * Bidiagonal, random entries
701 *
702  temp1 = -two*log( ulp )
703  DO 90 j = 1, mnmin
704  bd( j ) = exp( temp1*dlarnd( 2, iseed ) )
705  IF( j.LT.mnmin )
706  $ be( j ) = exp( temp1*dlarnd( 2, iseed ) )
707  90 continue
708 *
709  iinfo = 0
710  bidiag = .true.
711  IF( m.GE.n ) THEN
712  uplo = 'U'
713  ELSE
714  uplo = 'L'
715  END IF
716  ELSE
717  iinfo = 1
718  END IF
719 *
720  IF( iinfo.EQ.0 ) THEN
721 *
722 * Generate Right-Hand Side
723 *
724  IF( bidiag ) THEN
725  CALL zlatmr( mnmin, nrhs, 'S', iseed, 'N', work, 6,
726  $ one, cone, 'T', 'N', work( mnmin+1 ), 1,
727  $ one, work( 2*mnmin+1 ), 1, one, 'N',
728  $ iwork, mnmin, nrhs, zero, one, 'NO', y,
729  $ ldx, iwork, iinfo )
730  ELSE
731  CALL zlatmr( m, nrhs, 'S', iseed, 'N', work, 6, one,
732  $ cone, 'T', 'N', work( m+1 ), 1, one,
733  $ work( 2*m+1 ), 1, one, 'N', iwork, m,
734  $ nrhs, zero, one, 'NO', x, ldx, iwork,
735  $ iinfo )
736  END IF
737  END IF
738 *
739 * Error Exit
740 *
741  IF( iinfo.NE.0 ) THEN
742  WRITE( nout, fmt = 9998 )'Generator', iinfo, m, n,
743  $ jtype, ioldsd
744  info = abs( iinfo )
745  return
746  END IF
747 *
748  100 continue
749 *
750 * Call ZGEBRD and ZUNGBR to compute B, Q, and P, do tests.
751 *
752  IF( .NOT.bidiag ) THEN
753 *
754 * Compute transformations to reduce A to bidiagonal form:
755 * B := Q' * A * P.
756 *
757  CALL zlacpy( ' ', m, n, a, lda, q, ldq )
758  CALL zgebrd( m, n, q, ldq, bd, be, work, work( mnmin+1 ),
759  $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
760 *
761 * Check error code from ZGEBRD.
762 *
763  IF( iinfo.NE.0 ) THEN
764  WRITE( nout, fmt = 9998 )'ZGEBRD', iinfo, m, n,
765  $ jtype, ioldsd
766  info = abs( iinfo )
767  return
768  END IF
769 *
770  CALL zlacpy( ' ', m, n, q, ldq, pt, ldpt )
771  IF( m.GE.n ) THEN
772  uplo = 'U'
773  ELSE
774  uplo = 'L'
775  END IF
776 *
777 * Generate Q
778 *
779  mq = m
780  IF( nrhs.LE.0 )
781  $ mq = mnmin
782  CALL zungbr( 'Q', m, mq, n, q, ldq, work,
783  $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
784 *
785 * Check error code from ZUNGBR.
786 *
787  IF( iinfo.NE.0 ) THEN
788  WRITE( nout, fmt = 9998 )'ZUNGBR(Q)', iinfo, m, n,
789  $ jtype, ioldsd
790  info = abs( iinfo )
791  return
792  END IF
793 *
794 * Generate P'
795 *
796  CALL zungbr( 'P', mnmin, n, m, pt, ldpt, work( mnmin+1 ),
797  $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
798 *
799 * Check error code from ZUNGBR.
800 *
801  IF( iinfo.NE.0 ) THEN
802  WRITE( nout, fmt = 9998 )'ZUNGBR(P)', iinfo, m, n,
803  $ jtype, ioldsd
804  info = abs( iinfo )
805  return
806  END IF
807 *
808 * Apply Q' to an M by NRHS matrix X: Y := Q' * X.
809 *
810  CALL zgemm( 'Conjugate transpose', 'No transpose', m,
811  $ nrhs, m, cone, q, ldq, x, ldx, czero, y,
812  $ ldx )
813 *
814 * Test 1: Check the decomposition A := Q * B * PT
815 * 2: Check the orthogonality of Q
816 * 3: Check the orthogonality of PT
817 *
818  CALL zbdt01( m, n, 1, a, lda, q, ldq, bd, be, pt, ldpt,
819  $ work, rwork, result( 1 ) )
820  CALL zunt01( 'Columns', m, mq, q, ldq, work, lwork,
821  $ rwork, result( 2 ) )
822  CALL zunt01( 'Rows', mnmin, n, pt, ldpt, work, lwork,
823  $ rwork, result( 3 ) )
824  END IF
825 *
826 * Use ZBDSQR to form the SVD of the bidiagonal matrix B:
827 * B := U * S1 * VT, and compute Z = U' * Y.
828 *
829  CALL dcopy( mnmin, bd, 1, s1, 1 )
830  IF( mnmin.GT.0 )
831  $ CALL dcopy( mnmin-1, be, 1, rwork, 1 )
832  CALL zlacpy( ' ', m, nrhs, y, ldx, z, ldx )
833  CALL zlaset( 'Full', mnmin, mnmin, czero, cone, u, ldpt )
834  CALL zlaset( 'Full', mnmin, mnmin, czero, cone, vt, ldpt )
835 *
836  CALL zbdsqr( uplo, mnmin, mnmin, mnmin, nrhs, s1, rwork, vt,
837  $ ldpt, u, ldpt, z, ldx, rwork( mnmin+1 ),
838  $ iinfo )
839 *
840 * Check error code from ZBDSQR.
841 *
842  IF( iinfo.NE.0 ) THEN
843  WRITE( nout, fmt = 9998 )'ZBDSQR(vects)', iinfo, m, n,
844  $ jtype, ioldsd
845  info = abs( iinfo )
846  IF( iinfo.LT.0 ) THEN
847  return
848  ELSE
849  result( 4 ) = ulpinv
850  go to 150
851  END IF
852  END IF
853 *
854 * Use ZBDSQR to compute only the singular values of the
855 * bidiagonal matrix B; U, VT, and Z should not be modified.
856 *
857  CALL dcopy( mnmin, bd, 1, s2, 1 )
858  IF( mnmin.GT.0 )
859  $ CALL dcopy( mnmin-1, be, 1, rwork, 1 )
860 *
861  CALL zbdsqr( uplo, mnmin, 0, 0, 0, s2, rwork, vt, ldpt, u,
862  $ ldpt, z, ldx, rwork( mnmin+1 ), iinfo )
863 *
864 * Check error code from ZBDSQR.
865 *
866  IF( iinfo.NE.0 ) THEN
867  WRITE( nout, fmt = 9998 )'ZBDSQR(values)', iinfo, m, n,
868  $ jtype, ioldsd
869  info = abs( iinfo )
870  IF( iinfo.LT.0 ) THEN
871  return
872  ELSE
873  result( 9 ) = ulpinv
874  go to 150
875  END IF
876  END IF
877 *
878 * Test 4: Check the decomposition B := U * S1 * VT
879 * 5: Check the computation Z := U' * Y
880 * 6: Check the orthogonality of U
881 * 7: Check the orthogonality of VT
882 *
883  CALL zbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
884  $ work, result( 4 ) )
885  CALL zbdt02( mnmin, nrhs, y, ldx, z, ldx, u, ldpt, work,
886  $ rwork, result( 5 ) )
887  CALL zunt01( 'Columns', mnmin, mnmin, u, ldpt, work, lwork,
888  $ rwork, result( 6 ) )
889  CALL zunt01( 'Rows', mnmin, mnmin, vt, ldpt, work, lwork,
890  $ rwork, result( 7 ) )
891 *
892 * Test 8: Check that the singular values are sorted in
893 * non-increasing order and are non-negative
894 *
895  result( 8 ) = zero
896  DO 110 i = 1, mnmin - 1
897  IF( s1( i ).LT.s1( i+1 ) )
898  $ result( 8 ) = ulpinv
899  IF( s1( i ).LT.zero )
900  $ result( 8 ) = ulpinv
901  110 continue
902  IF( mnmin.GE.1 ) THEN
903  IF( s1( mnmin ).LT.zero )
904  $ result( 8 ) = ulpinv
905  END IF
906 *
907 * Test 9: Compare ZBDSQR with and without singular vectors
908 *
909  temp2 = zero
910 *
911  DO 120 j = 1, mnmin
912  temp1 = abs( s1( j )-s2( j ) ) /
913  $ max( sqrt( unfl )*max( s1( 1 ), one ),
914  $ ulp*max( abs( s1( j ) ), abs( s2( j ) ) ) )
915  temp2 = max( temp1, temp2 )
916  120 continue
917 *
918  result( 9 ) = temp2
919 *
920 * Test 10: Sturm sequence test of singular values
921 * Go up by factors of two until it succeeds
922 *
923  temp1 = thresh*( half-ulp )
924 *
925  DO 130 j = 0, log2ui
926  CALL dsvdch( mnmin, bd, be, s1, temp1, iinfo )
927  IF( iinfo.EQ.0 )
928  $ go to 140
929  temp1 = temp1*two
930  130 continue
931 *
932  140 continue
933  result( 10 ) = temp1
934 *
935 * Use ZBDSQR to form the decomposition A := (QU) S (VT PT)
936 * from the bidiagonal form A := Q B PT.
937 *
938  IF( .NOT.bidiag ) THEN
939  CALL dcopy( mnmin, bd, 1, s2, 1 )
940  IF( mnmin.GT.0 )
941  $ CALL dcopy( mnmin-1, be, 1, rwork, 1 )
942 *
943  CALL zbdsqr( uplo, mnmin, n, m, nrhs, s2, rwork, pt,
944  $ ldpt, q, ldq, y, ldx, rwork( mnmin+1 ),
945  $ iinfo )
946 *
947 * Test 11: Check the decomposition A := Q*U * S2 * VT*PT
948 * 12: Check the computation Z := U' * Q' * X
949 * 13: Check the orthogonality of Q*U
950 * 14: Check the orthogonality of VT*PT
951 *
952  CALL zbdt01( m, n, 0, a, lda, q, ldq, s2, dumma, pt,
953  $ ldpt, work, rwork, result( 11 ) )
954  CALL zbdt02( m, nrhs, x, ldx, y, ldx, q, ldq, work,
955  $ rwork, result( 12 ) )
956  CALL zunt01( 'Columns', m, mq, q, ldq, work, lwork,
957  $ rwork, result( 13 ) )
958  CALL zunt01( 'Rows', mnmin, n, pt, ldpt, work, lwork,
959  $ rwork, result( 14 ) )
960  END IF
961 *
962 * End of Loop -- Check for RESULT(j) > THRESH
963 *
964  150 continue
965  DO 160 j = 1, 14
966  IF( result( j ).GE.thresh ) THEN
967  IF( nfail.EQ.0 )
968  $ CALL dlahd2( nout, path )
969  WRITE( nout, fmt = 9999 )m, n, jtype, ioldsd, j,
970  $ result( j )
971  nfail = nfail + 1
972  END IF
973  160 continue
974  IF( .NOT.bidiag ) THEN
975  ntest = ntest + 14
976  ELSE
977  ntest = ntest + 5
978  END IF
979 *
980  170 continue
981  180 continue
982 *
983 * Summary
984 *
985  CALL alasum( path, nout, nfail, ntest, 0 )
986 *
987  return
988 *
989 * End of ZCHKBD
990 *
991  9999 format( ' M=', i5, ', N=', i5, ', type ', i2, ', seed=',
992  $ 4( i4, ',' ), ' test(', i2, ')=', g11.4 )
993  9998 format( ' ZCHKBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
994  $ i6, ', N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
995  $ i5, ')' )
996 *
997  END