LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
schkbd.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine schkbd (NSIZES, MVAL, NVAL, NTYPES, DOTYPE, NRHS, ISEED, THRESH, A, LDA, BD, BE, S1, S2, X, LDX, Y, Z, Q, LDQ, PT, LDPT, U, VT, WORK, LWORK, IWORK, NOUT, INFO)
 SCHKBD More...
 

Function/Subroutine Documentation

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.

 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.
 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) )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 438 of file schkbd.f.

438 *
439 * -- LAPACK test routine (version 3.4.0) --
440 * -- LAPACK is a software package provided by Univ. of Tennessee, --
441 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
442 * November 2011
443 *
444 * .. Scalar Arguments ..
445  INTEGER info, lda, ldpt, ldq, ldx, lwork, nout, nrhs,
446  $ nsizes, ntypes
447  REAL thresh
448 * ..
449 * .. Array Arguments ..
450  LOGICAL dotype( * )
451  INTEGER iseed( 4 ), iwork( * ), mval( * ), nval( * )
452  REAL a( lda, * ), bd( * ), be( * ), pt( ldpt, * ),
453  $ q( ldq, * ), s1( * ), s2( * ), u( ldpt, * ),
454  $ vt( ldpt, * ), work( * ), x( ldx, * ),
455  $ y( ldx, * ), z( ldx, * )
456 * ..
457 *
458 * ======================================================================
459 *
460 * .. Parameters ..
461  REAL zero, one, two, half
462  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0,
463  $ half = 0.5e0 )
464  INTEGER maxtyp
465  parameter( maxtyp = 16 )
466 * ..
467 * .. Local Scalars ..
468  LOGICAL badmm, badnn, bidiag
469  CHARACTER uplo
470  CHARACTER*3 path
471  INTEGER i, iinfo, imode, itype, j, jcol, jsize, jtype,
472  $ log2ui, m, minwrk, mmax, mnmax, mnmin, mq,
473  $ mtypes, n, nfail, nmax, ntest
474  REAL amninv, anorm, cond, ovfl, rtovfl, rtunfl,
475  $ temp1, temp2, ulp, ulpinv, unfl
476 * ..
477 * .. Local Arrays ..
478  INTEGER idum( 1 ), ioldsd( 4 ), kmagn( maxtyp ),
479  $ kmode( maxtyp ), ktype( maxtyp )
480  REAL dum( 1 ), dumma( 1 ), result( 19 )
481 * ..
482 * .. External Functions ..
483  REAL slamch, slarnd
484  EXTERNAL slamch, slarnd
485 * ..
486 * .. External Subroutines ..
487  EXTERNAL alasum, sbdsdc, sbdsqr, sbdt01, sbdt02, sbdt03,
490 * ..
491 * .. Intrinsic Functions ..
492  INTRINSIC abs, exp, int, log, max, min, sqrt
493 * ..
494 * .. Scalars in Common ..
495  LOGICAL lerr, ok
496  CHARACTER*32 srnamt
497  INTEGER infot, nunit
498 * ..
499 * .. Common blocks ..
500  COMMON / infoc / infot, nunit, ok, lerr
501  COMMON / srnamc / srnamt
502 * ..
503 * .. Data statements ..
504  DATA ktype / 1, 2, 5*4, 5*6, 3*9, 10 /
505  DATA kmagn / 2*1, 3*1, 2, 3, 3*1, 2, 3, 1, 2, 3, 0 /
506  DATA kmode / 2*0, 4, 3, 1, 4, 4, 4, 3, 1, 4, 4, 0,
507  $ 0, 0, 0 /
508 * ..
509 * .. Executable Statements ..
510 *
511 * Check for errors
512 *
513  info = 0
514 *
515  badmm = .false.
516  badnn = .false.
517  mmax = 1
518  nmax = 1
519  mnmax = 1
520  minwrk = 1
521  DO 10 j = 1, nsizes
522  mmax = max( mmax, mval( j ) )
523  IF( mval( j ).LT.0 )
524  $ badmm = .true.
525  nmax = max( nmax, nval( j ) )
526  IF( nval( j ).LT.0 )
527  $ badnn = .true.
528  mnmax = max( mnmax, min( mval( j ), nval( j ) ) )
529  minwrk = max( minwrk, 3*( mval( j )+nval( j ) ),
530  $ mval( j )*( mval( j )+max( mval( j ), nval( j ),
531  $ nrhs )+1 )+nval( j )*min( nval( j ), mval( j ) ) )
532  10 CONTINUE
533 *
534 * Check for errors
535 *
536  IF( nsizes.LT.0 ) THEN
537  info = -1
538  ELSE IF( badmm ) THEN
539  info = -2
540  ELSE IF( badnn ) THEN
541  info = -3
542  ELSE IF( ntypes.LT.0 ) THEN
543  info = -4
544  ELSE IF( nrhs.LT.0 ) THEN
545  info = -6
546  ELSE IF( lda.LT.mmax ) THEN
547  info = -11
548  ELSE IF( ldx.LT.mmax ) THEN
549  info = -17
550  ELSE IF( ldq.LT.mmax ) THEN
551  info = -21
552  ELSE IF( ldpt.LT.mnmax ) THEN
553  info = -23
554  ELSE IF( minwrk.GT.lwork ) THEN
555  info = -27
556  END IF
557 *
558  IF( info.NE.0 ) THEN
559  CALL xerbla( 'SCHKBD', -info )
560  RETURN
561  END IF
562 *
563 * Initialize constants
564 *
565  path( 1: 1 ) = 'Single precision'
566  path( 2: 3 ) = 'BD'
567  nfail = 0
568  ntest = 0
569  unfl = slamch( 'Safe minimum' )
570  ovfl = slamch( 'Overflow' )
571  CALL slabad( unfl, ovfl )
572  ulp = slamch( 'Precision' )
573  ulpinv = one / ulp
574  log2ui = int( log( ulpinv ) / log( two ) )
575  rtunfl = sqrt( unfl )
576  rtovfl = sqrt( ovfl )
577  infot = 0
578 *
579 * Loop over sizes, types
580 *
581  DO 200 jsize = 1, nsizes
582  m = mval( jsize )
583  n = nval( jsize )
584  mnmin = min( m, n )
585  amninv = one / max( m, n, 1 )
586 *
587  IF( nsizes.NE.1 ) THEN
588  mtypes = min( maxtyp, ntypes )
589  ELSE
590  mtypes = min( maxtyp+1, ntypes )
591  END IF
592 *
593  DO 190 jtype = 1, mtypes
594  IF( .NOT.dotype( jtype ) )
595  $ GO TO 190
596 *
597  DO 20 j = 1, 4
598  ioldsd( j ) = iseed( j )
599  20 CONTINUE
600 *
601  DO 30 j = 1, 14
602  result( j ) = -one
603  30 CONTINUE
604 *
605  uplo = ' '
606 *
607 * Compute "A"
608 *
609 * Control parameters:
610 *
611 * KMAGN KMODE KTYPE
612 * =1 O(1) clustered 1 zero
613 * =2 large clustered 2 identity
614 * =3 small exponential (none)
615 * =4 arithmetic diagonal, (w/ eigenvalues)
616 * =5 random symmetric, w/ eigenvalues
617 * =6 nonsymmetric, w/ singular values
618 * =7 random diagonal
619 * =8 random symmetric
620 * =9 random nonsymmetric
621 * =10 random bidiagonal (log. distrib.)
622 *
623  IF( mtypes.GT.maxtyp )
624  $ GO TO 100
625 *
626  itype = ktype( jtype )
627  imode = kmode( jtype )
628 *
629 * Compute norm
630 *
631  GO TO ( 40, 50, 60 )kmagn( jtype )
632 *
633  40 CONTINUE
634  anorm = one
635  GO TO 70
636 *
637  50 CONTINUE
638  anorm = ( rtovfl*ulp )*amninv
639  GO TO 70
640 *
641  60 CONTINUE
642  anorm = rtunfl*max( m, n )*ulpinv
643  GO TO 70
644 *
645  70 CONTINUE
646 *
647  CALL slaset( 'Full', lda, n, zero, zero, a, lda )
648  iinfo = 0
649  cond = ulpinv
650 *
651  bidiag = .false.
652  IF( itype.EQ.1 ) THEN
653 *
654 * Zero matrix
655 *
656  iinfo = 0
657 *
658  ELSE IF( itype.EQ.2 ) THEN
659 *
660 * Identity
661 *
662  DO 80 jcol = 1, mnmin
663  a( jcol, jcol ) = anorm
664  80 CONTINUE
665 *
666  ELSE IF( itype.EQ.4 ) THEN
667 *
668 * Diagonal Matrix, [Eigen]values Specified
669 *
670  CALL slatms( mnmin, mnmin, 'S', iseed, 'N', work, imode,
671  $ cond, anorm, 0, 0, 'N', a, lda,
672  $ work( mnmin+1 ), iinfo )
673 *
674  ELSE IF( itype.EQ.5 ) THEN
675 *
676 * Symmetric, eigenvalues specified
677 *
678  CALL slatms( mnmin, mnmin, 'S', iseed, 'S', work, imode,
679  $ cond, anorm, m, n, 'N', a, lda,
680  $ work( mnmin+1 ), iinfo )
681 *
682  ELSE IF( itype.EQ.6 ) THEN
683 *
684 * Nonsymmetric, singular values specified
685 *
686  CALL slatms( m, n, 'S', iseed, 'N', work, imode, cond,
687  $ anorm, m, n, 'N', a, lda, work( mnmin+1 ),
688  $ iinfo )
689 *
690  ELSE IF( itype.EQ.7 ) THEN
691 *
692 * Diagonal, random entries
693 *
694  CALL slatmr( mnmin, mnmin, 'S', iseed, 'N', work, 6, one,
695  $ one, 'T', 'N', work( mnmin+1 ), 1, one,
696  $ work( 2*mnmin+1 ), 1, one, 'N', iwork, 0, 0,
697  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
698 *
699  ELSE IF( itype.EQ.8 ) THEN
700 *
701 * Symmetric, random entries
702 *
703  CALL slatmr( mnmin, mnmin, 'S', iseed, 'S', work, 6, one,
704  $ one, 'T', 'N', work( mnmin+1 ), 1, one,
705  $ work( m+mnmin+1 ), 1, one, 'N', iwork, m, n,
706  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
707 *
708  ELSE IF( itype.EQ.9 ) THEN
709 *
710 * Nonsymmetric, random entries
711 *
712  CALL slatmr( m, n, 'S', iseed, 'N', work, 6, one, one,
713  $ 'T', 'N', work( mnmin+1 ), 1, one,
714  $ work( m+mnmin+1 ), 1, one, 'N', iwork, m, n,
715  $ zero, anorm, 'NO', a, lda, iwork, iinfo )
716 *
717  ELSE IF( itype.EQ.10 ) THEN
718 *
719 * Bidiagonal, random entries
720 *
721  temp1 = -two*log( ulp )
722  DO 90 j = 1, mnmin
723  bd( j ) = exp( temp1*slarnd( 2, iseed ) )
724  IF( j.LT.mnmin )
725  $ be( j ) = exp( temp1*slarnd( 2, iseed ) )
726  90 CONTINUE
727 *
728  iinfo = 0
729  bidiag = .true.
730  IF( m.GE.n ) THEN
731  uplo = 'U'
732  ELSE
733  uplo = 'L'
734  END IF
735  ELSE
736  iinfo = 1
737  END IF
738 *
739  IF( iinfo.EQ.0 ) THEN
740 *
741 * Generate Right-Hand Side
742 *
743  IF( bidiag ) THEN
744  CALL slatmr( mnmin, nrhs, 'S', iseed, 'N', work, 6,
745  $ one, one, 'T', 'N', work( mnmin+1 ), 1,
746  $ one, work( 2*mnmin+1 ), 1, one, 'N',
747  $ iwork, mnmin, nrhs, zero, one, 'NO', y,
748  $ ldx, iwork, iinfo )
749  ELSE
750  CALL slatmr( m, nrhs, 'S', iseed, 'N', work, 6, one,
751  $ one, 'T', 'N', work( m+1 ), 1, one,
752  $ work( 2*m+1 ), 1, one, 'N', iwork, m,
753  $ nrhs, zero, one, 'NO', x, ldx, iwork,
754  $ iinfo )
755  END IF
756  END IF
757 *
758 * Error Exit
759 *
760  IF( iinfo.NE.0 ) THEN
761  WRITE( nout, fmt = 9998 )'Generator', iinfo, m, n,
762  $ jtype, ioldsd
763  info = abs( iinfo )
764  RETURN
765  END IF
766 *
767  100 CONTINUE
768 *
769 * Call SGEBRD and SORGBR to compute B, Q, and P, do tests.
770 *
771  IF( .NOT.bidiag ) THEN
772 *
773 * Compute transformations to reduce A to bidiagonal form:
774 * B := Q' * A * P.
775 *
776  CALL slacpy( ' ', m, n, a, lda, q, ldq )
777  CALL sgebrd( m, n, q, ldq, bd, be, work, work( mnmin+1 ),
778  $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
779 *
780 * Check error code from SGEBRD.
781 *
782  IF( iinfo.NE.0 ) THEN
783  WRITE( nout, fmt = 9998 )'SGEBRD', iinfo, m, n,
784  $ jtype, ioldsd
785  info = abs( iinfo )
786  RETURN
787  END IF
788 *
789  CALL slacpy( ' ', m, n, q, ldq, pt, ldpt )
790  IF( m.GE.n ) THEN
791  uplo = 'U'
792  ELSE
793  uplo = 'L'
794  END IF
795 *
796 * Generate Q
797 *
798  mq = m
799  IF( nrhs.LE.0 )
800  $ mq = mnmin
801  CALL sorgbr( 'Q', m, mq, n, q, ldq, work,
802  $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
803 *
804 * Check error code from SORGBR.
805 *
806  IF( iinfo.NE.0 ) THEN
807  WRITE( nout, fmt = 9998 )'SORGBR(Q)', iinfo, m, n,
808  $ jtype, ioldsd
809  info = abs( iinfo )
810  RETURN
811  END IF
812 *
813 * Generate P'
814 *
815  CALL sorgbr( 'P', mnmin, n, m, pt, ldpt, work( mnmin+1 ),
816  $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
817 *
818 * Check error code from SORGBR.
819 *
820  IF( iinfo.NE.0 ) THEN
821  WRITE( nout, fmt = 9998 )'SORGBR(P)', iinfo, m, n,
822  $ jtype, ioldsd
823  info = abs( iinfo )
824  RETURN
825  END IF
826 *
827 * Apply Q' to an M by NRHS matrix X: Y := Q' * X.
828 *
829  CALL sgemm( 'Transpose', 'No transpose', m, nrhs, m, one,
830  $ q, ldq, x, ldx, zero, y, ldx )
831 *
832 * Test 1: Check the decomposition A := Q * B * PT
833 * 2: Check the orthogonality of Q
834 * 3: Check the orthogonality of PT
835 *
836  CALL sbdt01( m, n, 1, a, lda, q, ldq, bd, be, pt, ldpt,
837  $ work, result( 1 ) )
838  CALL sort01( 'Columns', m, mq, q, ldq, work, lwork,
839  $ result( 2 ) )
840  CALL sort01( 'Rows', mnmin, n, pt, ldpt, work, lwork,
841  $ result( 3 ) )
842  END IF
843 *
844 * Use SBDSQR to form the SVD of the bidiagonal matrix B:
845 * B := U * S1 * VT, and compute Z = U' * Y.
846 *
847  CALL scopy( mnmin, bd, 1, s1, 1 )
848  IF( mnmin.GT.0 )
849  $ CALL scopy( mnmin-1, be, 1, work, 1 )
850  CALL slacpy( ' ', m, nrhs, y, ldx, z, ldx )
851  CALL slaset( 'Full', mnmin, mnmin, zero, one, u, ldpt )
852  CALL slaset( 'Full', mnmin, mnmin, zero, one, vt, ldpt )
853 *
854  CALL sbdsqr( uplo, mnmin, mnmin, mnmin, nrhs, s1, work, vt,
855  $ ldpt, u, ldpt, z, ldx, work( mnmin+1 ), iinfo )
856 *
857 * Check error code from SBDSQR.
858 *
859  IF( iinfo.NE.0 ) THEN
860  WRITE( nout, fmt = 9998 )'SBDSQR(vects)', iinfo, m, n,
861  $ jtype, ioldsd
862  info = abs( iinfo )
863  IF( iinfo.LT.0 ) THEN
864  RETURN
865  ELSE
866  result( 4 ) = ulpinv
867  GO TO 170
868  END IF
869  END IF
870 *
871 * Use SBDSQR to compute only the singular values of the
872 * bidiagonal matrix B; U, VT, and Z should not be modified.
873 *
874  CALL scopy( mnmin, bd, 1, s2, 1 )
875  IF( mnmin.GT.0 )
876  $ CALL scopy( mnmin-1, be, 1, work, 1 )
877 *
878  CALL sbdsqr( uplo, mnmin, 0, 0, 0, s2, work, vt, ldpt, u,
879  $ ldpt, z, ldx, work( mnmin+1 ), iinfo )
880 *
881 * Check error code from SBDSQR.
882 *
883  IF( iinfo.NE.0 ) THEN
884  WRITE( nout, fmt = 9998 )'SBDSQR(values)', iinfo, m, n,
885  $ jtype, ioldsd
886  info = abs( iinfo )
887  IF( iinfo.LT.0 ) THEN
888  RETURN
889  ELSE
890  result( 9 ) = ulpinv
891  GO TO 170
892  END IF
893  END IF
894 *
895 * Test 4: Check the decomposition B := U * S1 * VT
896 * 5: Check the computation Z := U' * Y
897 * 6: Check the orthogonality of U
898 * 7: Check the orthogonality of VT
899 *
900  CALL sbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
901  $ work, result( 4 ) )
902  CALL sbdt02( mnmin, nrhs, y, ldx, z, ldx, u, ldpt, work,
903  $ result( 5 ) )
904  CALL sort01( 'Columns', mnmin, mnmin, u, ldpt, work, lwork,
905  $ result( 6 ) )
906  CALL sort01( 'Rows', mnmin, mnmin, vt, ldpt, work, lwork,
907  $ result( 7 ) )
908 *
909 * Test 8: Check that the singular values are sorted in
910 * non-increasing order and are non-negative
911 *
912  result( 8 ) = zero
913  DO 110 i = 1, mnmin - 1
914  IF( s1( i ).LT.s1( i+1 ) )
915  $ result( 8 ) = ulpinv
916  IF( s1( i ).LT.zero )
917  $ result( 8 ) = ulpinv
918  110 CONTINUE
919  IF( mnmin.GE.1 ) THEN
920  IF( s1( mnmin ).LT.zero )
921  $ result( 8 ) = ulpinv
922  END IF
923 *
924 * Test 9: Compare SBDSQR with and without singular vectors
925 *
926  temp2 = zero
927 *
928  DO 120 j = 1, mnmin
929  temp1 = abs( s1( j )-s2( j ) ) /
930  $ max( sqrt( unfl )*max( s1( 1 ), one ),
931  $ ulp*max( abs( s1( j ) ), abs( s2( j ) ) ) )
932  temp2 = max( temp1, temp2 )
933  120 CONTINUE
934 *
935  result( 9 ) = temp2
936 *
937 * Test 10: Sturm sequence test of singular values
938 * Go up by factors of two until it succeeds
939 *
940  temp1 = thresh*( half-ulp )
941 *
942  DO 130 j = 0, log2ui
943 * CALL SSVDCH( MNMIN, BD, BE, S1, TEMP1, IINFO )
944  IF( iinfo.EQ.0 )
945  $ GO TO 140
946  temp1 = temp1*two
947  130 CONTINUE
948 *
949  140 CONTINUE
950  result( 10 ) = temp1
951 *
952 * Use SBDSQR to form the decomposition A := (QU) S (VT PT)
953 * from the bidiagonal form A := Q B PT.
954 *
955  IF( .NOT.bidiag ) THEN
956  CALL scopy( mnmin, bd, 1, s2, 1 )
957  IF( mnmin.GT.0 )
958  $ CALL scopy( mnmin-1, be, 1, work, 1 )
959 *
960  CALL sbdsqr( uplo, mnmin, n, m, nrhs, s2, work, pt, ldpt,
961  $ q, ldq, y, ldx, work( mnmin+1 ), iinfo )
962 *
963 * Test 11: Check the decomposition A := Q*U * S2 * VT*PT
964 * 12: Check the computation Z := U' * Q' * X
965 * 13: Check the orthogonality of Q*U
966 * 14: Check the orthogonality of VT*PT
967 *
968  CALL sbdt01( m, n, 0, a, lda, q, ldq, s2, dumma, pt,
969  $ ldpt, work, result( 11 ) )
970  CALL sbdt02( m, nrhs, x, ldx, y, ldx, q, ldq, work,
971  $ result( 12 ) )
972  CALL sort01( 'Columns', m, mq, q, ldq, work, lwork,
973  $ result( 13 ) )
974  CALL sort01( 'Rows', mnmin, n, pt, ldpt, work, lwork,
975  $ result( 14 ) )
976  END IF
977 *
978 * Use SBDSDC to form the SVD of the bidiagonal matrix B:
979 * B := U * S1 * VT
980 *
981  CALL scopy( mnmin, bd, 1, s1, 1 )
982  IF( mnmin.GT.0 )
983  $ CALL scopy( mnmin-1, be, 1, work, 1 )
984  CALL slaset( 'Full', mnmin, mnmin, zero, one, u, ldpt )
985  CALL slaset( 'Full', mnmin, mnmin, zero, one, vt, ldpt )
986 *
987  CALL sbdsdc( uplo, 'I', mnmin, s1, work, u, ldpt, vt, ldpt,
988  $ dum, idum, work( mnmin+1 ), iwork, iinfo )
989 *
990 * Check error code from SBDSDC.
991 *
992  IF( iinfo.NE.0 ) THEN
993  WRITE( nout, fmt = 9998 )'SBDSDC(vects)', iinfo, m, n,
994  $ jtype, ioldsd
995  info = abs( iinfo )
996  IF( iinfo.LT.0 ) THEN
997  RETURN
998  ELSE
999  result( 15 ) = ulpinv
1000  GO TO 170
1001  END IF
1002  END IF
1003 *
1004 * Use SBDSDC to compute only the singular values of the
1005 * bidiagonal matrix B; U and VT should not be modified.
1006 *
1007  CALL scopy( mnmin, bd, 1, s2, 1 )
1008  IF( mnmin.GT.0 )
1009  $ CALL scopy( mnmin-1, be, 1, work, 1 )
1010 *
1011  CALL sbdsdc( uplo, 'N', mnmin, s2, work, dum, 1, dum, 1,
1012  $ dum, idum, work( mnmin+1 ), iwork, iinfo )
1013 *
1014 * Check error code from SBDSDC.
1015 *
1016  IF( iinfo.NE.0 ) THEN
1017  WRITE( nout, fmt = 9998 )'SBDSDC(values)', iinfo, m, n,
1018  $ jtype, ioldsd
1019  info = abs( iinfo )
1020  IF( iinfo.LT.0 ) THEN
1021  RETURN
1022  ELSE
1023  result( 18 ) = ulpinv
1024  GO TO 170
1025  END IF
1026  END IF
1027 *
1028 * Test 15: Check the decomposition B := U * S1 * VT
1029 * 16: Check the orthogonality of U
1030 * 17: Check the orthogonality of VT
1031 *
1032  CALL sbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
1033  $ work, result( 15 ) )
1034  CALL sort01( 'Columns', mnmin, mnmin, u, ldpt, work, lwork,
1035  $ result( 16 ) )
1036  CALL sort01( 'Rows', mnmin, mnmin, vt, ldpt, work, lwork,
1037  $ result( 17 ) )
1038 *
1039 * Test 18: Check that the singular values are sorted in
1040 * non-increasing order and are non-negative
1041 *
1042  result( 18 ) = zero
1043  DO 150 i = 1, mnmin - 1
1044  IF( s1( i ).LT.s1( i+1 ) )
1045  $ result( 18 ) = ulpinv
1046  IF( s1( i ).LT.zero )
1047  $ result( 18 ) = ulpinv
1048  150 CONTINUE
1049  IF( mnmin.GE.1 ) THEN
1050  IF( s1( mnmin ).LT.zero )
1051  $ result( 18 ) = ulpinv
1052  END IF
1053 *
1054 * Test 19: Compare SBDSQR with and without singular vectors
1055 *
1056  temp2 = zero
1057 *
1058  DO 160 j = 1, mnmin
1059  temp1 = abs( s1( j )-s2( j ) ) /
1060  $ max( sqrt( unfl )*max( s1( 1 ), one ),
1061  $ ulp*max( abs( s1( 1 ) ), abs( s2( 1 ) ) ) )
1062  temp2 = max( temp1, temp2 )
1063  160 CONTINUE
1064 *
1065  result( 19 ) = temp2
1066 *
1067 * End of Loop -- Check for RESULT(j) > THRESH
1068 *
1069  170 CONTINUE
1070  DO 180 j = 1, 19
1071  IF( result( j ).GE.thresh ) THEN
1072  IF( nfail.EQ.0 )
1073  $ CALL slahd2( nout, path )
1074  WRITE( nout, fmt = 9999 )m, n, jtype, ioldsd, j,
1075  $ result( j )
1076  nfail = nfail + 1
1077  END IF
1078  180 CONTINUE
1079  IF( .NOT.bidiag ) THEN
1080  ntest = ntest + 19
1081  ELSE
1082  ntest = ntest + 5
1083  END IF
1084 *
1085  190 CONTINUE
1086  200 CONTINUE
1087 *
1088 * Summary
1089 *
1090  CALL alasum( path, nout, nfail, ntest, 0 )
1091 *
1092  RETURN
1093 *
1094 * End of SCHKBD
1095 *
1096  9999 FORMAT( ' M=', i5, ', N=', i5, ', type ', i2, ', seed=',
1097  $ 4( i4, ',' ), ' test(', i2, ')=', g11.4 )
1098  9998 FORMAT( ' SCHKBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1099  $ i6, ', N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
1100  $ i5, ')' )
1101 *
subroutine sgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
SGEBRD
Definition: sgebrd.f:207
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:112
subroutine sbdt02(M, N, B, LDB, C, LDC, U, LDU, WORK, RESID)
SBDT02
Definition: sbdt02.f:113
subroutine sbdt03(UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK, RESID)
SBDT03
Definition: sbdt03.f:137
subroutine slahd2(IOUNIT, PATH)
SLAHD2
Definition: slahd2.f:67
subroutine slatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
SLATMS
Definition: slatms.f:323
subroutine slatmr(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, RSIGN, GRADE, DL, MODEL, CONDL, DR, MODER, CONDR, PIVTNG, IPIVOT, KL, KU, SPARSE, ANORM, PACK, A, LDA, IWORK, INFO)
SLATMR
Definition: slatmr.f:473
subroutine sorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGBR
Definition: sorgbr.f:159
subroutine sbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
SBDSQR
Definition: sbdsqr.f:232
subroutine sbdsdc(UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
SBDSDC
Definition: sbdsdc.f:207
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
subroutine alasum(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASUM
Definition: alasum.f:75
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine sort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
SORT01
Definition: sort01.f:118
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
real function slarnd(IDIST, ISEED)
SLARND
Definition: slarnd.f:75
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine sbdt01(M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK, RESID)
SBDT01
Definition: sbdt01.f:142

Here is the call graph for this function:

Here is the caller graph for this function: