LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
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 *> 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 SGEBRD and SORGBR
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 SBDSQR 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 *> SSVDCH)
106 *>
107 *> Test SBDSQR 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 SBDSDC 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) SGEBRD 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, SCHKBD
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 SBDSQR. 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 SCHKBD to continue the same random
240 *> number sequence.
241 *> \endverbatim
242 *>
243 *> \param[in] THRESH
244 *> \verbatim
245 *> THRESH is REAL
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 REAL 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 REAL array, dimension
269 *> (max(min(MVAL(j),NVAL(j))))
270 *> \endverbatim
271 *>
272 *> \param[out] BE
273 *> \verbatim
274 *> BE is REAL array, dimension
275 *> (max(min(MVAL(j),NVAL(j))))
276 *> \endverbatim
277 *>
278 *> \param[out] S1
279 *> \verbatim
280 *> S1 is REAL array, dimension
281 *> (max(min(MVAL(j),NVAL(j))))
282 *> \endverbatim
283 *>
284 *> \param[out] S2
285 *> \verbatim
286 *> S2 is REAL array, dimension
287 *> (max(min(MVAL(j),NVAL(j))))
288 *> \endverbatim
289 *>
290 *> \param[out] X
291 *> \verbatim
292 *> X is REAL 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 REAL array, dimension (LDX,NRHS)
305 *> \endverbatim
306 *>
307 *> \param[out] Z
308 *> \verbatim
309 *> Z is REAL array, dimension (LDX,NRHS)
310 *> \endverbatim
311 *>
312 *> \param[out] Q
313 *> \verbatim
314 *> Q is REAL 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 REAL 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 REAL array, dimension
338 *> (LDPT,max(min(MVAL(j),NVAL(j))))
339 *> \endverbatim
340 *>
341 *> \param[out] VT
342 *> \verbatim
343 *> VT is REAL array, dimension
344 *> (LDPT,max(min(MVAL(j),NVAL(j))))
345 *> \endverbatim
346 *>
347 *> \param[out] WORK
348 *> \verbatim
349 *> WORK is REAL 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 SLATMR, SLATMS, SGEBRD, SORGBR, or SBDSQR,
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 single_eig
432 *
433 * =====================================================================
434  SUBROUTINE schkbd( 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  REAL thresh
448 * ..
449 * .. Array Arguments ..
450  LOGICAL dotype( * )
451  INTEGER iseed( 4 ), iwork( * ), mval( * ), nval( * )
452  REAL 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  REAL zero, one, two, half
462  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
463  $ half = 0.5e0 )
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  REAL 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  REAL dum( 1 ), dumma( 1 ), result( 19 )
481 * ..
482 * .. External Functions ..
483  REAL slamch, slarnd
484  EXTERNAL slamch, slarnd
485 * ..
486 * .. External Subroutines ..
487  EXTERNAL alasum, sbdsdc, sbdsqr, sbdt01, sbdt02, sbdt03,
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( 'SCHKBD', -info )
560  RETURN
561  END IF
562 *
563 * Initialize constants
564 *
565  path( 1: 1 ) = 'Single precision'
566  path( 2: 3 ) = 'BD'
567  nfail = 0
568  ntest = 0
569  unfl = slamch( 'Safe minimum' )
570  ovfl = slamch( 'Overflow' )
571  CALL slabad( unfl, ovfl )
572  ulp = slamch( '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 slaset( '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 slatms( 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 slatms( 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 slatms( 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 slatmr( 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 slatmr( 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 slatmr( 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*slarnd( 2, iseed ) )
724  IF( j.LT.mnmin )
725  $ be( j ) = exp( temp1*slarnd( 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 slatmr( 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 slatmr( 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 SGEBRD and SORGBR 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 slacpy( ' ', m, n, a, lda, q, ldq )
777  CALL sgebrd( 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 SGEBRD.
781 *
782  IF( iinfo.NE.0 ) THEN
783  WRITE( nout, fmt = 9998 )'SGEBRD', iinfo, m, n,
784  $ jtype, ioldsd
785  info = abs( iinfo )
786  RETURN
787  END IF
788 *
789  CALL slacpy( ' ', 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 sorgbr( 'Q', m, mq, n, q, ldq, work,
802  $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
803 *
804 * Check error code from SORGBR.
805 *
806  IF( iinfo.NE.0 ) THEN
807  WRITE( nout, fmt = 9998 )'SORGBR(Q)', iinfo, m, n,
808  $ jtype, ioldsd
809  info = abs( iinfo )
810  RETURN
811  END IF
812 *
813 * Generate P'
814 *
815  CALL sorgbr( 'P', mnmin, n, m, pt, ldpt, work( mnmin+1 ),
816  $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
817 *
818 * Check error code from SORGBR.
819 *
820  IF( iinfo.NE.0 ) THEN
821  WRITE( nout, fmt = 9998 )'SORGBR(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 sgemm( '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 sbdt01( m, n, 1, a, lda, q, ldq, bd, be, pt, ldpt,
837  $ work, result( 1 ) )
838  CALL sort01( 'Columns', m, mq, q, ldq, work, lwork,
839  $ result( 2 ) )
840  CALL sort01( 'Rows', mnmin, n, pt, ldpt, work, lwork,
841  $ result( 3 ) )
842  END IF
843 *
844 * Use SBDSQR to form the SVD of the bidiagonal matrix B:
845 * B := U * S1 * VT, and compute Z = U' * Y.
846 *
847  CALL scopy( mnmin, bd, 1, s1, 1 )
848  IF( mnmin.GT.0 )
849  $ CALL scopy( mnmin-1, be, 1, work, 1 )
850  CALL slacpy( ' ', m, nrhs, y, ldx, z, ldx )
851  CALL slaset( 'Full', mnmin, mnmin, zero, one, u, ldpt )
852  CALL slaset( 'Full', mnmin, mnmin, zero, one, vt, ldpt )
853 *
854  CALL sbdsqr( uplo, mnmin, mnmin, mnmin, nrhs, s1, work, vt,
855  $ ldpt, u, ldpt, z, ldx, work( mnmin+1 ), iinfo )
856 *
857 * Check error code from SBDSQR.
858 *
859  IF( iinfo.NE.0 ) THEN
860  WRITE( nout, fmt = 9998 )'SBDSQR(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 SBDSQR to compute only the singular values of the
872 * bidiagonal matrix B; U, VT, and Z should not be modified.
873 *
874  CALL scopy( mnmin, bd, 1, s2, 1 )
875  IF( mnmin.GT.0 )
876  $ CALL scopy( mnmin-1, be, 1, work, 1 )
877 *
878  CALL sbdsqr( uplo, mnmin, 0, 0, 0, s2, work, vt, ldpt, u,
879  $ ldpt, z, ldx, work( mnmin+1 ), iinfo )
880 *
881 * Check error code from SBDSQR.
882 *
883  IF( iinfo.NE.0 ) THEN
884  WRITE( nout, fmt = 9998 )'SBDSQR(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 sbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
901  $ work, result( 4 ) )
902  CALL sbdt02( mnmin, nrhs, y, ldx, z, ldx, u, ldpt, work,
903  $ result( 5 ) )
904  CALL sort01( 'Columns', mnmin, mnmin, u, ldpt, work, lwork,
905  $ result( 6 ) )
906  CALL sort01( '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 SBDSQR 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 SSVDCH( 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 SBDSQR 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 scopy( mnmin, bd, 1, s2, 1 )
957  IF( mnmin.GT.0 )
958  $ CALL scopy( mnmin-1, be, 1, work, 1 )
959 *
960  CALL sbdsqr( 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 sbdt01( m, n, 0, a, lda, q, ldq, s2, dumma, pt,
969  $ ldpt, work, result( 11 ) )
970  CALL sbdt02( m, nrhs, x, ldx, y, ldx, q, ldq, work,
971  $ result( 12 ) )
972  CALL sort01( 'Columns', m, mq, q, ldq, work, lwork,
973  $ result( 13 ) )
974  CALL sort01( 'Rows', mnmin, n, pt, ldpt, work, lwork,
975  $ result( 14 ) )
976  END IF
977 *
978 * Use SBDSDC to form the SVD of the bidiagonal matrix B:
979 * B := U * S1 * VT
980 *
981  CALL scopy( mnmin, bd, 1, s1, 1 )
982  IF( mnmin.GT.0 )
983  $ CALL scopy( mnmin-1, be, 1, work, 1 )
984  CALL slaset( 'Full', mnmin, mnmin, zero, one, u, ldpt )
985  CALL slaset( 'Full', mnmin, mnmin, zero, one, vt, ldpt )
986 *
987  CALL sbdsdc( uplo, 'I', mnmin, s1, work, u, ldpt, vt, ldpt,
988  $ dum, idum, work( mnmin+1 ), iwork, iinfo )
989 *
990 * Check error code from SBDSDC.
991 *
992  IF( iinfo.NE.0 ) THEN
993  WRITE( nout, fmt = 9998 )'SBDSDC(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 SBDSDC to compute only the singular values of the
1005 * bidiagonal matrix B; U and VT should not be modified.
1006 *
1007  CALL scopy( mnmin, bd, 1, s2, 1 )
1008  IF( mnmin.GT.0 )
1009  $ CALL scopy( mnmin-1, be, 1, work, 1 )
1010 *
1011  CALL sbdsdc( uplo, 'N', mnmin, s2, work, dum, 1, dum, 1,
1012  $ dum, idum, work( mnmin+1 ), iwork, iinfo )
1013 *
1014 * Check error code from SBDSDC.
1015 *
1016  IF( iinfo.NE.0 ) THEN
1017  WRITE( nout, fmt = 9998 )'SBDSDC(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 sbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
1033  $ work, result( 15 ) )
1034  CALL sort01( 'Columns', mnmin, mnmin, u, ldpt, work, lwork,
1035  $ result( 16 ) )
1036  CALL sort01( '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 SBDSQR 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 slahd2( 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 SCHKBD
1095 *
1096  9999 FORMAT( ' M=', i5, ', N=', i5, ', type ', i2, ', seed=',
1097  $ 4( i4, ',' ), ' test(', i2, ')=', g11.4 )
1098  9998 FORMAT( ' SCHKBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1099  $ i6, ', N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
1100  $ i5, ')' )
1101 *
1102  END
subroutine sgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
SGEBRD
Definition: sgebrd.f:205
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:111
subroutine sbdt02(M, N, B, LDB, C, LDC, U, LDU, WORK, RESID)
SBDT02
Definition: sbdt02.f:112
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:434
subroutine sbdt03(UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK, RESID)
SBDT03
Definition: sbdt03.f:135
subroutine slahd2(IOUNIT, PATH)
SLAHD2
Definition: slahd2.f:66
subroutine slatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
SLATMS
Definition: slatms.f:321
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:469
subroutine sorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGBR
Definition: sorgbr.f:158
subroutine sbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SBDSQR
Definition: sbdsqr.f:230
subroutine sbdsdc(UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
SBDSDC
Definition: sbdsdc.f:205
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:61
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:52
subroutine alasum(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASUM
Definition: alasum.f:74
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:104
subroutine sort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
SORT01
Definition: sort01.f:117
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:188
real function slarnd(IDIST, ISEED)
SLARND
Definition: slarnd.f:74
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
Definition: xerbla-fortran:9
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:75
subroutine sbdt01(M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK, RESID)
SBDT01
Definition: sbdt01.f:140