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