LAPACK 3.12.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ schkbd()

 subroutine schkbd ( integer nsizes, integer, dimension( * ) mval, integer, dimension( * ) nval, integer ntypes, logical, dimension( * ) dotype, integer nrhs, integer, dimension( 4 ) iseed, real thresh, real, dimension( lda, * ) a, integer lda, real, dimension( * ) bd, real, dimension( * ) be, real, dimension( * ) s1, real, dimension( * ) s2, real, dimension( ldx, * ) x, integer ldx, real, dimension( ldx, * ) y, real, dimension( ldx, * ) z, real, dimension( ldq, * ) q, integer ldq, real, dimension( ldpt, * ) pt, integer ldpt, real, dimension( ldpt, * ) u, real, dimension( ldpt, * ) vt, real, dimension( * ) work, integer lwork, integer, dimension( * ) iwork, integer nout, integer info )

SCHKBD

Purpose:
``` SCHKBD checks the singular value decomposition (SVD) routines.

SGEBRD reduces a real general m by n matrix A to upper or lower
bidiagonal form B by an orthogonal transformation:  Q' * A * P = B
(or A = Q * B * P').  The matrix B is upper bidiagonal if m >= n
and lower bidiagonal if m < n.

SORGBR generates the orthogonal matrices Q and P' from SGEBRD.
Note that Q and P are not necessarily square.

SBDSQR computes the singular value decomposition of the bidiagonal
matrix B as B = U S V'.  It is called three times to compute
1)  B = U S1 V', where S1 is the diagonal matrix of singular
values and the columns of the matrices U and V are the left
and right singular vectors, respectively, of B.
2)  Same as 1), but the singular values are stored in S2 and the
singular vectors are not computed.
3)  A = (UQ) S (P'V'), the SVD of the original matrix A.
In addition, SBDSQR has an option to apply the left orthogonal matrix
U to a matrix X, useful in least squares applications.

SBDSDC computes the singular value decomposition of the bidiagonal
matrix B as B = U S V' using divide-and-conquer. It is called twice
to compute
1) B = U S1 V', where S1 is the diagonal matrix of singular
values and the columns of the matrices U and V are the left
and right singular vectors, respectively, of B.
2) Same as 1), but the singular values are stored in S2 and the
singular vectors are not computed.

SBDSVDX computes the singular value decomposition of the bidiagonal
matrix B as B = U S V' using bisection and inverse iteration. It is
called six times to compute
1) B = U S1 V', RANGE='A', where S1 is the diagonal matrix of singular
values and the columns of the matrices U and V are the left
and right singular vectors, respectively, of B.
2) Same as 1), but the singular values are stored in S2 and the
singular vectors are not computed.
3) B = U S1 V', RANGE='I', with where S1 is the diagonal matrix of singular
values and the columns of the matrices U and V are the left
and right singular vectors, respectively, of B
4) Same as 3), but the singular values are stored in S2 and the
singular vectors are not computed.
5) B = U S1 V', RANGE='V', with where S1 is the diagonal matrix of singular
values and the columns of the matrices U and V are the left
and right singular vectors, respectively, of B
6) Same as 5), but the singular values are stored in S2 and the
singular vectors are not computed.

For each pair of matrix dimensions (M,N) and each selected matrix
type, an M by N matrix A and an M by NRHS matrix X are generated.
The problem dimensions are as follows
A:          M x N
Q:          M x min(M,N) (but M x M if NRHS > 0)
P:          min(M,N) x N
B:          min(M,N) x min(M,N)
U, V:       min(M,N) x min(M,N)
S1, S2      diagonal, order min(M,N)
X:          M x NRHS

For each generated matrix, 14 tests are performed:

Test SGEBRD and SORGBR

(1)   | A - Q B PT | / ( |A| max(M,N) ulp ), PT = P'

(2)   | I - Q' Q | / ( M ulp )

(3)   | I - PT PT' | / ( N ulp )

Test SBDSQR on bidiagonal matrix B

(4)   | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V'

(5)   | Y - U Z | / ( |Y| max(min(M,N),k) ulp ), where Y = Q' X
and   Z = U' Y.
(6)   | I - U' U | / ( min(M,N) ulp )

(7)   | I - VT VT' | / ( min(M,N) ulp )

(8)   S1 contains min(M,N) nonnegative values in decreasing order.
(Return 0 if true, 1/ULP if false.)

(9)   | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
computing U and V.

(10)  0 if the true singular values of B are within THRESH of
those in S1.  2*THRESH if they are not.  (Tested using
SSVDCH)

Test SBDSQR on matrix A

(11)  | A - (QU) S (VT PT) | / ( |A| max(M,N) ulp )

(12)  | X - (QU) Z | / ( |X| max(M,k) ulp )

(13)  | I - (QU)'(QU) | / ( M ulp )

(14)  | I - (VT PT) (PT'VT') | / ( N ulp )

Test SBDSDC on bidiagonal matrix B

(15)  | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V'

(16)  | I - U' U | / ( min(M,N) ulp )

(17)  | I - VT VT' | / ( min(M,N) ulp )

(18)  S1 contains min(M,N) nonnegative values in decreasing order.
(Return 0 if true, 1/ULP if false.)

(19)  | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
computing U and V.
Test SBDSVDX on bidiagonal matrix B

(20)  | B - U S1 VT | / ( |B| min(M,N) ulp ), VT = V'

(21)  | I - U' U | / ( min(M,N) ulp )

(22)  | I - VT VT' | / ( min(M,N) ulp )

(23)  S1 contains min(M,N) nonnegative values in decreasing order.
(Return 0 if true, 1/ULP if false.)

(24)  | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
computing U and V.

(25)  | S1 - U' B VT' | / ( |S| n ulp )    SBDSVDX('V', 'I')

(26)  | I - U' U | / ( min(M,N) ulp )

(27)  | I - VT VT' | / ( min(M,N) ulp )

(28)  S1 contains min(M,N) nonnegative values in decreasing order.
(Return 0 if true, 1/ULP if false.)

(29)  | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
computing U and V.

(30)  | S1 - U' B VT' | / ( |S1| n ulp )   SBDSVDX('V', 'V')

(31)  | I - U' U | / ( min(M,N) ulp )

(32)  | I - VT VT' | / ( min(M,N) ulp )

(33)  S1 contains min(M,N) nonnegative values in decreasing order.
(Return 0 if true, 1/ULP if false.)

(34)  | S1 - S2 | / ( |S1| ulp ), where S2 is computed without
computing U and V.

The possible matrix types are

(1)  The zero matrix.
(2)  The identity matrix.

(3)  A diagonal matrix with evenly spaced entries
1, ..., ULP  and random signs.
(ULP = (first number larger than 1) - 1 )
(4)  A diagonal matrix with geometrically spaced entries
1, ..., ULP  and random signs.
(5)  A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
and random signs.

(6)  Same as (3), but multiplied by SQRT( overflow threshold )
(7)  Same as (3), but multiplied by SQRT( underflow threshold )

(8)  A matrix of the form  U D V, where U and V are orthogonal and
D has evenly spaced entries 1, ..., ULP with random signs
on the diagonal.

(9)  A matrix of the form  U D V, where U and V are orthogonal and
D has geometrically spaced entries 1, ..., ULP with random
signs on the diagonal.

(10) A matrix of the form  U D V, where U and V are orthogonal and
D has "clustered" entries 1, ULP,..., ULP with random
signs on the diagonal.

(11) Same as (8), but multiplied by SQRT( overflow threshold )
(12) Same as (8), but multiplied by SQRT( underflow threshold )

(13) Rectangular matrix with random entries chosen from (-1,1).
(14) Same as (13), but multiplied by SQRT( overflow threshold )
(15) Same as (13), but multiplied by SQRT( underflow threshold )

Special case:
(16) A bidiagonal matrix with random entries chosen from a
logarithmic distribution on [ulp^2,ulp^(-2)]  (I.e., each
entry is  e^x, where x is chosen uniformly on
[ 2 log(ulp), -2 log(ulp) ] .)  For *this* type:
(a) SGEBRD is not called to reduce it to bidiagonal form.
(b) the bidiagonal is  min(M,N) x min(M,N); if M<N, the
matrix will be lower bidiagonal, otherwise upper.
(c) only tests 5--8 and 14 are performed.

A subset of the full set of matrix types may be selected through
the logical array DOTYPE.```
Parameters
 [in] NSIZES ``` NSIZES is INTEGER The number of values of M and N contained in the vectors MVAL and NVAL. The matrix sizes are used in pairs (M,N).``` [in] MVAL ``` MVAL is INTEGER array, dimension (NM) The values of the matrix row dimension M.``` [in] NVAL ``` NVAL is INTEGER array, dimension (NM) The values of the matrix column dimension N.``` [in] NTYPES ``` NTYPES is INTEGER The number of elements in DOTYPE. If it is zero, SCHKBD does nothing. It must be at least zero. If it is MAXTYP+1 and NSIZES is 1, then an additional type, MAXTYP+1 is defined, which is to use whatever matrices are in A and B. This is only useful if DOTYPE(1:MAXTYP) is .FALSE. and DOTYPE(MAXTYP+1) is .TRUE. .``` [in] DOTYPE ``` DOTYPE is LOGICAL array, dimension (NTYPES) If DOTYPE(j) is .TRUE., then for each size (m,n), a matrix of type j will be generated. If NTYPES is smaller than the maximum number of types defined (PARAMETER MAXTYP), then types NTYPES+1 through MAXTYP will not be generated. If NTYPES is larger than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES) will be ignored.``` [in] NRHS ``` NRHS is INTEGER The number of columns in the "right-hand side" matrices X, Y, and Z, used in testing SBDSQR. If NRHS = 0, then the operations on the right-hand side will not be tested. NRHS must be at least 0.``` [in,out] ISEED ``` ISEED is INTEGER array, dimension (4) On entry ISEED specifies the seed of the random number generator. The array elements should be between 0 and 4095; if not they will be reduced mod 4096. Also, ISEED(4) must be odd. The values of ISEED are changed on exit, and can be used in the next call to SCHKBD to continue the same random number sequence.``` [in] THRESH ``` THRESH is REAL The threshold value for the test ratios. A result is included in the output file if RESULT >= THRESH. To have every test ratio printed, use THRESH = 0. Note that the expected value of the test ratios is O(1), so THRESH should be a reasonably small multiple of 1, e.g., 10 or 100.``` [out] A ``` A is REAL array, dimension (LDA,NMAX) where NMAX is the maximum value of N in NVAL.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,MMAX), where MMAX is the maximum value of M in MVAL.``` [out] BD ``` BD is REAL array, dimension (max(min(MVAL(j),NVAL(j))))``` [out] BE ``` BE is REAL array, dimension (max(min(MVAL(j),NVAL(j))))``` [out] S1 ``` S1 is REAL array, dimension (max(min(MVAL(j),NVAL(j))))``` [out] S2 ``` S2 is REAL array, dimension (max(min(MVAL(j),NVAL(j))))``` [out] X ` X is REAL array, dimension (LDX,NRHS)` [in] LDX ``` LDX is INTEGER The leading dimension of the arrays X, Y, and Z. LDX >= max(1,MMAX)``` [out] Y ` Y is REAL array, dimension (LDX,NRHS)` [out] Z ` Z is REAL array, dimension (LDX,NRHS)` [out] Q ` Q is REAL array, dimension (LDQ,MMAX)` [in] LDQ ``` LDQ is INTEGER The leading dimension of the array Q. LDQ >= max(1,MMAX).``` [out] PT ` PT is REAL array, dimension (LDPT,NMAX)` [in] LDPT ``` LDPT is INTEGER The leading dimension of the arrays PT, U, and V. LDPT >= max(1, max(min(MVAL(j),NVAL(j)))).``` [out] U ``` U is REAL array, dimension (LDPT,max(min(MVAL(j),NVAL(j))))``` [out] VT ``` VT is REAL array, dimension (LDPT,max(min(MVAL(j),NVAL(j))))``` [out] WORK ` WORK is REAL array, dimension (LWORK)` [in] LWORK ``` LWORK is INTEGER The number of entries in WORK. This must be at least 3(M+N) and M(M + max(M,N,k) + 1) + N*min(M,N) for all pairs (M,N)=(MM(j),NN(j))``` [out] IWORK ` IWORK is INTEGER array, dimension at least 8*min(M,N)` [in] NOUT ``` NOUT is INTEGER The FORTRAN unit number for printing out error messages (e.g., if a routine returns IINFO not equal to 0.)``` [out] INFO ``` INFO is INTEGER If 0, then everything ran OK. -1: NSIZES < 0 -2: Some MM(j) < 0 -3: Some NN(j) < 0 -4: NTYPES < 0 -6: NRHS < 0 -8: THRESH < 0 -11: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ). -17: LDB < 1 or LDB < MMAX. -21: LDQ < 1 or LDQ < MMAX. -23: LDPT< 1 or LDPT< MNMAX. -27: LWORK too small. If SLATMR, SLATMS, SGEBRD, SORGBR, or SBDSQR, returns an error code, the absolute value of it is returned. ----------------------------------------------------------------------- Some Local Variables and Parameters: ---- ----- --------- --- ---------- ZERO, ONE Real 0 and 1. MAXTYP The number of types defined. NTEST The number of tests performed, or which can be performed so far, for the current matrix. MMAX Largest value in NN. NMAX Largest value in NN. MNMIN min(MM(j), NN(j)) (the dimension of the bidiagonal matrix.) MNMAX The maximum value of MNMIN for j=1,...,NSIZES. NFAIL The number of tests which have exceeded THRESH COND, IMODE Values to be passed to the matrix generators. ANORM Norm of A; passed to matrix generators. OVFL, UNFL Overflow and underflow thresholds. RTOVFL, RTUNFL Square roots of the previous 2 values. ULP, ULPINV Finest relative precision and its inverse. The following four arrays decode JTYPE: KTYPE(j) The general type (1-10) for type "j". KMODE(j) The MODE value to be passed to the matrix generator for type "j". KMAGN(j) The order of magnitude ( O(1), O(overflow^(1/2) ), O(underflow^(1/2) )```

