LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
schkbd.f
Go to the documentation of this file.
1 *> \brief \b SCHKBD
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 SCHKBD( 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 * REAL THRESH
20 * ..
21 * .. Array Arguments ..
22 * LOGICAL DOTYPE( * )
23 * INTEGER ISEED( 4 ), IWORK( * ), MVAL( * ), NVAL( * )
24 * REAL 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 *> SCHKBD checks the singular value decomposition (SVD) routines.
37 *>
38 *> SGEBRD 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 *> SORGBR generates the orthogonal matrices Q and P' from SGEBRD.
44 *> Note that Q and P are not necessarily square.
45 *>
46 *> SBDSQR 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, SBDSQR has an option to apply the left orthogonal matrix
55 *> U to a matrix X, useful in least squares applications.
56 *>
57 *> SBDSDC 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 *> SBDSVDX computes the singular value decomposition of the bidiagonal
67 *> matrix B as B = U S V' using bisection and inverse iteration. It is
68 *> called six times to compute
69 *> 1) B = U S1 V', RANGE='A', where S1 is the diagonal matrix of singular
70 *> values and the columns of the matrices U and V are the left
71 *> and right singular vectors, respectively, of B.
72 *> 2) Same as 1), but the singular values are stored in S2 and the
73 *> singular vectors are not computed.
74 *> 3) B = U S1 V', RANGE='I', with where S1 is the diagonal matrix of singular
75 *> values and the columns of the matrices U and V are the left
76 *> and right singular vectors, respectively, of B
77 *> 4) Same as 3), but the singular values are stored in S2 and the
78 *> singular vectors are not computed.
79 *> 5) B = U S1 V', RANGE='V', with where S1 is the diagonal matrix of singular
80 *> values and the columns of the matrices U and V are the left
81 *> and right singular vectors, respectively, of B
82 *> 6) Same as 5), but the singular values are stored in S2 and the
83 *> singular vectors are not computed.
84 *>
85 *> For each pair of matrix dimensions (M,N) and each selected matrix
86 *> type, an M by N matrix A and an M by NRHS matrix X are generated.
87 *> The problem dimensions are as follows
88 *> A: M x N
89 *> Q: M x min(M,N) (but M x M if NRHS > 0)
90 *> P: min(M,N) x N
91 *> B: min(M,N) x min(M,N)
92 *> U, V: min(M,N) x min(M,N)
93 *> S1, S2 diagonal, order min(M,N)
94 *> X: M x NRHS
95 *>
96 *> For each generated matrix, 14 tests are performed:
97 *>
98 *> Test SGEBRD and SORGBR
99 *>
100 *> (1) | A - Q B PT | / ( |A| max(M,N) ulp ), PT = P'
101 *>
102 *> (2) | I - Q' Q | / ( M ulp )
103 *>
104 *> (3) | I - PT PT' | / ( N ulp )
105 *>
106 *> Test SBDSQR on bidiagonal matrix B
107 *>
108 *> (4) | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V'
109 *>
110 *> (5) | Y - U Z | / ( |Y| max(min(M,N),k) ulp ), where Y = Q' X
111 *> and Z = U' Y.
112 *> (6) | I - U' U | / ( min(M,N) ulp )
113 *>
114 *> (7) | I - VT VT' | / ( min(M,N) ulp )
115 *>
116 *> (8) S1 contains min(M,N) nonnegative values in decreasing order.
117 *> (Return 0 if true, 1/ULP if false.)
118 *>
119 *> (9) | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
120 *> computing U and V.
121 *>
122 *> (10) 0 if the true singular values of B are within THRESH of
123 *> those in S1. 2*THRESH if they are not. (Tested using
124 *> SSVDCH)
125 *>
126 *> Test SBDSQR on matrix A
127 *>
128 *> (11) | A - (QU) S (VT PT) | / ( |A| max(M,N) ulp )
129 *>
130 *> (12) | X - (QU) Z | / ( |X| max(M,k) ulp )
131 *>
132 *> (13) | I - (QU)'(QU) | / ( M ulp )
133 *>
134 *> (14) | I - (VT PT) (PT'VT') | / ( N ulp )
135 *>
136 *> Test SBDSDC on bidiagonal matrix B
137 *>
138 *> (15) | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V'
139 *>
140 *> (16) | I - U' U | / ( min(M,N) ulp )
141 *>
142 *> (17) | I - VT VT' | / ( min(M,N) ulp )
143 *>
144 *> (18) S1 contains min(M,N) nonnegative values in decreasing order.
145 *> (Return 0 if true, 1/ULP if false.)
146 *>
147 *> (19) | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
148 *> computing U and V.
149 *> Test SBDSVDX on bidiagonal matrix B
150 *>
151 *> (20) | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V'
152 *>
153 *> (21) | I - U' U | / ( min(M,N) ulp )
154 *>
155 *> (22) | I - VT VT' | / ( min(M,N) ulp )
156 *>
157 *> (23) S1 contains min(M,N) nonnegative values in decreasing order.
158 *> (Return 0 if true, 1/ULP if false.)
159 *>
160 *> (24) | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
161 *> computing U and V.
162 *>
163 *> (25) | S1 - U' B VT' | / ( |S| n ulp ) SBDSVDX('V', 'I')
164 *>
165 *> (26) | I - U' U | / ( min(M,N) ulp )
166 *>
167 *> (27) | I - VT VT' | / ( min(M,N) ulp )
168 *>
169 *> (28) S1 contains min(M,N) nonnegative values in decreasing order.
170 *> (Return 0 if true, 1/ULP if false.)
171 *>
172 *> (29) | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
173 *> computing U and V.
174 *>
175 *> (30) | S1 - U' B VT' | / ( |S1| n ulp ) SBDSVDX('V', 'V')
176 *>
177 *> (31) | I - U' U | / ( min(M,N) ulp )
178 *>
179 *> (32) | I - VT VT' | / ( min(M,N) ulp )
180 *>
181 *> (33) S1 contains min(M,N) nonnegative values in decreasing order.
182 *> (Return 0 if true, 1/ULP if false.)
183 *>
184 *> (34) | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
185 *> computing U and V.
186 *>
187 *> The possible matrix types are
188 *>
189 *> (1) The zero matrix.
190 *> (2) The identity matrix.
191 *>
192 *> (3) A diagonal matrix with evenly spaced entries
193 *> 1, ..., ULP and random signs.
194 *> (ULP = (first number larger than 1) - 1 )
195 *> (4) A diagonal matrix with geometrically spaced entries
196 *> 1, ..., ULP and random signs.
197 *> (5) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
198 *> and random signs.
199 *>
200 *> (6) Same as (3), but multiplied by SQRT( overflow threshold )
201 *> (7) Same as (3), but multiplied by SQRT( underflow threshold )
202 *>
203 *> (8) A matrix of the form U D V, where U and V are orthogonal and
204 *> D has evenly spaced entries 1, ..., ULP with random signs
205 *> on the diagonal.
206 *>
207 *> (9) A matrix of the form U D V, where U and V are orthogonal and
208 *> D has geometrically spaced entries 1, ..., ULP with random
209 *> signs on the diagonal.
210 *>
211 *> (10) A matrix of the form U D V, where U and V are orthogonal and
212 *> D has "clustered" entries 1, ULP,..., ULP with random
213 *> signs on the diagonal.
214 *>
215 *> (11) Same as (8), but multiplied by SQRT( overflow threshold )
216 *> (12) Same as (8), but multiplied by SQRT( underflow threshold )
217 *>
218 *> (13) Rectangular matrix with random entries chosen from (-1,1).
219 *> (14) Same as (13), but multiplied by SQRT( overflow threshold )
220 *> (15) Same as (13), but multiplied by SQRT( underflow threshold )
221 *>
222 *> Special case:
223 *> (16) A bidiagonal matrix with random entries chosen from a
224 *> logarithmic distribution on [ulp^2,ulp^(-2)] (I.e., each
225 *> entry is e^x, where x is chosen uniformly on
226 *> [ 2 log(ulp), -2 log(ulp) ] .) For *this* type:
227 *> (a) SGEBRD is not called to reduce it to bidiagonal form.
228 *> (b) the bidiagonal is min(M,N) x min(M,N); if M<N, the
229 *> matrix will be lower bidiagonal, otherwise upper.
230 *> (c) only tests 5--8 and 14 are performed.
231 *>
232 *> A subset of the full set of matrix types may be selected through
233 *> the logical array DOTYPE.
234 *> \endverbatim
235 *
236 * Arguments:
237 * ==========
238 *
239 *> \param[in] NSIZES
240 *> \verbatim
241 *> NSIZES is INTEGER
242 *> The number of values of M and N contained in the vectors
243 *> MVAL and NVAL. The matrix sizes are used in pairs (M,N).
244 *> \endverbatim
245 *>
246 *> \param[in] MVAL
247 *> \verbatim
248 *> MVAL is INTEGER array, dimension (NM)
249 *> The values of the matrix row dimension M.
250 *> \endverbatim
251 *>
252 *> \param[in] NVAL
253 *> \verbatim
254 *> NVAL is INTEGER array, dimension (NM)
255 *> The values of the matrix column dimension N.
256 *> \endverbatim
257 *>
258 *> \param[in] NTYPES
259 *> \verbatim
260 *> NTYPES is INTEGER
261 *> The number of elements in DOTYPE. If it is zero, SCHKBD
262 *> does nothing. It must be at least zero. If it is MAXTYP+1
263 *> and NSIZES is 1, then an additional type, MAXTYP+1 is
264 *> defined, which is to use whatever matrices are in A and B.
265 *> This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
266 *> DOTYPE(MAXTYP+1) is .TRUE. .
267 *> \endverbatim
268 *>
269 *> \param[in] DOTYPE
270 *> \verbatim
271 *> DOTYPE is LOGICAL array, dimension (NTYPES)
272 *> If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix
273 *> of type j will be generated. If NTYPES is smaller than the
274 *> maximum number of types defined (PARAMETER MAXTYP), then
275 *> types NTYPES+1 through MAXTYP will not be generated. If
276 *> NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through
277 *> DOTYPE(NTYPES) will be ignored.
278 *> \endverbatim
279 *>
280 *> \param[in] NRHS
281 *> \verbatim
282 *> NRHS is INTEGER
283 *> The number of columns in the "right-hand side" matrices X, Y,
284 *> and Z, used in testing SBDSQR. If NRHS = 0, then the
285 *> operations on the right-hand side will not be tested.
286 *> NRHS must be at least 0.
287 *> \endverbatim
288 *>
289 *> \param[in,out] ISEED
290 *> \verbatim
291 *> ISEED is INTEGER array, dimension (4)
292 *> On entry ISEED specifies the seed of the random number
293 *> generator. The array elements should be between 0 and 4095;
294 *> if not they will be reduced mod 4096. Also, ISEED(4) must
295 *> be odd. The values of ISEED are changed on exit, and can be
296 *> used in the next call to SCHKBD to continue the same random
297 *> number sequence.
298 *> \endverbatim
299 *>
300 *> \param[in] THRESH
301 *> \verbatim
302 *> THRESH is REAL
303 *> The threshold value for the test ratios. A result is
304 *> included in the output file if RESULT >= THRESH. To have
305 *> every test ratio printed, use THRESH = 0. Note that the
306 *> expected value of the test ratios is O(1), so THRESH should
307 *> be a reasonably small multiple of 1, e.g., 10 or 100.
308 *> \endverbatim
309 *>
310 *> \param[out] A
311 *> \verbatim
312 *> A is REAL array, dimension (LDA,NMAX)
313 *> where NMAX is the maximum value of N in NVAL.
314 *> \endverbatim
315 *>
316 *> \param[in] LDA
317 *> \verbatim
318 *> LDA is INTEGER
319 *> The leading dimension of the array A. LDA >= max(1,MMAX),
320 *> where MMAX is the maximum value of M in MVAL.
321 *> \endverbatim
322 *>
323 *> \param[out] BD
324 *> \verbatim
325 *> BD is REAL array, dimension
326 *> (max(min(MVAL(j),NVAL(j))))
327 *> \endverbatim
328 *>
329 *> \param[out] BE
330 *> \verbatim
331 *> BE is REAL array, dimension
332 *> (max(min(MVAL(j),NVAL(j))))
333 *> \endverbatim
334 *>
335 *> \param[out] S1
336 *> \verbatim
337 *> S1 is REAL array, dimension
338 *> (max(min(MVAL(j),NVAL(j))))
339 *> \endverbatim
340 *>
341 *> \param[out] S2
342 *> \verbatim
343 *> S2 is REAL array, dimension
344 *> (max(min(MVAL(j),NVAL(j))))
345 *> \endverbatim
346 *>
347 *> \param[out] X
348 *> \verbatim
349 *> X is REAL array, dimension (LDX,NRHS)
350 *> \endverbatim
351 *>
352 *> \param[in] LDX
353 *> \verbatim
354 *> LDX is INTEGER
355 *> The leading dimension of the arrays X, Y, and Z.
356 *> LDX >= max(1,MMAX)
357 *> \endverbatim
358 *>
359 *> \param[out] Y
360 *> \verbatim
361 *> Y is REAL array, dimension (LDX,NRHS)
362 *> \endverbatim
363 *>
364 *> \param[out] Z
365 *> \verbatim
366 *> Z is REAL array, dimension (LDX,NRHS)
367 *> \endverbatim
368 *>
369 *> \param[out] Q
370 *> \verbatim
371 *> Q is REAL array, dimension (LDQ,MMAX)
372 *> \endverbatim
373 *>
374 *> \param[in] LDQ
375 *> \verbatim
376 *> LDQ is INTEGER
377 *> The leading dimension of the array Q. LDQ >= max(1,MMAX).
378 *> \endverbatim
379 *>
380 *> \param[out] PT
381 *> \verbatim
382 *> PT is REAL array, dimension (LDPT,NMAX)
383 *> \endverbatim
384 *>
385 *> \param[in] LDPT
386 *> \verbatim
387 *> LDPT is INTEGER
388 *> The leading dimension of the arrays PT, U, and V.
389 *> LDPT >= max(1, max(min(MVAL(j),NVAL(j)))).
390 *> \endverbatim
391 *>
392 *> \param[out] U
393 *> \verbatim
394 *> U is REAL array, dimension
395 *> (LDPT,max(min(MVAL(j),NVAL(j))))
396 *> \endverbatim
397 *>
398 *> \param[out] VT
399 *> \verbatim
400 *> VT is REAL array, dimension
401 *> (LDPT,max(min(MVAL(j),NVAL(j))))
402 *> \endverbatim
403 *>
404 *> \param[out] WORK
405 *> \verbatim
406 *> WORK is REAL array, dimension (LWORK)
407 *> \endverbatim
408 *>
409 *> \param[in] LWORK
410 *> \verbatim
411 *> LWORK is INTEGER
412 *> The number of entries in WORK. This must be at least
413 *> 3(M+N) and M(M + max(M,N,k) + 1) + N*min(M,N) for all
414 *> pairs (M,N)=(MM(j),NN(j))
415 *> \endverbatim
416 *>
417 *> \param[out] IWORK
418 *> \verbatim
419 *> IWORK is INTEGER array, dimension at least 8*min(M,N)
420 *> \endverbatim
421 *>
422 *> \param[in] NOUT
423 *> \verbatim
424 *> NOUT is INTEGER
425 *> The FORTRAN unit number for printing out error messages
426 *> (e.g., if a routine returns IINFO not equal to 0.)
427 *> \endverbatim
428 *>
429 *> \param[out] INFO
430 *> \verbatim
431 *> INFO is INTEGER
432 *> If 0, then everything ran OK.
433 *> -1: NSIZES < 0
434 *> -2: Some MM(j) < 0
435 *> -3: Some NN(j) < 0
436 *> -4: NTYPES < 0
437 *> -6: NRHS < 0
438 *> -8: THRESH < 0
439 *> -11: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ).
440 *> -17: LDB < 1 or LDB < MMAX.
441 *> -21: LDQ < 1 or LDQ < MMAX.
442 *> -23: LDPT< 1 or LDPT< MNMAX.
443 *> -27: LWORK too small.
444 *> If SLATMR, SLATMS, SGEBRD, SORGBR, or SBDSQR,
445 *> returns an error code, the
446 *> absolute value of it is returned.
447 *>
448 *>-----------------------------------------------------------------------
449 *>
450 *> Some Local Variables and Parameters:
451 *> ---- ----- --------- --- ----------
452 *>
453 *> ZERO, ONE Real 0 and 1.
454 *> MAXTYP The number of types defined.
455 *> NTEST The number of tests performed, or which can
456 *> be performed so far, for the current matrix.
457 *> MMAX Largest value in NN.
458 *> NMAX Largest value in NN.
459 *> MNMIN min(MM(j), NN(j)) (the dimension of the bidiagonal
460 *> matrix.)
461 *> MNMAX The maximum value of MNMIN for j=1,...,NSIZES.
462 *> NFAIL The number of tests which have exceeded THRESH
463 *> COND, IMODE Values to be passed to the matrix generators.
464 *> ANORM Norm of A; passed to matrix generators.
465 *>
466 *> OVFL, UNFL Overflow and underflow thresholds.
467 *> RTOVFL, RTUNFL Square roots of the previous 2 values.
468 *> ULP, ULPINV Finest relative precision and its inverse.
469 *>
470 *> The following four arrays decode JTYPE:
471 *> KTYPE(j) The general type (1-10) for type "j".
472 *> KMODE(j) The MODE value to be passed to the matrix
473 *> generator for type "j".
474 *> KMAGN(j) The order of magnitude ( O(1),
475 *> O(overflow^(1/2) ), O(underflow^(1/2) )
476 *> \endverbatim
477 *
478 * Authors:
479 * ========
480 *
481 *> \author Univ. of Tennessee
482 *> \author Univ. of California Berkeley
483 *> \author Univ. of Colorado Denver
484 *> \author NAG Ltd.
485 *
486 *> \date June 2016
487 *
488 *> \ingroup single_eig
489 *
490 * =====================================================================
491  SUBROUTINE schkbd( NSIZES, MVAL, NVAL, NTYPES, DOTYPE, NRHS,
492  $ iseed, thresh, a, lda, bd, be, s1, s2, x, ldx,
493  $ y, z, q, ldq, pt, ldpt, u, vt, work, lwork,
494  $ iwork, nout, info )
495 *
496 * -- LAPACK test routine (version 3.6.1) --
497 * -- LAPACK is a software package provided by Univ. of Tennessee, --
498 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
499 * June 2016
500 *
501 * .. Scalar Arguments ..
502  INTEGER INFO, LDA, LDPT, LDQ, LDX, LWORK, NOUT, NRHS,
503  $ nsizes, ntypes
504  REAL THRESH
505 * ..
506 * .. Array Arguments ..
507  LOGICAL DOTYPE( * )
508  INTEGER ISEED( 4 ), IWORK( * ), MVAL( * ), NVAL( * )
509  REAL A( lda, * ), BD( * ), BE( * ), PT( ldpt, * ),
510  $ q( ldq, * ), s1( * ), s2( * ), u( ldpt, * ),
511  $ vt( ldpt, * ), work( * ), x( ldx, * ),
512  $ y( ldx, * ), z( ldx, * )
513 * ..
514 *
515 * ======================================================================
516 *
517 * .. Parameters ..
518  REAL ZERO, ONE, TWO, HALF
519  parameter ( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
520  $ half = 0.5e0 )
521  INTEGER MAXTYP
522  parameter ( maxtyp = 16 )
523 * ..
524 * .. Local Scalars ..
525  LOGICAL BADMM, BADNN, BIDIAG
526  CHARACTER UPLO
527  CHARACTER*3 PATH
528  INTEGER I, IINFO, IL, IMODE, ITEMP, ITYPE, IU, IWBD,
529  $ iwbe, iwbs, iwbz, iwwork, j, jcol, jsize,
530  $ jtype, log2ui, m, minwrk, mmax, mnmax, mnmin,
531  $ mnmin2, mq, mtypes, n, nfail, nmax,
532  $ ns1, ns2, ntest
533  REAL ABSTOL, AMNINV, ANORM, COND, OVFL, RTOVFL,
534  $ rtunfl, temp1, temp2, temp3, ulp, ulpinv,
535  $ unfl, vl, vu
536 * ..
537 * .. Local Arrays ..
538  INTEGER IDUM( 1 ), IOLDSD( 4 ), ISEED2( 4 ),
539  $ kmagn( maxtyp ), kmode( maxtyp ),
540  $ ktype( maxtyp )
541  REAL DUM( 1 ), DUMMA( 1 ), RESULT( 40 )
542 * ..
543 * .. External Functions ..
544  REAL SLAMCH, SLARND, SSXT1
545  EXTERNAL slamch, slarnd, ssxt1
546 * ..
547 * .. External Subroutines ..
548  EXTERNAL alasum, sbdsdc, sbdsqr, sbdsvdx, sbdt01,
552 * ..
553 * .. Intrinsic Functions ..
554  INTRINSIC abs, exp, int, log, max, min, sqrt
555 * ..
556 * .. Scalars in Common ..
557  LOGICAL LERR, OK
558  CHARACTER*32 SRNAMT
559  INTEGER INFOT, NUNIT
560 * ..
561 * .. Common blocks ..
562  COMMON / infoc / infot, nunit, ok, lerr
563  COMMON / srnamc / srnamt
564 * ..
565 * .. Data statements ..
566  DATA ktype / 1, 2, 5*4, 5*6, 3*9, 10 /
567  DATA kmagn / 2*1, 3*1, 2, 3, 3*1, 2, 3, 1, 2, 3, 0 /
568  DATA kmode / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
569  $ 0, 0, 0 /
570 * ..
571 * .. Executable Statements ..
572 *
573 * Check for errors
574 *
575  info = 0
576 *
577  badmm = .false.
578  badnn = .false.
579  mmax = 1
580  nmax = 1
581  mnmax = 1
582  minwrk = 1
583  DO 10 j = 1, nsizes
584  mmax = max( mmax, mval( j ) )
585  IF( mval( j ).LT.0 )
586  $ badmm = .true.
587  nmax = max( nmax, nval( j ) )
588  IF( nval( j ).LT.0 )
589  $ badnn = .true.
590  mnmax = max( mnmax, min( mval( j ), nval( j ) ) )
591  minwrk = max( minwrk, 3*( mval( j )+nval( j ) ),
592  $ mval( j )*( mval( j )+max( mval( j ), nval( j ),
593  $ nrhs )+1 )+nval( j )*min( nval( j ), mval( j ) ) )
594  10 CONTINUE
595 *
596 * Check for errors
597 *
598  IF( nsizes.LT.0 ) THEN
599  info = -1
600  ELSE IF( badmm ) THEN
601  info = -2
602  ELSE IF( badnn ) THEN
603  info = -3
604  ELSE IF( ntypes.LT.0 ) THEN
605  info = -4
606  ELSE IF( nrhs.LT.0 ) THEN
607  info = -6
608  ELSE IF( lda.LT.mmax ) THEN
609  info = -11
610  ELSE IF( ldx.LT.mmax ) THEN
611  info = -17
612  ELSE IF( ldq.LT.mmax ) THEN
613  info = -21
614  ELSE IF( ldpt.LT.mnmax ) THEN
615  info = -23
616  ELSE IF( minwrk.GT.lwork ) THEN
617  info = -27
618  END IF
619 *
620  IF( info.NE.0 ) THEN
621  CALL xerbla( 'SCHKBD', -info )
622  RETURN
623  END IF
624 *
625 * Initialize constants
626 *
627  path( 1: 1 ) = 'Single precision'
628  path( 2: 3 ) = 'BD'
629  nfail = 0
630  ntest = 0
631  unfl = slamch( 'Safe minimum' )
632  ovfl = slamch( 'Overflow' )
633  CALL slabad( unfl, ovfl )
634  ulp = slamch( 'Precision' )
635  ulpinv = one / ulp
636  log2ui = int( log( ulpinv ) / log( two ) )
637  rtunfl = sqrt( unfl )
638  rtovfl = sqrt( ovfl )
639  infot = 0
640  abstol = 2*unfl
641 *
642 * Loop over sizes, types
643 *
644  DO 300 jsize = 1, nsizes
645  m = mval( jsize )
646  n = nval( jsize )
647  mnmin = min( m, n )
648  amninv = one / max( m, n, 1 )
649 *
650  IF( nsizes.NE.1 ) THEN
651  mtypes = min( maxtyp, ntypes )
652  ELSE
653  mtypes = min( maxtyp+1, ntypes )
654  END IF
655 *
656  DO 290 jtype = 1, mtypes
657  IF( .NOT.dotype( jtype ) )
658  $ GO TO 290
659 *
660  DO 20 j = 1, 4
661  ioldsd( j ) = iseed( j )
662  20 CONTINUE
663 *
664  DO 30 j = 1, 34
665  result( j ) = -one
666  30 CONTINUE
667 *
668  uplo = ' '
669 *
670 * Compute "A"
671 *
672 * Control parameters:
673 *
674 * KMAGN KMODE KTYPE
675 * =1 O(1) clustered 1 zero
676 * =2 large clustered 2 identity
677 * =3 small exponential (none)
678 * =4 arithmetic diagonal, (w/ eigenvalues)
679 * =5 random symmetric, w/ eigenvalues
680 * =6 nonsymmetric, w/ singular values
681 * =7 random diagonal
682 * =8 random symmetric
683 * =9 random nonsymmetric
684 * =10 random bidiagonal (log. distrib.)
685 *
686  IF( mtypes.GT.maxtyp )
687  $ GO TO 100
688 *
689  itype = ktype( jtype )
690  imode = kmode( jtype )
691 *
692 * Compute norm
693 *
694  GO TO ( 40, 50, 60 )kmagn( jtype )
695 *
696  40 CONTINUE
697  anorm = one
698  GO TO 70
699 *
700  50 CONTINUE
701  anorm = ( rtovfl*ulp )*amninv
702  GO TO 70
703 *
704  60 CONTINUE
705  anorm = rtunfl*max( m, n )*ulpinv
706  GO TO 70
707 *
708  70 CONTINUE
709 *
710  CALL slaset( 'Full', lda, n, zero, zero, a, lda )
711  iinfo = 0
712  cond = ulpinv
713 *
714  bidiag = .false.
715  IF( itype.EQ.1 ) THEN
716 *
717 * Zero matrix
718 *
719  iinfo = 0
720 *
721  ELSE IF( itype.EQ.2 ) THEN
722 *
723 * Identity
724 *
725  DO 80 jcol = 1, mnmin
726  a( jcol, jcol ) = anorm
727  80 CONTINUE
728 *
729  ELSE IF( itype.EQ.4 ) THEN
730 *
731 * Diagonal Matrix, [Eigen]values Specified
732 *
733  CALL slatms( mnmin, mnmin, 'S', iseed, 'N', work, imode,
734  $ cond, anorm, 0, 0, 'N', a, lda,
735  $ work( mnmin+1 ), iinfo )
736 *
737  ELSE IF( itype.EQ.5 ) THEN
738 *
739 * Symmetric, eigenvalues specified
740 *
741  CALL slatms( mnmin, mnmin, 'S', iseed, 'S', work, imode,
742  $ cond, anorm, m, n, 'N', a, lda,
743  $ work( mnmin+1 ), iinfo )
744 *
745  ELSE IF( itype.EQ.6 ) THEN
746 *
747 * Nonsymmetric, singular values specified
748 *
749  CALL slatms( m, n, 'S', iseed, 'N', work, imode, cond,
750  $ anorm, m, n, 'N', a, lda, work( mnmin+1 ),
751  $ iinfo )
752 *
753  ELSE IF( itype.EQ.7 ) THEN
754 *
755 * Diagonal, random entries
756 *
757  CALL slatmr( mnmin, mnmin, 'S', iseed, 'N', work, 6, one,
758  $ one, 'T', 'N', work( mnmin+1 ), 1, one,
759  $ work( 2*mnmin+1 ), 1, one, 'N', iwork, 0, 0,
760  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
761 *
762  ELSE IF( itype.EQ.8 ) THEN
763 *
764 * Symmetric, random entries
765 *
766  CALL slatmr( mnmin, mnmin, 'S', iseed, 'S', work, 6, one,
767  $ one, 'T', 'N', work( mnmin+1 ), 1, one,
768  $ work( m+mnmin+1 ), 1, one, 'N', iwork, m, n,
769  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
770 *
771  ELSE IF( itype.EQ.9 ) THEN
772 *
773 * Nonsymmetric, random entries
774 *
775  CALL slatmr( m, n, 'S', iseed, 'N', work, 6, one, one,
776  $ 'T', 'N', work( mnmin+1 ), 1, one,
777  $ work( m+mnmin+1 ), 1, one, 'N', iwork, m, n,
778  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
779 *
780  ELSE IF( itype.EQ.10 ) THEN
781 *
782 * Bidiagonal, random entries
783 *
784  temp1 = -two*log( ulp )
785  DO 90 j = 1, mnmin
786  bd( j ) = exp( temp1*slarnd( 2, iseed ) )
787  IF( j.LT.mnmin )
788  $ be( j ) = exp( temp1*slarnd( 2, iseed ) )
789  90 CONTINUE
790 *
791  iinfo = 0
792  bidiag = .true.
793  IF( m.GE.n ) THEN
794  uplo = 'U'
795  ELSE
796  uplo = 'L'
797  END IF
798  ELSE
799  iinfo = 1
800  END IF
801 *
802  IF( iinfo.EQ.0 ) THEN
803 *
804 * Generate Right-Hand Side
805 *
806  IF( bidiag ) THEN
807  CALL slatmr( mnmin, nrhs, 'S', iseed, 'N', work, 6,
808  $ one, one, 'T', 'N', work( mnmin+1 ), 1,
809  $ one, work( 2*mnmin+1 ), 1, one, 'N',
810  $ iwork, mnmin, nrhs, zero, one, 'NO', y,
811  $ ldx, iwork, iinfo )
812  ELSE
813  CALL slatmr( m, nrhs, 'S', iseed, 'N', work, 6, one,
814  $ one, 'T', 'N', work( m+1 ), 1, one,
815  $ work( 2*m+1 ), 1, one, 'N', iwork, m,
816  $ nrhs, zero, one, 'NO', x, ldx, iwork,
817  $ iinfo )
818  END IF
819  END IF
820 *
821 * Error Exit
822 *
823  IF( iinfo.NE.0 ) THEN
824  WRITE( nout, fmt = 9998 )'Generator', iinfo, m, n,
825  $ jtype, ioldsd
826  info = abs( iinfo )
827  RETURN
828  END IF
829 *
830  100 CONTINUE
831 *
832 * Call SGEBRD and SORGBR to compute B, Q, and P, do tests.
833 *
834  IF( .NOT.bidiag ) THEN
835 *
836 * Compute transformations to reduce A to bidiagonal form:
837 * B := Q' * A * P.
838 *
839  CALL slacpy( ' ', m, n, a, lda, q, ldq )
840  CALL sgebrd( m, n, q, ldq, bd, be, work, work( mnmin+1 ),
841  $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
842 *
843 * Check error code from SGEBRD.
844 *
845  IF( iinfo.NE.0 ) THEN
846  WRITE( nout, fmt = 9998 )'SGEBRD', iinfo, m, n,
847  $ jtype, ioldsd
848  info = abs( iinfo )
849  RETURN
850  END IF
851 *
852  CALL slacpy( ' ', m, n, q, ldq, pt, ldpt )
853  IF( m.GE.n ) THEN
854  uplo = 'U'
855  ELSE
856  uplo = 'L'
857  END IF
858 *
859 * Generate Q
860 *
861  mq = m
862  IF( nrhs.LE.0 )
863  $ mq = mnmin
864  CALL sorgbr( 'Q', m, mq, n, q, ldq, work,
865  $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
866 *
867 * Check error code from SORGBR.
868 *
869  IF( iinfo.NE.0 ) THEN
870  WRITE( nout, fmt = 9998 )'SORGBR(Q)', iinfo, m, n,
871  $ jtype, ioldsd
872  info = abs( iinfo )
873  RETURN
874  END IF
875 *
876 * Generate P'
877 *
878  CALL sorgbr( 'P', mnmin, n, m, pt, ldpt, work( mnmin+1 ),
879  $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
880 *
881 * Check error code from SORGBR.
882 *
883  IF( iinfo.NE.0 ) THEN
884  WRITE( nout, fmt = 9998 )'SORGBR(P)', iinfo, m, n,
885  $ jtype, ioldsd
886  info = abs( iinfo )
887  RETURN
888  END IF
889 *
890 * Apply Q' to an M by NRHS matrix X: Y := Q' * X.
891 *
892  CALL sgemm( 'Transpose', 'No transpose', m, nrhs, m, one,
893  $ q, ldq, x, ldx, zero, y, ldx )
894 *
895 * Test 1: Check the decomposition A := Q * B * PT
896 * 2: Check the orthogonality of Q
897 * 3: Check the orthogonality of PT
898 *
899  CALL sbdt01( m, n, 1, a, lda, q, ldq, bd, be, pt, ldpt,
900  $ work, result( 1 ) )
901  CALL sort01( 'Columns', m, mq, q, ldq, work, lwork,
902  $ result( 2 ) )
903  CALL sort01( 'Rows', mnmin, n, pt, ldpt, work, lwork,
904  $ result( 3 ) )
905  END IF
906 *
907 * Use SBDSQR to form the SVD of the bidiagonal matrix B:
908 * B := U * S1 * VT, and compute Z = U' * Y.
909 *
910  CALL scopy( mnmin, bd, 1, s1, 1 )
911  IF( mnmin.GT.0 )
912  $ CALL scopy( mnmin-1, be, 1, work, 1 )
913  CALL slacpy( ' ', m, nrhs, y, ldx, z, ldx )
914  CALL slaset( 'Full', mnmin, mnmin, zero, one, u, ldpt )
915  CALL slaset( 'Full', mnmin, mnmin, zero, one, vt, ldpt )
916 *
917  CALL sbdsqr( uplo, mnmin, mnmin, mnmin, nrhs, s1, work, vt,
918  $ ldpt, u, ldpt, z, ldx, work( mnmin+1 ), iinfo )
919 *
920 * Check error code from SBDSQR.
921 *
922  IF( iinfo.NE.0 ) THEN
923  WRITE( nout, fmt = 9998 )'SBDSQR(vects)', iinfo, m, n,
924  $ jtype, ioldsd
925  info = abs( iinfo )
926  IF( iinfo.LT.0 ) THEN
927  RETURN
928  ELSE
929  result( 4 ) = ulpinv
930  GO TO 270
931  END IF
932  END IF
933 *
934 * Use SBDSQR to compute only the singular values of the
935 * bidiagonal matrix B; U, VT, and Z should not be modified.
936 *
937  CALL scopy( mnmin, bd, 1, s2, 1 )
938  IF( mnmin.GT.0 )
939  $ CALL scopy( mnmin-1, be, 1, work, 1 )
940 *
941  CALL sbdsqr( uplo, mnmin, 0, 0, 0, s2, work, vt, ldpt, u,
942  $ ldpt, z, ldx, work( mnmin+1 ), iinfo )
943 *
944 * Check error code from SBDSQR.
945 *
946  IF( iinfo.NE.0 ) THEN
947  WRITE( nout, fmt = 9998 )'SBDSQR(values)', iinfo, m, n,
948  $ jtype, ioldsd
949  info = abs( iinfo )
950  IF( iinfo.LT.0 ) THEN
951  RETURN
952  ELSE
953  result( 9 ) = ulpinv
954  GO TO 270
955  END IF
956  END IF
957 *
958 * Test 4: Check the decomposition B := U * S1 * VT
959 * 5: Check the computation Z := U' * Y
960 * 6: Check the orthogonality of U
961 * 7: Check the orthogonality of VT
962 *
963  CALL sbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
964  $ work, result( 4 ) )
965  CALL sbdt02( mnmin, nrhs, y, ldx, z, ldx, u, ldpt, work,
966  $ result( 5 ) )
967  CALL sort01( 'Columns', mnmin, mnmin, u, ldpt, work, lwork,
968  $ result( 6 ) )
969  CALL sort01( 'Rows', mnmin, mnmin, vt, ldpt, work, lwork,
970  $ result( 7 ) )
971 *
972 * Test 8: Check that the singular values are sorted in
973 * non-increasing order and are non-negative
974 *
975  result( 8 ) = zero
976  DO 110 i = 1, mnmin - 1
977  IF( s1( i ).LT.s1( i+1 ) )
978  $ result( 8 ) = ulpinv
979  IF( s1( i ).LT.zero )
980  $ result( 8 ) = ulpinv
981  110 CONTINUE
982  IF( mnmin.GE.1 ) THEN
983  IF( s1( mnmin ).LT.zero )
984  $ result( 8 ) = ulpinv
985  END IF
986 *
987 * Test 9: Compare SBDSQR with and without singular vectors
988 *
989  temp2 = zero
990 *
991  DO 120 j = 1, mnmin
992  temp1 = abs( s1( j )-s2( j ) ) /
993  $ max( sqrt( unfl )*max( s1( 1 ), one ),
994  $ ulp*max( abs( s1( j ) ), abs( s2( j ) ) ) )
995  temp2 = max( temp1, temp2 )
996  120 CONTINUE
997 *
998  result( 9 ) = temp2
999 *
1000 * Test 10: Sturm sequence test of singular values
1001 * Go up by factors of two until it succeeds
1002 *
1003  temp1 = thresh*( half-ulp )
1004 *
1005  DO 130 j = 0, log2ui
1006 * CALL SSVDCH( MNMIN, BD, BE, S1, TEMP1, IINFO )
1007  IF( iinfo.EQ.0 )
1008  $ GO TO 140
1009  temp1 = temp1*two
1010  130 CONTINUE
1011 *
1012  140 CONTINUE
1013  result( 10 ) = temp1
1014 *
1015 * Use SBDSQR to form the decomposition A := (QU) S (VT PT)
1016 * from the bidiagonal form A := Q B PT.
1017 *
1018  IF( .NOT.bidiag ) THEN
1019  CALL scopy( mnmin, bd, 1, s2, 1 )
1020  IF( mnmin.GT.0 )
1021  $ CALL scopy( mnmin-1, be, 1, work, 1 )
1022 *
1023  CALL sbdsqr( uplo, mnmin, n, m, nrhs, s2, work, pt, ldpt,
1024  $ q, ldq, y, ldx, work( mnmin+1 ), iinfo )
1025 *
1026 * Test 11: Check the decomposition A := Q*U * S2 * VT*PT
1027 * 12: Check the computation Z := U' * Q' * X
1028 * 13: Check the orthogonality of Q*U
1029 * 14: Check the orthogonality of VT*PT
1030 *
1031  CALL sbdt01( m, n, 0, a, lda, q, ldq, s2, dumma, pt,
1032  $ ldpt, work, result( 11 ) )
1033  CALL sbdt02( m, nrhs, x, ldx, y, ldx, q, ldq, work,
1034  $ result( 12 ) )
1035  CALL sort01( 'Columns', m, mq, q, ldq, work, lwork,
1036  $ result( 13 ) )
1037  CALL sort01( 'Rows', mnmin, n, pt, ldpt, work, lwork,
1038  $ result( 14 ) )
1039  END IF
1040 *
1041 * Use SBDSDC to form the SVD of the bidiagonal matrix B:
1042 * B := U * S1 * VT
1043 *
1044  CALL scopy( mnmin, bd, 1, s1, 1 )
1045  IF( mnmin.GT.0 )
1046  $ CALL scopy( mnmin-1, be, 1, work, 1 )
1047  CALL slaset( 'Full', mnmin, mnmin, zero, one, u, ldpt )
1048  CALL slaset( 'Full', mnmin, mnmin, zero, one, vt, ldpt )
1049 *
1050  CALL sbdsdc( uplo, 'I', mnmin, s1, work, u, ldpt, vt, ldpt,
1051  $ dum, idum, work( mnmin+1 ), iwork, iinfo )
1052 *
1053 * Check error code from SBDSDC.
1054 *
1055  IF( iinfo.NE.0 ) THEN
1056  WRITE( nout, fmt = 9998 )'SBDSDC(vects)', iinfo, m, n,
1057  $ jtype, ioldsd
1058  info = abs( iinfo )
1059  IF( iinfo.LT.0 ) THEN
1060  RETURN
1061  ELSE
1062  result( 15 ) = ulpinv
1063  GO TO 270
1064  END IF
1065  END IF
1066 *
1067 * Use SBDSDC to compute only the singular values of the
1068 * bidiagonal matrix B; U and VT should not be modified.
1069 *
1070  CALL scopy( mnmin, bd, 1, s2, 1 )
1071  IF( mnmin.GT.0 )
1072  $ CALL scopy( mnmin-1, be, 1, work, 1 )
1073 *
1074  CALL sbdsdc( uplo, 'N', mnmin, s2, work, dum, 1, dum, 1,
1075  $ dum, idum, work( mnmin+1 ), iwork, iinfo )
1076 *
1077 * Check error code from SBDSDC.
1078 *
1079  IF( iinfo.NE.0 ) THEN
1080  WRITE( nout, fmt = 9998 )'SBDSDC(values)', iinfo, m, n,
1081  $ jtype, ioldsd
1082  info = abs( iinfo )
1083  IF( iinfo.LT.0 ) THEN
1084  RETURN
1085  ELSE
1086  result( 18 ) = ulpinv
1087  GO TO 270
1088  END IF
1089  END IF
1090 *
1091 * Test 15: Check the decomposition B := U * S1 * VT
1092 * 16: Check the orthogonality of U
1093 * 17: Check the orthogonality of VT
1094 *
1095  CALL sbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
1096  $ work, result( 15 ) )
1097  CALL sort01( 'Columns', mnmin, mnmin, u, ldpt, work, lwork,
1098  $ result( 16 ) )
1099  CALL sort01( 'Rows', mnmin, mnmin, vt, ldpt, work, lwork,
1100  $ result( 17 ) )
1101 *
1102 * Test 18: Check that the singular values are sorted in
1103 * non-increasing order and are non-negative
1104 *
1105  result( 18 ) = zero
1106  DO 150 i = 1, mnmin - 1
1107  IF( s1( i ).LT.s1( i+1 ) )
1108  $ result( 18 ) = ulpinv
1109  IF( s1( i ).LT.zero )
1110  $ result( 18 ) = ulpinv
1111  150 CONTINUE
1112  IF( mnmin.GE.1 ) THEN
1113  IF( s1( mnmin ).LT.zero )
1114  $ result( 18 ) = ulpinv
1115  END IF
1116 *
1117 * Test 19: Compare SBDSQR with and without singular vectors
1118 *
1119  temp2 = zero
1120 *
1121  DO 160 j = 1, mnmin
1122  temp1 = abs( s1( j )-s2( j ) ) /
1123  $ max( sqrt( unfl )*max( s1( 1 ), one ),
1124  $ ulp*max( abs( s1( 1 ) ), abs( s2( 1 ) ) ) )
1125  temp2 = max( temp1, temp2 )
1126  160 CONTINUE
1127 *
1128  result( 19 ) = temp2
1129 *
1130 *
1131 * Use SBDSVDX to compute the SVD of the bidiagonal matrix B:
1132 * B := U * S1 * VT
1133 *
1134  IF( jtype.EQ.10 .OR. jtype.EQ.16 ) THEN
1135 * =================================
1136 * Matrix types temporarily disabled
1137 * =================================
1138  result( 20:34 ) = zero
1139  GO TO 270
1140  END IF
1141 *
1142  iwbs = 1
1143  iwbd = iwbs + mnmin
1144  iwbe = iwbd + mnmin
1145  iwbz = iwbe + mnmin
1146  iwwork = iwbz + 2*mnmin*(mnmin+1)
1147  mnmin2 = max( 1,mnmin*2 )
1148 *
1149  CALL scopy( mnmin, bd, 1, work( iwbd ), 1 )
1150  IF( mnmin.GT.0 )
1151  $ CALL scopy( mnmin-1, be, 1, work( iwbe ), 1 )
1152 *
1153  CALL sbdsvdx( uplo, 'V', 'A', mnmin, work( iwbd ),
1154  $ work( iwbe ), zero, zero, 0, 0, ns1, s1,
1155  $ work( iwbz ), mnmin2, work( iwwork ),
1156  $ iwork, iinfo)
1157 *
1158 * Check error code from SBDSVDX.
1159 *
1160  IF( iinfo.NE.0 ) THEN
1161  WRITE( nout, fmt = 9998 )'SBDSVDX(vects,A)', iinfo, m, n,
1162  $ jtype, ioldsd
1163  info = abs( iinfo )
1164  IF( iinfo.LT.0 ) THEN
1165  RETURN
1166  ELSE
1167  result( 20 ) = ulpinv
1168  GO TO 270
1169  END IF
1170  END IF
1171 *
1172  j = iwbz
1173  DO 170 i = 1, ns1
1174  CALL scopy( mnmin, work( j ), 1, u( 1,i ), 1 )
1175  j = j + mnmin
1176  CALL scopy( mnmin, work( j ), 1, vt( i,1 ), ldpt )
1177  j = j + mnmin
1178  170 CONTINUE
1179 *
1180 * Use SBDSVDX to compute only the singular values of the
1181 * bidiagonal matrix B; U and VT should not be modified.
1182 *
1183  IF( jtype.EQ.9 ) THEN
1184 * =================================
1185 * Matrix types temporarily disabled
1186 * =================================
1187  result( 24 ) = zero
1188  GO TO 270
1189  END IF
1190 *
1191  CALL scopy( mnmin, bd, 1, work( iwbd ), 1 )
1192  IF( mnmin.GT.0 )
1193  $ CALL scopy( mnmin-1, be, 1, work( iwbe ), 1 )
1194 *
1195  CALL sbdsvdx( uplo, 'N', 'A', mnmin, work( iwbd ),
1196  $ work( iwbe ), zero, zero, 0, 0, ns2, s2,
1197  $ work( iwbz ), mnmin2, work( iwwork ),
1198  $ iwork, iinfo )
1199 *
1200 * Check error code from SBDSVDX.
1201 *
1202  IF( iinfo.NE.0 ) THEN
1203  WRITE( nout, fmt = 9998 )'SBDSVDX(values,A)', iinfo,
1204  $ m, n, jtype, ioldsd
1205  info = abs( iinfo )
1206  IF( iinfo.LT.0 ) THEN
1207  RETURN
1208  ELSE
1209  result( 24 ) = ulpinv
1210  GO TO 270
1211  END IF
1212  END IF
1213 *
1214 * Save S1 for tests 30-34.
1215 *
1216  CALL scopy( mnmin, s1, 1, work( iwbs ), 1 )
1217 *
1218 * Test 20: Check the decomposition B := U * S1 * VT
1219 * 21: Check the orthogonality of U
1220 * 22: Check the orthogonality of VT
1221 * 23: Check that the singular values are sorted in
1222 * non-increasing order and are non-negative
1223 * 24: Compare SBDSVDX with and without singular vectors
1224 *
1225  CALL sbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt,
1226  $ ldpt, work( iwbs+mnmin ), result( 20 ) )
1227  CALL sort01( 'Columns', mnmin, mnmin, u, ldpt,
1228  $ work( iwbs+mnmin ), lwork-mnmin,
1229  $ result( 21 ) )
1230  CALL sort01( 'Rows', mnmin, mnmin, vt, ldpt,
1231  $ work( iwbs+mnmin ), lwork-mnmin,
1232  $ result( 22) )
1233 *
1234  result( 23 ) = zero
1235  DO 180 i = 1, mnmin - 1
1236  IF( s1( i ).LT.s1( i+1 ) )
1237  $ result( 23 ) = ulpinv
1238  IF( s1( i ).LT.zero )
1239  $ result( 23 ) = ulpinv
1240  180 CONTINUE
1241  IF( mnmin.GE.1 ) THEN
1242  IF( s1( mnmin ).LT.zero )
1243  $ result( 23 ) = ulpinv
1244  END IF
1245 *
1246  temp2 = zero
1247  DO 190 j = 1, mnmin
1248  temp1 = abs( s1( j )-s2( j ) ) /
1249  $ max( sqrt( unfl )*max( s1( 1 ), one ),
1250  $ ulp*max( abs( s1( 1 ) ), abs( s2( 1 ) ) ) )
1251  temp2 = max( temp1, temp2 )
1252  190 CONTINUE
1253  result( 24 ) = temp2
1254  anorm = s1( 1 )
1255 *
1256 * Use SBDSVDX with RANGE='I': choose random values for IL and
1257 * IU, and ask for the IL-th through IU-th singular values
1258 * and corresponding vectors.
1259 *
1260  DO 200 i = 1, 4
1261  iseed2( i ) = iseed( i )
1262  200 CONTINUE
1263  IF( mnmin.LE.1 ) THEN
1264  il = 1
1265  iu = mnmin
1266  ELSE
1267  il = 1 + int( ( mnmin-1 )*slarnd( 1, iseed2 ) )
1268  iu = 1 + int( ( mnmin-1 )*slarnd( 1, iseed2 ) )
1269  IF( iu.LT.il ) THEN
1270  itemp = iu
1271  iu = il
1272  il = itemp
1273  END IF
1274  END IF
1275 *
1276  CALL scopy( mnmin, bd, 1, work( iwbd ), 1 )
1277  IF( mnmin.GT.0 )
1278  $ CALL scopy( mnmin-1, be, 1, work( iwbe ), 1 )
1279 *
1280  CALL sbdsvdx( uplo, 'V', 'I', mnmin, work( iwbd ),
1281  $ work( iwbe ), zero, zero, il, iu, ns1, s1,
1282  $ work( iwbz ), mnmin2, work( iwwork ),
1283  $ iwork, iinfo)
1284 *
1285 * Check error code from SBDSVDX.
1286 *
1287  IF( iinfo.NE.0 ) THEN
1288  WRITE( nout, fmt = 9998 )'SBDSVDX(vects,I)', iinfo,
1289  $ m, n, jtype, ioldsd
1290  info = abs( iinfo )
1291  IF( iinfo.LT.0 ) THEN
1292  RETURN
1293  ELSE
1294  result( 25 ) = ulpinv
1295  GO TO 270
1296  END IF
1297  END IF
1298 *
1299  j = iwbz
1300  DO 210 i = 1, ns1
1301  CALL scopy( mnmin, work( j ), 1, u( 1,i ), 1 )
1302  j = j + mnmin
1303  CALL scopy( mnmin, work( j ), 1, vt( i,1 ), ldpt )
1304  j = j + mnmin
1305  210 CONTINUE
1306 *
1307 * Use SBDSVDX to compute only the singular values of the
1308 * bidiagonal matrix B; U and VT should not be modified.
1309 *
1310  CALL scopy( mnmin, bd, 1, work( iwbd ), 1 )
1311  IF( mnmin.GT.0 )
1312  $ CALL scopy( mnmin-1, be, 1, work( iwbe ), 1 )
1313 *
1314  CALL sbdsvdx( uplo, 'N', 'I', mnmin, work( iwbd ),
1315  $ work( iwbe ), zero, zero, il, iu, ns2, s2,
1316  $ work( iwbz ), mnmin2, work( iwwork ),
1317  $ iwork, iinfo )
1318 *
1319 * Check error code from SBDSVDX.
1320 *
1321  IF( iinfo.NE.0 ) THEN
1322  WRITE( nout, fmt = 9998 )'SBDSVDX(values,I)', iinfo,
1323  $ m, n, jtype, ioldsd
1324  info = abs( iinfo )
1325  IF( iinfo.LT.0 ) THEN
1326  RETURN
1327  ELSE
1328  result( 29 ) = ulpinv
1329  GO TO 270
1330  END IF
1331  END IF
1332 *
1333 * Test 25: Check S1 - U' * B * VT'
1334 * 26: Check the orthogonality of U
1335 * 27: Check the orthogonality of VT
1336 * 28: Check that the singular values are sorted in
1337 * non-increasing order and are non-negative
1338 * 29: Compare SBDSVDX with and without singular vectors
1339 *
1340  CALL sbdt04( uplo, mnmin, bd, be, s1, ns1, u,
1341  $ ldpt, vt, ldpt, work( iwbs+mnmin ),
1342  $ result( 25 ) )
1343  CALL sort01( 'Columns', mnmin, ns1, u, ldpt,
1344  $ work( iwbs+mnmin ), lwork-mnmin,
1345  $ result( 26 ) )
1346  CALL sort01( 'Rows', ns1, mnmin, vt, ldpt,
1347  $ work( iwbs+mnmin ), lwork-mnmin,
1348  $ result( 27 ) )
1349 *
1350  result( 28 ) = zero
1351  DO 220 i = 1, ns1 - 1
1352  IF( s1( i ).LT.s1( i+1 ) )
1353  $ result( 28 ) = ulpinv
1354  IF( s1( i ).LT.zero )
1355  $ result( 28 ) = ulpinv
1356  220 CONTINUE
1357  IF( ns1.GE.1 ) THEN
1358  IF( s1( ns1 ).LT.zero )
1359  $ result( 28 ) = ulpinv
1360  END IF
1361 *
1362  temp2 = zero
1363  DO 230 j = 1, ns1
1364  temp1 = abs( s1( j )-s2( j ) ) /
1365  $ max( sqrt( unfl )*max( s1( 1 ), one ),
1366  $ ulp*max( abs( s1( 1 ) ), abs( s2( 1 ) ) ) )
1367  temp2 = max( temp1, temp2 )
1368  230 CONTINUE
1369  result( 29 ) = temp2
1370 *
1371 * Use SBDSVDX with RANGE='V': determine the values VL and VU
1372 * of the IL-th and IU-th singular values and ask for all
1373 * singular values in this range.
1374 *
1375  CALL scopy( mnmin, work( iwbs ), 1, s1, 1 )
1376 *
1377  IF( mnmin.GT.0 ) THEN
1378  IF( il.NE.1 ) THEN
1379  vu = s1( il ) + max( half*abs( s1( il )-s1( il-1 ) ),
1380  $ ulp*anorm, two*rtunfl )
1381  ELSE
1382  vu = s1( 1 ) + max( half*abs( s1( mnmin )-s1( 1 ) ),
1383  $ ulp*anorm, two*rtunfl )
1384  END IF
1385  IF( iu.NE.ns1 ) THEN
1386  vl = s1( iu ) - max( ulp*anorm, two*rtunfl,
1387  $ half*abs( s1( iu+1 )-s1( iu ) ) )
1388  ELSE
1389  vl = s1( ns1 ) - max( ulp*anorm, two*rtunfl,
1390  $ half*abs( s1( mnmin )-s1( 1 ) ) )
1391  END IF
1392  vl = max( vl,zero )
1393  vu = max( vu,zero )
1394  IF( vl.GE.vu ) vu = max( vu*2, vu+vl+half )
1395  ELSE
1396  vl = zero
1397  vu = one
1398  END IF
1399 *
1400  CALL scopy( mnmin, bd, 1, work( iwbd ), 1 )
1401  IF( mnmin.GT.0 )
1402  $ CALL scopy( mnmin-1, be, 1, work( iwbe ), 1 )
1403 *
1404  CALL sbdsvdx( uplo, 'V', 'V', mnmin, work( iwbd ),
1405  $ work( iwbe ), vl, vu, 0, 0, ns1, s1,
1406  $ work( iwbz ), mnmin2, work( iwwork ),
1407  $ iwork, iinfo )
1408 *
1409 * Check error code from SBDSVDX.
1410 *
1411  IF( iinfo.NE.0 ) THEN
1412  WRITE( nout, fmt = 9998 )'SBDSVDX(vects,V)', iinfo,
1413  $ m, n, jtype, ioldsd
1414  info = abs( iinfo )
1415  IF( iinfo.LT.0 ) THEN
1416  RETURN
1417  ELSE
1418  result( 30 ) = ulpinv
1419  GO TO 270
1420  END IF
1421  END IF
1422 *
1423  j = iwbz
1424  DO 240 i = 1, ns1
1425  CALL scopy( mnmin, work( j ), 1, u( 1,i ), 1 )
1426  j = j + mnmin
1427  CALL scopy( mnmin, work( j ), 1, vt( i,1 ), ldpt )
1428  j = j + mnmin
1429  240 CONTINUE
1430 *
1431 * Use SBDSVDX to compute only the singular values of the
1432 * bidiagonal matrix B; U and VT should not be modified.
1433 *
1434  CALL scopy( mnmin, bd, 1, work( iwbd ), 1 )
1435  IF( mnmin.GT.0 )
1436  $ CALL scopy( mnmin-1, be, 1, work( iwbe ), 1 )
1437 *
1438  CALL sbdsvdx( uplo, 'N', 'V', mnmin, work( iwbd ),
1439  $ work( iwbe ), vl, vu, 0, 0, ns2, s2,
1440  $ work( iwbz ), mnmin2, work( iwwork ),
1441  $ iwork, iinfo )
1442 *
1443 * Check error code from SBDSVDX.
1444 *
1445  IF( iinfo.NE.0 ) THEN
1446  WRITE( nout, fmt = 9998 )'SBDSVDX(values,V)', iinfo,
1447  $ m, n, jtype, ioldsd
1448  info = abs( iinfo )
1449  IF( iinfo.LT.0 ) THEN
1450  RETURN
1451  ELSE
1452  result( 34 ) = ulpinv
1453  GO TO 270
1454  END IF
1455  END IF
1456 *
1457 * Test 30: Check S1 - U' * B * VT'
1458 * 31: Check the orthogonality of U
1459 * 32: Check the orthogonality of VT
1460 * 33: Check that the singular values are sorted in
1461 * non-increasing order and are non-negative
1462 * 34: Compare SBDSVDX with and without singular vectors
1463 *
1464  CALL sbdt04( uplo, mnmin, bd, be, s1, ns1, u,
1465  $ ldpt, vt, ldpt, work( iwbs+mnmin ),
1466  $ result( 30 ) )
1467  CALL sort01( 'Columns', mnmin, ns1, u, ldpt,
1468  $ work( iwbs+mnmin ), lwork-mnmin,
1469  $ result( 31 ) )
1470  CALL sort01( 'Rows', ns1, mnmin, vt, ldpt,
1471  $ work( iwbs+mnmin ), lwork-mnmin,
1472  $ result( 32 ) )
1473 *
1474  result( 33 ) = zero
1475  DO 250 i = 1, ns1 - 1
1476  IF( s1( i ).LT.s1( i+1 ) )
1477  $ result( 28 ) = ulpinv
1478  IF( s1( i ).LT.zero )
1479  $ result( 28 ) = ulpinv
1480  250 CONTINUE
1481  IF( ns1.GE.1 ) THEN
1482  IF( s1( ns1 ).LT.zero )
1483  $ result( 28 ) = ulpinv
1484  END IF
1485 *
1486  temp2 = zero
1487  DO 260 j = 1, ns1
1488  temp1 = abs( s1( j )-s2( j ) ) /
1489  $ max( sqrt( unfl )*max( s1( 1 ), one ),
1490  $ ulp*max( abs( s1( 1 ) ), abs( s2( 1 ) ) ) )
1491  temp2 = max( temp1, temp2 )
1492  260 CONTINUE
1493  result( 34 ) = temp2
1494 *
1495 * End of Loop -- Check for RESULT(j) > THRESH
1496 *
1497  270 CONTINUE
1498 *
1499  DO 280 j = 1, 34
1500  IF( result( j ).GE.thresh ) THEN
1501  IF( nfail.EQ.0 )
1502  $ CALL slahd2( nout, path )
1503  WRITE( nout, fmt = 9999 )m, n, jtype, ioldsd, j,
1504  $ result( j )
1505  nfail = nfail + 1
1506  END IF
1507  280 CONTINUE
1508  IF( .NOT.bidiag ) THEN
1509  ntest = ntest + 34
1510  ELSE
1511  ntest = ntest + 30
1512  END IF
1513 *
1514  290 CONTINUE
1515  300 CONTINUE
1516 *
1517 * Summary
1518 *
1519  CALL alasum( path, nout, nfail, ntest, 0 )
1520 *
1521  RETURN
1522 *
1523 * End of SCHKBD
1524 *
1525  9999 FORMAT( ' M=', i5, ', N=', i5, ', type ', i2, ', seed=',
1526  $ 4( i4, ',' ), ' test(', i2, ')=', g11.4 )
1527  9998 FORMAT( ' SCHKBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1528  $ i6, ', N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
1529  $ i5, ')' )
1530 *
1531  END
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine slatmr(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)
SLATMR
Definition: slatmr.f:473
subroutine sbdt04(UPLO, N, D, E, S, NS, U, LDU, VT, LDVT, WORK, RESID)
Definition: sbdt04.f:132
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
subroutine slahd2(IOUNIT, PATH)
SLAHD2
Definition: slahd2.f:67
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine sort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
SORT01
Definition: sort01.f:118
subroutine slatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
SLATMS
Definition: slatms.f:323
subroutine sbdt01(M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK, RESID)
SBDT01
Definition: sbdt01.f:142
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:112
subroutine sbdsvdx(UPLO, JOBZ, RANGE, N, D, E, VL, VU, IL, IU, NS, S, Z, LDZ, WORK, IWORK, INFO)
SBDSVDX
Definition: sbdsvdx.f:228
subroutine sbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SBDSQR
Definition: sbdsqr.f:232
subroutine sbdsdc(UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
SBDSDC
Definition: sbdsdc.f:207
subroutine schkbd(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)
SCHKBD
Definition: schkbd.f:495
subroutine sgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
SGEBRD
Definition: sgebrd.f:207
subroutine sbdt03(UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK, RESID)
SBDT03
Definition: sbdt03.f:137
subroutine sorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGBR
Definition: sorgbr.f:159
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
subroutine sbdt02(M, N, B, LDB, C, LDC, U, LDU, WORK, RESID)
SBDT02
Definition: sbdt02.f:113
subroutine alasum(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASUM
Definition: alasum.f:75