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