Definition at line 489 of file schkbd.f.

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 REAL THRESH
502* ..
503* .. Array Arguments ..
504 LOGICAL DOTYPE( * )
505 INTEGER ISEED( 4 ), IWORK( * ), MVAL( * ), NVAL( * )
506 REAL 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 REAL ZERO, ONE, TWO, HALF
516 parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
517 \$ half = 0.5e0 )
518 INTEGER MAXTYP
519 parameter( maxtyp = 16 )
520* ..
521* .. Local Scalars ..
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 REAL 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 REAL DUM( 1 ), DUMMA( 1 ), RESULT( 40 )
539* ..
540* .. External Functions ..
541 REAL SLAMCH, SLARND, SSXT1
542 EXTERNAL slamch, slarnd, ssxt1
543* ..
544* .. External Subroutines ..
545 EXTERNAL alasum, sbdsdc, sbdsqr, sbdsvdx, sbdt01,
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*
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 )
584 nmax = max( nmax, nval( j ) )
585 IF( nval( j ).LT.0 )
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( 'SCHKBD', -info )
619 RETURN
620 END IF
621*
622* Initialize constants
623*
624 path( 1: 1 ) = 'Single precision'
625 path( 2: 3 ) = 'BD'
626 nfail = 0
627 ntest = 0
628 unfl = slamch( 'Safe minimum' )
629 ovfl = slamch( 'Overflow' )
630 ulp = slamch( '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 slaset( '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 slatms( 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 slatms( 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 slatms( 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 slatmr( 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 slatmr( 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 slatmr( 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*slarnd( 2, iseed ) )
783 IF( j.LT.mnmin )
784 \$ be( j ) = exp( temp1*slarnd( 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 slatmr( 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 slatmr( 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 SGEBRD and SORGBR 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 slacpy( ' ', m, n, a, lda, q, ldq )
836 CALL sgebrd( 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 SGEBRD.
840*
841 IF( iinfo.NE.0 ) THEN
842 WRITE( nout, fmt = 9998 )'SGEBRD', iinfo, m, n,
843 \$ jtype, ioldsd
844 info = abs( iinfo )
845 RETURN
846 END IF
847*
848 CALL slacpy( ' ', 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 sorgbr( 'Q', m, mq, n, q, ldq, work,
861 \$ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
862*
863* Check error code from SORGBR.
864*
865 IF( iinfo.NE.0 ) THEN
866 WRITE( nout, fmt = 9998 )'SORGBR(Q)', iinfo, m, n,
867 \$ jtype, ioldsd
868 info = abs( iinfo )
869 RETURN
870 END IF
871*
872* Generate P'
873*
874 CALL sorgbr( 'P', mnmin, n, m, pt, ldpt, work( mnmin+1 ),
875 \$ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
876*
877* Check error code from SORGBR.
878*
879 IF( iinfo.NE.0 ) THEN
880 WRITE( nout, fmt = 9998 )'SORGBR(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 sgemm( '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 sbdt01( m, n, 1, a, lda, q, ldq, bd, be, pt, ldpt,
896 \$ work, result( 1 ) )
897 CALL sort01( 'Columns', m, mq, q, ldq, work, lwork,
898 \$ result( 2 ) )
899 CALL sort01( 'Rows', mnmin, n, pt, ldpt, work, lwork,
900 \$ result( 3 ) )
901 END IF
902*
903* Use SBDSQR to form the SVD of the bidiagonal matrix B:
904* B := U * S1 * VT, and compute Z = U' * Y.
905*
906 CALL scopy( mnmin, bd, 1, s1, 1 )
907 IF( mnmin.GT.0 )
908 \$ CALL scopy( mnmin-1, be, 1, work, 1 )
909 CALL slacpy( ' ', m, nrhs, y, ldx, z, ldx )
910 CALL slaset( 'Full', mnmin, mnmin, zero, one, u, ldpt )
911 CALL slaset( 'Full', mnmin, mnmin, zero, one, vt, ldpt )
912*
913 CALL sbdsqr( uplo, mnmin, mnmin, mnmin, nrhs, s1, work, vt,
914 \$ ldpt, u, ldpt, z, ldx, work( mnmin+1 ), iinfo )
915*
916* Check error code from SBDSQR.
917*
918 IF( iinfo.NE.0 ) THEN
919 WRITE( nout, fmt = 9998 )'SBDSQR(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 SBDSQR to compute only the singular values of the
931* bidiagonal matrix B; U, VT, and Z should not be modified.
932*
933 CALL scopy( mnmin, bd, 1, s2, 1 )
934 IF( mnmin.GT.0 )
935 \$ CALL scopy( mnmin-1, be, 1, work, 1 )
936*
937 CALL sbdsqr( uplo, mnmin, 0, 0, 0, s2, work, vt, ldpt, u,
938 \$ ldpt, z, ldx, work( mnmin+1 ), iinfo )
939*
940* Check error code from SBDSQR.
941*
942 IF( iinfo.NE.0 ) THEN
943 WRITE( nout, fmt = 9998 )'SBDSQR(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 sbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
960 \$ work, result( 4 ) )
961 CALL sbdt02( mnmin, nrhs, y, ldx, z, ldx, u, ldpt, work,
962 \$ result( 5 ) )
963 CALL sort01( 'Columns', mnmin, mnmin, u, ldpt, work, lwork,
964 \$ result( 6 ) )
965 CALL sort01( '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 SBDSQR 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 SSVDCH( 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 SBDSQR 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 scopy( mnmin, bd, 1, s2, 1 )
1016 IF( mnmin.GT.0 )
1017 \$ CALL scopy( mnmin-1, be, 1, work, 1 )
1018*
1019 CALL sbdsqr( 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 sbdt01( m, n, 0, a, lda, q, ldq, s2, dumma, pt,
1028 \$ ldpt, work, result( 11 ) )
1029 CALL sbdt02( m, nrhs, x, ldx, y, ldx, q, ldq, work,
1030 \$ result( 12 ) )
1031 CALL sort01( 'Columns', m, mq, q, ldq, work, lwork,
1032 \$ result( 13 ) )
1033 CALL sort01( 'Rows', mnmin, n, pt, ldpt, work, lwork,
1034 \$ result( 14 ) )
1035 END IF
1036*
1037* Use SBDSDC to form the SVD of the bidiagonal matrix B:
1038* B := U * S1 * VT
1039*
1040 CALL scopy( mnmin, bd, 1, s1, 1 )
1041 IF( mnmin.GT.0 )
1042 \$ CALL scopy( mnmin-1, be, 1, work, 1 )
1043 CALL slaset( 'Full', mnmin, mnmin, zero, one, u, ldpt )
1044 CALL slaset( 'Full', mnmin, mnmin, zero, one, vt, ldpt )
1045*
1046 CALL sbdsdc( uplo, 'I', mnmin, s1, work, u, ldpt, vt, ldpt,
1047 \$ dum, idum, work( mnmin+1 ), iwork, iinfo )
1048*
1049* Check error code from SBDSDC.
1050*
1051 IF( iinfo.NE.0 ) THEN
1052 WRITE( nout, fmt = 9998 )'SBDSDC(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 SBDSDC to compute only the singular values of the
1064* bidiagonal matrix B; U and VT should not be modified.
1065*
1066 CALL scopy( mnmin, bd, 1, s2, 1 )
1067 IF( mnmin.GT.0 )
1068 \$ CALL scopy( mnmin-1, be, 1, work, 1 )
1069*
1070 CALL sbdsdc( uplo, 'N', mnmin, s2, work, dum, 1, dum, 1,
1071 \$ dum, idum, work( mnmin+1 ), iwork, iinfo )
1072*
1073* Check error code from SBDSDC.
1074*
1075 IF( iinfo.NE.0 ) THEN
1076 WRITE( nout, fmt = 9998 )'SBDSDC(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 sbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
1092 \$ work, result( 15 ) )
1093 CALL sort01( 'Columns', mnmin, mnmin, u, ldpt, work, lwork,
1094 \$ result( 16 ) )
1095 CALL sort01( '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 SBDSQR 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 SBDSVDX 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 scopy( mnmin, bd, 1, work( iwbd ), 1 )
1146 IF( mnmin.GT.0 )
1147 \$ CALL scopy( mnmin-1, be, 1, work( iwbe ), 1 )
1148*
1149 CALL sbdsvdx( 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 SBDSVDX.
1155*
1156 IF( iinfo.NE.0 ) THEN
1157 WRITE( nout, fmt = 9998 )'SBDSVDX(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 scopy( mnmin, work( j ), 1, u( 1,i ), 1 )
1171 j = j + mnmin
1172 CALL scopy( mnmin, work( j ), 1, vt( i,1 ), ldpt )
1173 j = j + mnmin
1174 170 CONTINUE
1175*
1176* Use SBDSVDX 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 scopy( mnmin, bd, 1, work( iwbd ), 1 )
1188 IF( mnmin.GT.0 )
1189 \$ CALL scopy( mnmin-1, be, 1, work( iwbe ), 1 )
1190*
1191 CALL sbdsvdx( 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 SBDSVDX.
1197*
1198 IF( iinfo.NE.0 ) THEN
1199 WRITE( nout, fmt = 9998 )'SBDSVDX(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 scopy( 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 SBDSVDX with and without singular vectors
1220*
1221 CALL sbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt,
1222 \$ ldpt, work( iwbs+mnmin ), result( 20 ) )
1223 CALL sort01( 'Columns', mnmin, mnmin, u, ldpt,
1224 \$ work( iwbs+mnmin ), lwork-mnmin,
1225 \$ result( 21 ) )
1226 CALL sort01( '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 SBDSVDX 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 )*slarnd( 1, iseed2 ) )
1264 iu = 1 + int( ( mnmin-1 )*slarnd( 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 scopy( mnmin, bd, 1, work( iwbd ), 1 )
1273 IF( mnmin.GT.0 )
1274 \$ CALL scopy( mnmin-1, be, 1, work( iwbe ), 1 )
1275*
1276 CALL sbdsvdx( 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 SBDSVDX.
1282*
1283 IF( iinfo.NE.0 ) THEN
1284 WRITE( nout, fmt = 9998 )'SBDSVDX(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 scopy( mnmin, work( j ), 1, u( 1,i ), 1 )
1298 j = j + mnmin
1299 CALL scopy( mnmin, work( j ), 1, vt( i,1 ), ldpt )
1300 j = j + mnmin
1301 210 CONTINUE
1302*
1303* Use SBDSVDX to compute only the singular values of the
1304* bidiagonal matrix B; U and VT should not be modified.
1305*
1306 CALL scopy( mnmin, bd, 1, work( iwbd ), 1 )
1307 IF( mnmin.GT.0 )
1308 \$ CALL scopy( mnmin-1, be, 1, work( iwbe ), 1 )
1309*
1310 CALL sbdsvdx( 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 SBDSVDX.
1316*
1317 IF( iinfo.NE.0 ) THEN
1318 WRITE( nout, fmt = 9998 )'SBDSVDX(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 SBDSVDX with and without singular vectors
1335*
1336 CALL sbdt04( uplo, mnmin, bd, be, s1, ns1, u,
1337 \$ ldpt, vt, ldpt, work( iwbs+mnmin ),
1338 \$ result( 25 ) )
1339 CALL sort01( 'Columns', mnmin, ns1, u, ldpt,
1340 \$ work( iwbs+mnmin ), lwork-mnmin,
1341 \$ result( 26 ) )
1342 CALL sort01( '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 SBDSVDX 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 scopy( 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 scopy( mnmin, bd, 1, work( iwbd ), 1 )
1397 IF( mnmin.GT.0 )
1398 \$ CALL scopy( mnmin-1, be, 1, work( iwbe ), 1 )
1399*
1400 CALL sbdsvdx( 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 SBDSVDX.
1406*
1407 IF( iinfo.NE.0 ) THEN
1408 WRITE( nout, fmt = 9998 )'SBDSVDX(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 scopy( mnmin, work( j ), 1, u( 1,i ), 1 )
1422 j = j + mnmin
1423 CALL scopy( mnmin, work( j ), 1, vt( i,1 ), ldpt )
1424 j = j + mnmin
1425 240 CONTINUE
1426*
1427* Use SBDSVDX to compute only the singular values of the
1428* bidiagonal matrix B; U and VT should not be modified.
1429*
1430 CALL scopy( mnmin, bd, 1, work( iwbd ), 1 )
1431 IF( mnmin.GT.0 )
1432 \$ CALL scopy( mnmin-1, be, 1, work( iwbe ), 1 )
1433*
1434 CALL sbdsvdx( 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 SBDSVDX.
1440*
1441 IF( iinfo.NE.0 ) THEN
1442 WRITE( nout, fmt = 9998 )'SBDSVDX(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 SBDSVDX with and without singular vectors
1459*
1460 CALL sbdt04( uplo, mnmin, bd, be, s1, ns1, u,
1461 \$ ldpt, vt, ldpt, work( iwbs+mnmin ),
1462 \$ result( 30 ) )
1463 CALL sort01( 'Columns', mnmin, ns1, u, ldpt,
1464 \$ work( iwbs+mnmin ), lwork-mnmin,
1465 \$ result( 31 ) )
1466 CALL sort01( '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 slahd2( 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 SCHKBD
1520*
1521 9999 FORMAT( ' M=', i5, ', N=', i5, ', type ', i2, ', seed=',
1522 \$ 4( i4, ',' ), ' test(', i2, ')=', g11.4 )
1523 9998 FORMAT( ' SCHKBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1524 \$ i6, ', N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
1525 \$ i5, ')' )
1526*
subroutine alasum(type, nout, nfail, nrun, nerrs)
ALASUM
Definition alasum.f:73
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine sbdsdc(uplo, compq, n, d, e, u, ldu, vt, ldvt, q, iq, work, iwork, info)
SBDSDC
Definition sbdsdc.f:198
subroutine sbdsqr(uplo, n, ncvt, nru, ncc, d, e, vt, ldvt, u, ldu, c, ldc, work, info)
SBDSQR
Definition sbdsqr.f:240
subroutine sbdsvdx(uplo, jobz, range, n, d, e, vl, vu, il, iu, ns, s, z, ldz, work, iwork, info)
SBDSVDX
Definition sbdsvdx.f:226
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
subroutine sgebrd(m, n, a, lda, d, e, tauq, taup, work, lwork, info)
SGEBRD
Definition sgebrd.f:205
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
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:110
subroutine sorgbr(vect, m, n, k, a, lda, tau, work, lwork, info)
SORGBR
Definition sorgbr.f:157
subroutine sbdt01(m, n, kd, a, lda, q, ldq, d, e, pt, ldpt, work, resid)
SBDT01
Definition sbdt01.f:141
subroutine sbdt02(m, n, b, ldb, c, ldc, u, ldu, work, resid)
SBDT02
Definition sbdt02.f:112
subroutine sbdt03(uplo, n, kd, d, e, u, ldu, s, vt, ldvt, work, resid)
SBDT03
Definition sbdt03.f:135
subroutine sbdt04(uplo, n, d, e, s, ns, u, ldu, vt, ldvt, work, resid)
SBDT04
Definition sbdt04.f:131
subroutine slahd2(iounit, path)
SLAHD2
Definition slahd2.f:65
real function slarnd(idist, iseed)
SLARND
Definition slarnd.f:73
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:471
subroutine slatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
SLATMS
Definition slatms.f:321
subroutine sort01(rowcol, m, n, u, ldu, work, lwork, resid)
SORT01
Definition sort01.f:116
real function ssxt1(ijob, d1, n1, d2, n2, abstol, ulp, unfl)
SSXT1
Definition ssxt1.f:106
Here is the call graph for this function:
Here is the caller graph for this function: