LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> DBDSVDX 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 DGEBRD and DORGBR
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 DBDSQR 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*> DSVDCH)
125*>
126*> Test DBDSQR 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 DBDSDC 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 DBDSVDX 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 ) DBDSVDX('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 ) DBDSVDX('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) DGEBRD 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, DCHKBD
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 DBDSQR. 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 DCHKBD to continue the same random
297*> number sequence.
298*> \endverbatim
299*>
300*> \param[in] THRESH
301*> \verbatim
302*> THRESH is DOUBLE PRECISION
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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension
326*> (max(min(MVAL(j),NVAL(j))))
327*> \endverbatim
328*>
329*> \param[out] BE
330*> \verbatim
331*> BE is DOUBLE PRECISION array, dimension
332*> (max(min(MVAL(j),NVAL(j))))
333*> \endverbatim
334*>
335*> \param[out] S1
336*> \verbatim
337*> S1 is DOUBLE PRECISION array, dimension
338*> (max(min(MVAL(j),NVAL(j))))
339*> \endverbatim
340*>
341*> \param[out] S2
342*> \verbatim
343*> S2 is DOUBLE PRECISION array, dimension
344*> (max(min(MVAL(j),NVAL(j))))
345*> \endverbatim
346*>
347*> \param[out] X
348*> \verbatim
349*> X is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDX,NRHS)
362*> \endverbatim
363*>
364*> \param[out] Z
365*> \verbatim
366*> Z is DOUBLE PRECISION array, dimension (LDX,NRHS)
367*> \endverbatim
368*>
369*> \param[out] Q
370*> \verbatim
371*> Q is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension
395*> (LDPT,max(min(MVAL(j),NVAL(j))))
396*> \endverbatim
397*>
398*> \param[out] VT
399*> \verbatim
400*> VT is DOUBLE PRECISION array, dimension
401*> (LDPT,max(min(MVAL(j),NVAL(j))))
402*> \endverbatim
403*>
404*> \param[out] WORK
405*> \verbatim
406*> WORK is DOUBLE PRECISION 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 DLATMR, SLATMS, DGEBRD, DORGBR, or DBDSQR,
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*> \ingroup double_eig
487*
488* =====================================================================
489 SUBROUTINE dchkbd( NSIZES, MVAL, NVAL, NTYPES, DOTYPE, NRHS,
490 $ ISEED, THRESH, A, LDA, BD, BE, S1, S2, X, LDX,
491 $ Y, Z, Q, LDQ, PT, LDPT, U, VT, WORK, LWORK,
492 $ IWORK, NOUT, INFO )
493*
494* -- LAPACK test routine --
495* -- LAPACK is a software package provided by Univ. of Tennessee, --
496* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
497*
498* .. Scalar Arguments ..
499 INTEGER INFO, LDA, LDPT, LDQ, LDX, LWORK, NOUT, NRHS,
500 $ NSIZES, NTYPES
501 DOUBLE PRECISION THRESH
502* ..
503* .. Array Arguments ..
504 LOGICAL DOTYPE( * )
505 INTEGER ISEED( 4 ), IWORK( * ), MVAL( * ), NVAL( * )
506 DOUBLE PRECISION A( LDA, * ), BD( * ), BE( * ), PT( LDPT, * ),
507 $ q( ldq, * ), s1( * ), s2( * ), u( ldpt, * ),
508 $ vt( ldpt, * ), work( * ), x( ldx, * ),
509 $ y( ldx, * ), z( ldx, * )
510* ..
511*
512* ======================================================================
513*
514* .. Parameters ..
515 DOUBLE PRECISION ZERO, ONE, TWO, HALF
516 PARAMETER ( ZERO = 0.0d0, one = 1.0d0, two = 2.0d0,
517 $ half = 0.5d0 )
518 INTEGER MAXTYP
519 parameter( maxtyp = 16 )
520* ..
521* .. Local Scalars ..
522 LOGICAL BADMM, BADNN, BIDIAG
523 CHARACTER UPLO
524 CHARACTER*3 PATH
525 INTEGER I, IINFO, IL, IMODE, ITEMP, ITYPE, IU, IWBD,
526 $ iwbe, iwbs, iwbz, iwwork, j, jcol, jsize,
527 $ jtype, log2ui, m, minwrk, mmax, mnmax, mnmin,
528 $ mnmin2, mq, mtypes, n, nfail, nmax,
529 $ ns1, ns2, ntest
530 DOUBLE PRECISION ABSTOL, AMNINV, ANORM, COND, OVFL, RTOVFL,
531 $ RTUNFL, TEMP1, TEMP2, ULP, ULPINV, UNFL,
532 $ VL, VU
533* ..
534* .. Local Arrays ..
535 INTEGER IDUM( 1 ), IOLDSD( 4 ), ISEED2( 4 ),
536 $ KMAGN( MAXTYP ), KMODE( MAXTYP ),
537 $ KTYPE( MAXTYP )
538 DOUBLE PRECISION DUM( 1 ), DUMMA( 1 ), RESULT( 40 )
539* ..
540* .. External Functions ..
541 DOUBLE PRECISION DLAMCH, DLARND, DSXT1
542 EXTERNAL DLAMCH, DLARND, DSXT1
543* ..
544* .. External Subroutines ..
545 EXTERNAL alasum, dbdsdc, dbdsqr, dbdsvdx, dbdt01,
549* ..
550* .. Intrinsic Functions ..
551 INTRINSIC abs, exp, int, log, max, min, sqrt
552* ..
553* .. Scalars in Common ..
554 LOGICAL LERR, OK
555 CHARACTER*32 SRNAMT
556 INTEGER INFOT, NUNIT
557* ..
558* .. Common blocks ..
559 COMMON / infoc / infot, nunit, ok, lerr
560 COMMON / srnamc / srnamt
561* ..
562* .. Data statements ..
563 DATA ktype / 1, 2, 5*4, 5*6, 3*9, 10 /
564 DATA kmagn / 2*1, 3*1, 2, 3, 3*1, 2, 3, 1, 2, 3, 0 /
565 DATA kmode / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
566 $ 0, 0, 0 /
567* ..
568* .. Executable Statements ..
569*
570* Check for errors
571*
572 info = 0
573*
574 badmm = .false.
575 badnn = .false.
576 mmax = 1
577 nmax = 1
578 mnmax = 1
579 minwrk = 1
580 DO 10 j = 1, nsizes
581 mmax = max( mmax, mval( j ) )
582 IF( mval( j ).LT.0 )
583 $ badmm = .true.
584 nmax = max( nmax, nval( j ) )
585 IF( nval( j ).LT.0 )
586 $ badnn = .true.
587 mnmax = max( mnmax, min( mval( j ), nval( j ) ) )
588 minwrk = max( minwrk, 3*( mval( j )+nval( j ) ),
589 $ mval( j )*( mval( j )+max( mval( j ), nval( j ),
590 $ nrhs )+1 )+nval( j )*min( nval( j ), mval( j ) ) )
591 10 CONTINUE
592*
593* Check for errors
594*
595 IF( nsizes.LT.0 ) THEN
596 info = -1
597 ELSE IF( badmm ) THEN
598 info = -2
599 ELSE IF( badnn ) THEN
600 info = -3
601 ELSE IF( ntypes.LT.0 ) THEN
602 info = -4
603 ELSE IF( nrhs.LT.0 ) THEN
604 info = -6
605 ELSE IF( lda.LT.mmax ) THEN
606 info = -11
607 ELSE IF( ldx.LT.mmax ) THEN
608 info = -17
609 ELSE IF( ldq.LT.mmax ) THEN
610 info = -21
611 ELSE IF( ldpt.LT.mnmax ) THEN
612 info = -23
613 ELSE IF( minwrk.GT.lwork ) THEN
614 info = -27
615 END IF
616*
617 IF( info.NE.0 ) THEN
618 CALL xerbla( 'DCHKBD', -info )
619 RETURN
620 END IF
621*
622* Initialize constants
623*
624 path( 1: 1 ) = 'Double precision'
625 path( 2: 3 ) = 'BD'
626 nfail = 0
627 ntest = 0
628 unfl = dlamch( 'Safe minimum' )
629 ovfl = dlamch( 'Overflow' )
630 ulp = dlamch( 'Precision' )
631 ulpinv = one / ulp
632 log2ui = int( log( ulpinv ) / log( two ) )
633 rtunfl = sqrt( unfl )
634 rtovfl = sqrt( ovfl )
635 infot = 0
636 abstol = 2*unfl
637*
638* Loop over sizes, types
639*
640 DO 300 jsize = 1, nsizes
641 m = mval( jsize )
642 n = nval( jsize )
643 mnmin = min( m, n )
644 amninv = one / max( m, n, 1 )
645*
646 IF( nsizes.NE.1 ) THEN
647 mtypes = min( maxtyp, ntypes )
648 ELSE
649 mtypes = min( maxtyp+1, ntypes )
650 END IF
651*
652 DO 290 jtype = 1, mtypes
653 IF( .NOT.dotype( jtype ) )
654 $ GO TO 290
655*
656 DO 20 j = 1, 4
657 ioldsd( j ) = iseed( j )
658 20 CONTINUE
659*
660 DO 30 j = 1, 34
661 result( j ) = -one
662 30 CONTINUE
663*
664 uplo = ' '
665*
666* Compute "A"
667*
668* Control parameters:
669*
670* KMAGN KMODE KTYPE
671* =1 O(1) clustered 1 zero
672* =2 large clustered 2 identity
673* =3 small exponential (none)
674* =4 arithmetic diagonal, (w/ eigenvalues)
675* =5 random symmetric, w/ eigenvalues
676* =6 nonsymmetric, w/ singular values
677* =7 random diagonal
678* =8 random symmetric
679* =9 random nonsymmetric
680* =10 random bidiagonal (log. distrib.)
681*
682 IF( mtypes.GT.maxtyp )
683 $ GO TO 100
684*
685 itype = ktype( jtype )
686 imode = kmode( jtype )
687*
688* Compute norm
689*
690 GO TO ( 40, 50, 60 )kmagn( jtype )
691*
692 40 CONTINUE
693 anorm = one
694 GO TO 70
695*
696 50 CONTINUE
697 anorm = ( rtovfl*ulp )*amninv
698 GO TO 70
699*
700 60 CONTINUE
701 anorm = rtunfl*max( m, n )*ulpinv
702 GO TO 70
703*
704 70 CONTINUE
705*
706 CALL dlaset( 'Full', lda, n, zero, zero, a, lda )
707 iinfo = 0
708 cond = ulpinv
709*
710 bidiag = .false.
711 IF( itype.EQ.1 ) THEN
712*
713* Zero matrix
714*
715 iinfo = 0
716*
717 ELSE IF( itype.EQ.2 ) THEN
718*
719* Identity
720*
721 DO 80 jcol = 1, mnmin
722 a( jcol, jcol ) = anorm
723 80 CONTINUE
724*
725 ELSE IF( itype.EQ.4 ) THEN
726*
727* Diagonal Matrix, [Eigen]values Specified
728*
729 CALL dlatms( mnmin, mnmin, 'S', iseed, 'N', work, imode,
730 $ cond, anorm, 0, 0, 'N', a, lda,
731 $ work( mnmin+1 ), iinfo )
732*
733 ELSE IF( itype.EQ.5 ) THEN
734*
735* Symmetric, eigenvalues specified
736*
737 CALL dlatms( mnmin, mnmin, 'S', iseed, 'S', work, imode,
738 $ cond, anorm, m, n, 'N', a, lda,
739 $ work( mnmin+1 ), iinfo )
740*
741 ELSE IF( itype.EQ.6 ) THEN
742*
743* Nonsymmetric, singular values specified
744*
745 CALL dlatms( m, n, 'S', iseed, 'N', work, imode, cond,
746 $ anorm, m, n, 'N', a, lda, work( mnmin+1 ),
747 $ iinfo )
748*
749 ELSE IF( itype.EQ.7 ) THEN
750*
751* Diagonal, random entries
752*
753 CALL dlatmr( mnmin, mnmin, 'S', iseed, 'N', work, 6, one,
754 $ one, 'T', 'N', work( mnmin+1 ), 1, one,
755 $ work( 2*mnmin+1 ), 1, one, 'N', iwork, 0, 0,
756 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
757*
758 ELSE IF( itype.EQ.8 ) THEN
759*
760* Symmetric, random entries
761*
762 CALL dlatmr( mnmin, mnmin, 'S', iseed, 'S', work, 6, one,
763 $ one, 'T', 'N', work( mnmin+1 ), 1, one,
764 $ work( m+mnmin+1 ), 1, one, 'N', iwork, m, n,
765 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
766*
767 ELSE IF( itype.EQ.9 ) THEN
768*
769* Nonsymmetric, random entries
770*
771 CALL dlatmr( m, n, 'S', iseed, 'N', work, 6, one, one,
772 $ 'T', 'N', work( mnmin+1 ), 1, one,
773 $ work( m+mnmin+1 ), 1, one, 'N', iwork, m, n,
774 $ zero, anorm, 'NO', a, lda, iwork, iinfo )
775*
776 ELSE IF( itype.EQ.10 ) THEN
777*
778* Bidiagonal, random entries
779*
780 temp1 = -two*log( ulp )
781 DO 90 j = 1, mnmin
782 bd( j ) = exp( temp1*dlarnd( 2, iseed ) )
783 IF( j.LT.mnmin )
784 $ be( j ) = exp( temp1*dlarnd( 2, iseed ) )
785 90 CONTINUE
786*
787 iinfo = 0
788 bidiag = .true.
789 IF( m.GE.n ) THEN
790 uplo = 'U'
791 ELSE
792 uplo = 'L'
793 END IF
794 ELSE
795 iinfo = 1
796 END IF
797*
798 IF( iinfo.EQ.0 ) THEN
799*
800* Generate Right-Hand Side
801*
802 IF( bidiag ) THEN
803 CALL dlatmr( mnmin, nrhs, 'S', iseed, 'N', work, 6,
804 $ one, one, 'T', 'N', work( mnmin+1 ), 1,
805 $ one, work( 2*mnmin+1 ), 1, one, 'N',
806 $ iwork, mnmin, nrhs, zero, one, 'NO', y,
807 $ ldx, iwork, iinfo )
808 ELSE
809 CALL dlatmr( m, nrhs, 'S', iseed, 'N', work, 6, one,
810 $ one, 'T', 'N', work( m+1 ), 1, one,
811 $ work( 2*m+1 ), 1, one, 'N', iwork, m,
812 $ nrhs, zero, one, 'NO', x, ldx, iwork,
813 $ iinfo )
814 END IF
815 END IF
816*
817* Error Exit
818*
819 IF( iinfo.NE.0 ) THEN
820 WRITE( nout, fmt = 9998 )'Generator', iinfo, m, n,
821 $ jtype, ioldsd
822 info = abs( iinfo )
823 RETURN
824 END IF
825*
826 100 CONTINUE
827*
828* Call DGEBRD and DORGBR to compute B, Q, and P, do tests.
829*
830 IF( .NOT.bidiag ) THEN
831*
832* Compute transformations to reduce A to bidiagonal form:
833* B := Q' * A * P.
834*
835 CALL dlacpy( ' ', m, n, a, lda, q, ldq )
836 CALL dgebrd( m, n, q, ldq, bd, be, work, work( mnmin+1 ),
837 $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
838*
839* Check error code from DGEBRD.
840*
841 IF( iinfo.NE.0 ) THEN
842 WRITE( nout, fmt = 9998 )'DGEBRD', iinfo, m, n,
843 $ jtype, ioldsd
844 info = abs( iinfo )
845 RETURN
846 END IF
847*
848 CALL dlacpy( ' ', m, n, q, ldq, pt, ldpt )
849 IF( m.GE.n ) THEN
850 uplo = 'U'
851 ELSE
852 uplo = 'L'
853 END IF
854*
855* Generate Q
856*
857 mq = m
858 IF( nrhs.LE.0 )
859 $ mq = mnmin
860 CALL dorgbr( 'Q', m, mq, n, q, ldq, work,
861 $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
862*
863* Check error code from DORGBR.
864*
865 IF( iinfo.NE.0 ) THEN
866 WRITE( nout, fmt = 9998 )'DORGBR(Q)', iinfo, m, n,
867 $ jtype, ioldsd
868 info = abs( iinfo )
869 RETURN
870 END IF
871*
872* Generate P'
873*
874 CALL dorgbr( 'P', mnmin, n, m, pt, ldpt, work( mnmin+1 ),
875 $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
876*
877* Check error code from DORGBR.
878*
879 IF( iinfo.NE.0 ) THEN
880 WRITE( nout, fmt = 9998 )'DORGBR(P)', iinfo, m, n,
881 $ jtype, ioldsd
882 info = abs( iinfo )
883 RETURN
884 END IF
885*
886* Apply Q' to an M by NRHS matrix X: Y := Q' * X.
887*
888 CALL dgemm( 'Transpose', 'No transpose', m, nrhs, m, one,
889 $ q, ldq, x, ldx, zero, y, ldx )
890*
891* Test 1: Check the decomposition A := Q * B * PT
892* 2: Check the orthogonality of Q
893* 3: Check the orthogonality of PT
894*
895 CALL dbdt01( m, n, 1, a, lda, q, ldq, bd, be, pt, ldpt,
896 $ work, result( 1 ) )
897 CALL dort01( 'Columns', m, mq, q, ldq, work, lwork,
898 $ result( 2 ) )
899 CALL dort01( 'Rows', mnmin, n, pt, ldpt, work, lwork,
900 $ result( 3 ) )
901 END IF
902*
903* Use DBDSQR to form the SVD of the bidiagonal matrix B:
904* B := U * S1 * VT, and compute Z = U' * Y.
905*
906 CALL dcopy( mnmin, bd, 1, s1, 1 )
907 IF( mnmin.GT.0 )
908 $ CALL dcopy( mnmin-1, be, 1, work, 1 )
909 CALL dlacpy( ' ', m, nrhs, y, ldx, z, ldx )
910 CALL dlaset( 'Full', mnmin, mnmin, zero, one, u, ldpt )
911 CALL dlaset( 'Full', mnmin, mnmin, zero, one, vt, ldpt )
912*
913 CALL dbdsqr( uplo, mnmin, mnmin, mnmin, nrhs, s1, work, vt,
914 $ ldpt, u, ldpt, z, ldx, work( mnmin+1 ), iinfo )
915*
916* Check error code from DBDSQR.
917*
918 IF( iinfo.NE.0 ) THEN
919 WRITE( nout, fmt = 9998 )'DBDSQR(vects)', iinfo, m, n,
920 $ jtype, ioldsd
921 info = abs( iinfo )
922 IF( iinfo.LT.0 ) THEN
923 RETURN
924 ELSE
925 result( 4 ) = ulpinv
926 GO TO 270
927 END IF
928 END IF
929*
930* Use DBDSQR to compute only the singular values of the
931* bidiagonal matrix B; U, VT, and Z should not be modified.
932*
933 CALL dcopy( mnmin, bd, 1, s2, 1 )
934 IF( mnmin.GT.0 )
935 $ CALL dcopy( mnmin-1, be, 1, work, 1 )
936*
937 CALL dbdsqr( uplo, mnmin, 0, 0, 0, s2, work, vt, ldpt, u,
938 $ ldpt, z, ldx, work( mnmin+1 ), iinfo )
939*
940* Check error code from DBDSQR.
941*
942 IF( iinfo.NE.0 ) THEN
943 WRITE( nout, fmt = 9998 )'DBDSQR(values)', iinfo, m, n,
944 $ jtype, ioldsd
945 info = abs( iinfo )
946 IF( iinfo.LT.0 ) THEN
947 RETURN
948 ELSE
949 result( 9 ) = ulpinv
950 GO TO 270
951 END IF
952 END IF
953*
954* Test 4: Check the decomposition B := U * S1 * VT
955* 5: Check the computation Z := U' * Y
956* 6: Check the orthogonality of U
957* 7: Check the orthogonality of VT
958*
959 CALL dbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
960 $ work, result( 4 ) )
961 CALL dbdt02( mnmin, nrhs, y, ldx, z, ldx, u, ldpt, work,
962 $ result( 5 ) )
963 CALL dort01( 'Columns', mnmin, mnmin, u, ldpt, work, lwork,
964 $ result( 6 ) )
965 CALL dort01( 'Rows', mnmin, mnmin, vt, ldpt, work, lwork,
966 $ result( 7 ) )
967*
968* Test 8: Check that the singular values are sorted in
969* non-increasing order and are non-negative
970*
971 result( 8 ) = zero
972 DO 110 i = 1, mnmin - 1
973 IF( s1( i ).LT.s1( i+1 ) )
974 $ result( 8 ) = ulpinv
975 IF( s1( i ).LT.zero )
976 $ result( 8 ) = ulpinv
977 110 CONTINUE
978 IF( mnmin.GE.1 ) THEN
979 IF( s1( mnmin ).LT.zero )
980 $ result( 8 ) = ulpinv
981 END IF
982*
983* Test 9: Compare DBDSQR with and without singular vectors
984*
985 temp2 = zero
986*
987 DO 120 j = 1, mnmin
988 temp1 = abs( s1( j )-s2( j ) ) /
989 $ max( sqrt( unfl )*max( s1( 1 ), one ),
990 $ ulp*max( abs( s1( j ) ), abs( s2( j ) ) ) )
991 temp2 = max( temp1, temp2 )
992 120 CONTINUE
993*
994 result( 9 ) = temp2
995*
996* Test 10: Sturm sequence test of singular values
997* Go up by factors of two until it succeeds
998*
999 temp1 = thresh*( half-ulp )
1000*
1001 DO 130 j = 0, log2ui
1002* CALL DSVDCH( MNMIN, BD, BE, S1, TEMP1, IINFO )
1003 IF( iinfo.EQ.0 )
1004 $ GO TO 140
1005 temp1 = temp1*two
1006 130 CONTINUE
1007*
1008 140 CONTINUE
1009 result( 10 ) = temp1
1010*
1011* Use DBDSQR to form the decomposition A := (QU) S (VT PT)
1012* from the bidiagonal form A := Q B PT.
1013*
1014 IF( .NOT.bidiag ) THEN
1015 CALL dcopy( mnmin, bd, 1, s2, 1 )
1016 IF( mnmin.GT.0 )
1017 $ CALL dcopy( mnmin-1, be, 1, work, 1 )
1018*
1019 CALL dbdsqr( uplo, mnmin, n, m, nrhs, s2, work, pt, ldpt,
1020 $ q, ldq, y, ldx, work( mnmin+1 ), iinfo )
1021*
1022* Test 11: Check the decomposition A := Q*U * S2 * VT*PT
1023* 12: Check the computation Z := U' * Q' * X
1024* 13: Check the orthogonality of Q*U
1025* 14: Check the orthogonality of VT*PT
1026*
1027 CALL dbdt01( m, n, 0, a, lda, q, ldq, s2, dumma, pt,
1028 $ ldpt, work, result( 11 ) )
1029 CALL dbdt02( m, nrhs, x, ldx, y, ldx, q, ldq, work,
1030 $ result( 12 ) )
1031 CALL dort01( 'Columns', m, mq, q, ldq, work, lwork,
1032 $ result( 13 ) )
1033 CALL dort01( 'Rows', mnmin, n, pt, ldpt, work, lwork,
1034 $ result( 14 ) )
1035 END IF
1036*
1037* Use DBDSDC to form the SVD of the bidiagonal matrix B:
1038* B := U * S1 * VT
1039*
1040 CALL dcopy( mnmin, bd, 1, s1, 1 )
1041 IF( mnmin.GT.0 )
1042 $ CALL dcopy( mnmin-1, be, 1, work, 1 )
1043 CALL dlaset( 'Full', mnmin, mnmin, zero, one, u, ldpt )
1044 CALL dlaset( 'Full', mnmin, mnmin, zero, one, vt, ldpt )
1045*
1046 CALL dbdsdc( uplo, 'I', mnmin, s1, work, u, ldpt, vt, ldpt,
1047 $ dum, idum, work( mnmin+1 ), iwork, iinfo )
1048*
1049* Check error code from DBDSDC.
1050*
1051 IF( iinfo.NE.0 ) THEN
1052 WRITE( nout, fmt = 9998 )'DBDSDC(vects)', iinfo, m, n,
1053 $ jtype, ioldsd
1054 info = abs( iinfo )
1055 IF( iinfo.LT.0 ) THEN
1056 RETURN
1057 ELSE
1058 result( 15 ) = ulpinv
1059 GO TO 270
1060 END IF
1061 END IF
1062*
1063* Use DBDSDC to compute only the singular values of the
1064* bidiagonal matrix B; U and VT should not be modified.
1065*
1066 CALL dcopy( mnmin, bd, 1, s2, 1 )
1067 IF( mnmin.GT.0 )
1068 $ CALL dcopy( mnmin-1, be, 1, work, 1 )
1069*
1070 CALL dbdsdc( uplo, 'N', mnmin, s2, work, dum, 1, dum, 1,
1071 $ dum, idum, work( mnmin+1 ), iwork, iinfo )
1072*
1073* Check error code from DBDSDC.
1074*
1075 IF( iinfo.NE.0 ) THEN
1076 WRITE( nout, fmt = 9998 )'DBDSDC(values)', iinfo, m, n,
1077 $ jtype, ioldsd
1078 info = abs( iinfo )
1079 IF( iinfo.LT.0 ) THEN
1080 RETURN
1081 ELSE
1082 result( 18 ) = ulpinv
1083 GO TO 270
1084 END IF
1085 END IF
1086*
1087* Test 15: Check the decomposition B := U * S1 * VT
1088* 16: Check the orthogonality of U
1089* 17: Check the orthogonality of VT
1090*
1091 CALL dbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
1092 $ work, result( 15 ) )
1093 CALL dort01( 'Columns', mnmin, mnmin, u, ldpt, work, lwork,
1094 $ result( 16 ) )
1095 CALL dort01( 'Rows', mnmin, mnmin, vt, ldpt, work, lwork,
1096 $ result( 17 ) )
1097*
1098* Test 18: Check that the singular values are sorted in
1099* non-increasing order and are non-negative
1100*
1101 result( 18 ) = zero
1102 DO 150 i = 1, mnmin - 1
1103 IF( s1( i ).LT.s1( i+1 ) )
1104 $ result( 18 ) = ulpinv
1105 IF( s1( i ).LT.zero )
1106 $ result( 18 ) = ulpinv
1107 150 CONTINUE
1108 IF( mnmin.GE.1 ) THEN
1109 IF( s1( mnmin ).LT.zero )
1110 $ result( 18 ) = ulpinv
1111 END IF
1112*
1113* Test 19: Compare DBDSQR with and without singular vectors
1114*
1115 temp2 = zero
1116*
1117 DO 160 j = 1, mnmin
1118 temp1 = abs( s1( j )-s2( j ) ) /
1119 $ max( sqrt( unfl )*max( s1( 1 ), one ),
1120 $ ulp*max( abs( s1( 1 ) ), abs( s2( 1 ) ) ) )
1121 temp2 = max( temp1, temp2 )
1122 160 CONTINUE
1123*
1124 result( 19 ) = temp2
1125*
1126*
1127* Use DBDSVDX to compute the SVD of the bidiagonal matrix B:
1128* B := U * S1 * VT
1129*
1130 IF( jtype.EQ.10 .OR. jtype.EQ.16 ) THEN
1131* =================================
1132* Matrix types temporarily disabled
1133* =================================
1134 result( 20:34 ) = zero
1135 GO TO 270
1136 END IF
1137*
1138 iwbs = 1
1139 iwbd = iwbs + mnmin
1140 iwbe = iwbd + mnmin
1141 iwbz = iwbe + mnmin
1142 iwwork = iwbz + 2*mnmin*(mnmin+1)
1143 mnmin2 = max( 1,mnmin*2 )
1144*
1145 CALL dcopy( mnmin, bd, 1, work( iwbd ), 1 )
1146 IF( mnmin.GT.0 )
1147 $ CALL dcopy( mnmin-1, be, 1, work( iwbe ), 1 )
1148*
1149 CALL dbdsvdx( uplo, 'V', 'A', mnmin, work( iwbd ),
1150 $ work( iwbe ), zero, zero, 0, 0, ns1, s1,
1151 $ work( iwbz ), mnmin2, work( iwwork ),
1152 $ iwork, iinfo)
1153*
1154* Check error code from DBDSVDX.
1155*
1156 IF( iinfo.NE.0 ) THEN
1157 WRITE( nout, fmt = 9998 )'DBDSVDX(vects,A)', iinfo, m, n,
1158 $ jtype, ioldsd
1159 info = abs( iinfo )
1160 IF( iinfo.LT.0 ) THEN
1161 RETURN
1162 ELSE
1163 result( 20 ) = ulpinv
1164 GO TO 270
1165 END IF
1166 END IF
1167*
1168 j = iwbz
1169 DO 170 i = 1, ns1
1170 CALL dcopy( mnmin, work( j ), 1, u( 1,i ), 1 )
1171 j = j + mnmin
1172 CALL dcopy( mnmin, work( j ), 1, vt( i,1 ), ldpt )
1173 j = j + mnmin
1174 170 CONTINUE
1175*
1176* Use DBDSVDX to compute only the singular values of the
1177* bidiagonal matrix B; U and VT should not be modified.
1178*
1179 IF( jtype.EQ.9 ) THEN
1180* =================================
1181* Matrix types temporarily disabled
1182* =================================
1183 result( 24 ) = zero
1184 GO TO 270
1185 END IF
1186*
1187 CALL dcopy( mnmin, bd, 1, work( iwbd ), 1 )
1188 IF( mnmin.GT.0 )
1189 $ CALL dcopy( mnmin-1, be, 1, work( iwbe ), 1 )
1190*
1191 CALL dbdsvdx( uplo, 'N', 'A', mnmin, work( iwbd ),
1192 $ work( iwbe ), zero, zero, 0, 0, ns2, s2,
1193 $ work( iwbz ), mnmin2, work( iwwork ),
1194 $ iwork, iinfo )
1195*
1196* Check error code from DBDSVDX.
1197*
1198 IF( iinfo.NE.0 ) THEN
1199 WRITE( nout, fmt = 9998 )'DBDSVDX(values,A)', iinfo,
1200 $ m, n, jtype, ioldsd
1201 info = abs( iinfo )
1202 IF( iinfo.LT.0 ) THEN
1203 RETURN
1204 ELSE
1205 result( 24 ) = ulpinv
1206 GO TO 270
1207 END IF
1208 END IF
1209*
1210* Save S1 for tests 30-34.
1211*
1212 CALL dcopy( mnmin, s1, 1, work( iwbs ), 1 )
1213*
1214* Test 20: Check the decomposition B := U * S1 * VT
1215* 21: Check the orthogonality of U
1216* 22: Check the orthogonality of VT
1217* 23: Check that the singular values are sorted in
1218* non-increasing order and are non-negative
1219* 24: Compare DBDSVDX with and without singular vectors
1220*
1221 CALL dbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt,
1222 $ ldpt, work( iwbs+mnmin ), result( 20 ) )
1223 CALL dort01( 'Columns', mnmin, mnmin, u, ldpt,
1224 $ work( iwbs+mnmin ), lwork-mnmin,
1225 $ result( 21 ) )
1226 CALL dort01( 'Rows', mnmin, mnmin, vt, ldpt,
1227 $ work( iwbs+mnmin ), lwork-mnmin,
1228 $ result( 22) )
1229*
1230 result( 23 ) = zero
1231 DO 180 i = 1, mnmin - 1
1232 IF( s1( i ).LT.s1( i+1 ) )
1233 $ result( 23 ) = ulpinv
1234 IF( s1( i ).LT.zero )
1235 $ result( 23 ) = ulpinv
1236 180 CONTINUE
1237 IF( mnmin.GE.1 ) THEN
1238 IF( s1( mnmin ).LT.zero )
1239 $ result( 23 ) = ulpinv
1240 END IF
1241*
1242 temp2 = zero
1243 DO 190 j = 1, mnmin
1244 temp1 = abs( s1( j )-s2( j ) ) /
1245 $ max( sqrt( unfl )*max( s1( 1 ), one ),
1246 $ ulp*max( abs( s1( 1 ) ), abs( s2( 1 ) ) ) )
1247 temp2 = max( temp1, temp2 )
1248 190 CONTINUE
1249 result( 24 ) = temp2
1250 anorm = s1( 1 )
1251*
1252* Use DBDSVDX with RANGE='I': choose random values for IL and
1253* IU, and ask for the IL-th through IU-th singular values
1254* and corresponding vectors.
1255*
1256 DO 200 i = 1, 4
1257 iseed2( i ) = iseed( i )
1258 200 CONTINUE
1259 IF( mnmin.LE.1 ) THEN
1260 il = 1
1261 iu = mnmin
1262 ELSE
1263 il = 1 + int( ( mnmin-1 )*dlarnd( 1, iseed2 ) )
1264 iu = 1 + int( ( mnmin-1 )*dlarnd( 1, iseed2 ) )
1265 IF( iu.LT.il ) THEN
1266 itemp = iu
1267 iu = il
1268 il = itemp
1269 END IF
1270 END IF
1271*
1272 CALL dcopy( mnmin, bd, 1, work( iwbd ), 1 )
1273 IF( mnmin.GT.0 )
1274 $ CALL dcopy( mnmin-1, be, 1, work( iwbe ), 1 )
1275*
1276 CALL dbdsvdx( uplo, 'V', 'I', mnmin, work( iwbd ),
1277 $ work( iwbe ), zero, zero, il, iu, ns1, s1,
1278 $ work( iwbz ), mnmin2, work( iwwork ),
1279 $ iwork, iinfo)
1280*
1281* Check error code from DBDSVDX.
1282*
1283 IF( iinfo.NE.0 ) THEN
1284 WRITE( nout, fmt = 9998 )'DBDSVDX(vects,I)', iinfo,
1285 $ m, n, jtype, ioldsd
1286 info = abs( iinfo )
1287 IF( iinfo.LT.0 ) THEN
1288 RETURN
1289 ELSE
1290 result( 25 ) = ulpinv
1291 GO TO 270
1292 END IF
1293 END IF
1294*
1295 j = iwbz
1296 DO 210 i = 1, ns1
1297 CALL dcopy( mnmin, work( j ), 1, u( 1,i ), 1 )
1298 j = j + mnmin
1299 CALL dcopy( mnmin, work( j ), 1, vt( i,1 ), ldpt )
1300 j = j + mnmin
1301 210 CONTINUE
1302*
1303* Use DBDSVDX to compute only the singular values of the
1304* bidiagonal matrix B; U and VT should not be modified.
1305*
1306 CALL dcopy( mnmin, bd, 1, work( iwbd ), 1 )
1307 IF( mnmin.GT.0 )
1308 $ CALL dcopy( mnmin-1, be, 1, work( iwbe ), 1 )
1309*
1310 CALL dbdsvdx( uplo, 'N', 'I', mnmin, work( iwbd ),
1311 $ work( iwbe ), zero, zero, il, iu, ns2, s2,
1312 $ work( iwbz ), mnmin2, work( iwwork ),
1313 $ iwork, iinfo )
1314*
1315* Check error code from DBDSVDX.
1316*
1317 IF( iinfo.NE.0 ) THEN
1318 WRITE( nout, fmt = 9998 )'DBDSVDX(values,I)', iinfo,
1319 $ m, n, jtype, ioldsd
1320 info = abs( iinfo )
1321 IF( iinfo.LT.0 ) THEN
1322 RETURN
1323 ELSE
1324 result( 29 ) = ulpinv
1325 GO TO 270
1326 END IF
1327 END IF
1328*
1329* Test 25: Check S1 - U' * B * VT'
1330* 26: Check the orthogonality of U
1331* 27: Check the orthogonality of VT
1332* 28: Check that the singular values are sorted in
1333* non-increasing order and are non-negative
1334* 29: Compare DBDSVDX with and without singular vectors
1335*
1336 CALL dbdt04( uplo, mnmin, bd, be, s1, ns1, u,
1337 $ ldpt, vt, ldpt, work( iwbs+mnmin ),
1338 $ result( 25 ) )
1339 CALL dort01( 'Columns', mnmin, ns1, u, ldpt,
1340 $ work( iwbs+mnmin ), lwork-mnmin,
1341 $ result( 26 ) )
1342 CALL dort01( 'Rows', ns1, mnmin, vt, ldpt,
1343 $ work( iwbs+mnmin ), lwork-mnmin,
1344 $ result( 27 ) )
1345*
1346 result( 28 ) = zero
1347 DO 220 i = 1, ns1 - 1
1348 IF( s1( i ).LT.s1( i+1 ) )
1349 $ result( 28 ) = ulpinv
1350 IF( s1( i ).LT.zero )
1351 $ result( 28 ) = ulpinv
1352 220 CONTINUE
1353 IF( ns1.GE.1 ) THEN
1354 IF( s1( ns1 ).LT.zero )
1355 $ result( 28 ) = ulpinv
1356 END IF
1357*
1358 temp2 = zero
1359 DO 230 j = 1, ns1
1360 temp1 = abs( s1( j )-s2( j ) ) /
1361 $ max( sqrt( unfl )*max( s1( 1 ), one ),
1362 $ ulp*max( abs( s1( 1 ) ), abs( s2( 1 ) ) ) )
1363 temp2 = max( temp1, temp2 )
1364 230 CONTINUE
1365 result( 29 ) = temp2
1366*
1367* Use DBDSVDX with RANGE='V': determine the values VL and VU
1368* of the IL-th and IU-th singular values and ask for all
1369* singular values in this range.
1370*
1371 CALL dcopy( mnmin, work( iwbs ), 1, s1, 1 )
1372*
1373 IF( mnmin.GT.0 ) THEN
1374 IF( il.NE.1 ) THEN
1375 vu = s1( il ) + max( half*abs( s1( il )-s1( il-1 ) ),
1376 $ ulp*anorm, two*rtunfl )
1377 ELSE
1378 vu = s1( 1 ) + max( half*abs( s1( mnmin )-s1( 1 ) ),
1379 $ ulp*anorm, two*rtunfl )
1380 END IF
1381 IF( iu.NE.ns1 ) THEN
1382 vl = s1( iu ) - max( ulp*anorm, two*rtunfl,
1383 $ half*abs( s1( iu+1 )-s1( iu ) ) )
1384 ELSE
1385 vl = s1( ns1 ) - max( ulp*anorm, two*rtunfl,
1386 $ half*abs( s1( mnmin )-s1( 1 ) ) )
1387 END IF
1388 vl = max( vl,zero )
1389 vu = max( vu,zero )
1390 IF( vl.GE.vu ) vu = max( vu*2, vu+vl+half )
1391 ELSE
1392 vl = zero
1393 vu = one
1394 END IF
1395*
1396 CALL dcopy( mnmin, bd, 1, work( iwbd ), 1 )
1397 IF( mnmin.GT.0 )
1398 $ CALL dcopy( mnmin-1, be, 1, work( iwbe ), 1 )
1399*
1400 CALL dbdsvdx( uplo, 'V', 'V', mnmin, work( iwbd ),
1401 $ work( iwbe ), vl, vu, 0, 0, ns1, s1,
1402 $ work( iwbz ), mnmin2, work( iwwork ),
1403 $ iwork, iinfo )
1404*
1405* Check error code from DBDSVDX.
1406*
1407 IF( iinfo.NE.0 ) THEN
1408 WRITE( nout, fmt = 9998 )'DBDSVDX(vects,V)', iinfo,
1409 $ m, n, jtype, ioldsd
1410 info = abs( iinfo )
1411 IF( iinfo.LT.0 ) THEN
1412 RETURN
1413 ELSE
1414 result( 30 ) = ulpinv
1415 GO TO 270
1416 END IF
1417 END IF
1418*
1419 j = iwbz
1420 DO 240 i = 1, ns1
1421 CALL dcopy( mnmin, work( j ), 1, u( 1,i ), 1 )
1422 j = j + mnmin
1423 CALL dcopy( mnmin, work( j ), 1, vt( i,1 ), ldpt )
1424 j = j + mnmin
1425 240 CONTINUE
1426*
1427* Use DBDSVDX to compute only the singular values of the
1428* bidiagonal matrix B; U and VT should not be modified.
1429*
1430 CALL dcopy( mnmin, bd, 1, work( iwbd ), 1 )
1431 IF( mnmin.GT.0 )
1432 $ CALL dcopy( mnmin-1, be, 1, work( iwbe ), 1 )
1433*
1434 CALL dbdsvdx( uplo, 'N', 'V', mnmin, work( iwbd ),
1435 $ work( iwbe ), vl, vu, 0, 0, ns2, s2,
1436 $ work( iwbz ), mnmin2, work( iwwork ),
1437 $ iwork, iinfo )
1438*
1439* Check error code from DBDSVDX.
1440*
1441 IF( iinfo.NE.0 ) THEN
1442 WRITE( nout, fmt = 9998 )'DBDSVDX(values,V)', iinfo,
1443 $ m, n, jtype, ioldsd
1444 info = abs( iinfo )
1445 IF( iinfo.LT.0 ) THEN
1446 RETURN
1447 ELSE
1448 result( 34 ) = ulpinv
1449 GO TO 270
1450 END IF
1451 END IF
1452*
1453* Test 30: Check S1 - U' * B * VT'
1454* 31: Check the orthogonality of U
1455* 32: Check the orthogonality of VT
1456* 33: Check that the singular values are sorted in
1457* non-increasing order and are non-negative
1458* 34: Compare DBDSVDX with and without singular vectors
1459*
1460 CALL dbdt04( uplo, mnmin, bd, be, s1, ns1, u,
1461 $ ldpt, vt, ldpt, work( iwbs+mnmin ),
1462 $ result( 30 ) )
1463 CALL dort01( 'Columns', mnmin, ns1, u, ldpt,
1464 $ work( iwbs+mnmin ), lwork-mnmin,
1465 $ result( 31 ) )
1466 CALL dort01( 'Rows', ns1, mnmin, vt, ldpt,
1467 $ work( iwbs+mnmin ), lwork-mnmin,
1468 $ result( 32 ) )
1469*
1470 result( 33 ) = zero
1471 DO 250 i = 1, ns1 - 1
1472 IF( s1( i ).LT.s1( i+1 ) )
1473 $ result( 28 ) = ulpinv
1474 IF( s1( i ).LT.zero )
1475 $ result( 28 ) = ulpinv
1476 250 CONTINUE
1477 IF( ns1.GE.1 ) THEN
1478 IF( s1( ns1 ).LT.zero )
1479 $ result( 28 ) = ulpinv
1480 END IF
1481*
1482 temp2 = zero
1483 DO 260 j = 1, ns1
1484 temp1 = abs( s1( j )-s2( j ) ) /
1485 $ max( sqrt( unfl )*max( s1( 1 ), one ),
1486 $ ulp*max( abs( s1( 1 ) ), abs( s2( 1 ) ) ) )
1487 temp2 = max( temp1, temp2 )
1488 260 CONTINUE
1489 result( 34 ) = temp2
1490*
1491* End of Loop -- Check for RESULT(j) > THRESH
1492*
1493 270 CONTINUE
1494*
1495 DO 280 j = 1, 34
1496 IF( result( j ).GE.thresh ) THEN
1497 IF( nfail.EQ.0 )
1498 $ CALL dlahd2( nout, path )
1499 WRITE( nout, fmt = 9999 )m, n, jtype, ioldsd, j,
1500 $ result( j )
1501 nfail = nfail + 1
1502 END IF
1503 280 CONTINUE
1504 IF( .NOT.bidiag ) THEN
1505 ntest = ntest + 34
1506 ELSE
1507 ntest = ntest + 30
1508 END IF
1509*
1510 290 CONTINUE
1511 300 CONTINUE
1512*
1513* Summary
1514*
1515 CALL alasum( path, nout, nfail, ntest, 0 )
1516*
1517 RETURN
1518*
1519* End of DCHKBD
1520*
1521 9999 FORMAT( ' M=', i5, ', N=', i5, ', type ', i2, ', seed=',
1522 $ 4( i4, ',' ), ' test(', i2, ')=', g11.4 )
1523 9998 FORMAT( ' DCHKBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1524 $ i6, ', N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
1525 $ i5, ')' )
1526*
1527 END
subroutine alasum(type, nout, nfail, nrun, nerrs)
ALASUM
Definition alasum.f:73
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dbdt01(m, n, kd, a, lda, q, ldq, d, e, pt, ldpt, work, resid)
DBDT01
Definition dbdt01.f:141
subroutine dbdt02(m, n, b, ldb, c, ldc, u, ldu, work, resid)
DBDT02
Definition dbdt02.f:112
subroutine dbdt03(uplo, n, kd, d, e, u, ldu, s, vt, ldvt, work, resid)
DBDT03
Definition dbdt03.f:135
subroutine dbdt04(uplo, n, d, e, s, ns, u, ldu, vt, ldvt, work, resid)
DBDT04
Definition dbdt04.f:131
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:493
subroutine dlahd2(iounit, path)
DLAHD2
Definition dlahd2.f:65
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:471
subroutine dlatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
DLATMS
Definition dlatms.f:321
subroutine dort01(rowcol, m, n, u, ldu, work, lwork, resid)
DORT01
Definition dort01.f:116
subroutine dbdsdc(uplo, compq, n, d, e, u, ldu, vt, ldvt, q, iq, work, iwork, info)
DBDSDC
Definition dbdsdc.f:198
subroutine dbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, work, info)
DBDSQR
Definition dbdsqr.f:241
subroutine dbdsvdx(uplo, jobz, range, n, d, e, vl, vu, il, iu, ns, s, z, ldz, work, iwork, info)
DBDSVDX
Definition dbdsvdx.f:226
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
subroutine dgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
DGEBRD
Definition dgebrd.f:205
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
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:110
subroutine dorgbr(vect, m, n, k, a, lda, tau, work, lwork, info)
DORGBR
Definition dorgbr.f:157