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

Go to the source code of this file.

Functions/Subroutines

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

Function/Subroutine Documentation

subroutine dchkbd ( integer  NSIZES,
integer, dimension( * )  MVAL,
integer, dimension( * )  NVAL,
integer  NTYPES,
logical, dimension( * )  DOTYPE,
integer  NRHS,
integer, dimension( 4 )  ISEED,
double precision  THRESH,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  BD,
double precision, dimension( * )  BE,
double precision, dimension( * )  S1,
double precision, dimension( * )  S2,
double precision, dimension( ldx, * )  X,
integer  LDX,
double precision, dimension( ldx, * )  Y,
double precision, dimension( ldx, * )  Z,
double precision, dimension( ldq, * )  Q,
integer  LDQ,
double precision, dimension( ldpt, * )  PT,
integer  LDPT,
double precision, dimension( ldpt, * )  U,
double precision, dimension( ldpt, * )  VT,
double precision, dimension( * )  WORK,
integer  LWORK,
integer, dimension( * )  IWORK,
integer  NOUT,
integer  INFO 
)

DCHKBD

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

 DGEBRD 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.

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

 DBDSQR 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, DBDSQR has an option to apply the left orthogonal matrix
 U to a matrix X, useful in least squares applications.

 DBDSDC 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 DGEBRD and DORGBR

 (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 DBDSQR 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
       DSVDCH)

 Test DBDSQR 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 DBDSDC 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) DGEBRD 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, DCHKBD
          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 DBDSQR.  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 DCHKBD to continue the same random
          number sequence.
[in]THRESH
          THRESH is DOUBLE PRECISION
          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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension
                      (max(min(MVAL(j),NVAL(j))))
[out]BE
          BE is DOUBLE PRECISION array, dimension
                      (max(min(MVAL(j),NVAL(j))))
[out]S1
          S1 is DOUBLE PRECISION array, dimension
                      (max(min(MVAL(j),NVAL(j))))
[out]S2
          S2 is DOUBLE PRECISION array, dimension
                      (max(min(MVAL(j),NVAL(j))))
[out]X
          X is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (LDX,NRHS)
[out]Z
          Z is DOUBLE PRECISION array, dimension (LDX,NRHS)
[out]Q
          Q is DOUBLE PRECISION array, dimension (LDQ,MMAX)
[in]LDQ
          LDQ is INTEGER
          The leading dimension of the array Q.  LDQ >= max(1,MMAX).
[out]PT
          PT is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension
                      (LDPT,max(min(MVAL(j),NVAL(j))))
[out]VT
          VT is DOUBLE PRECISION array, dimension
                      (LDPT,max(min(MVAL(j),NVAL(j))))
[out]WORK
          WORK is DOUBLE PRECISION 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  DLATMR, SLATMS, DGEBRD, DORGBR, or DBDSQR,
              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 dchkbd.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  DOUBLE PRECISION thresh
448 * ..
449 * .. Array Arguments ..
450  LOGICAL dotype( * )
451  INTEGER iseed( 4 ), iwork( * ), mval( * ), nval( * )
452  DOUBLE PRECISION 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  DOUBLE PRECISION zero, one, two, half
462  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0,
463  $ half = 0.5d0 )
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  DOUBLE PRECISION 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  DOUBLE PRECISION dum( 1 ), dumma( 1 ), result( 19 )
481 * ..
482 * .. External Functions ..
483  DOUBLE PRECISION dlamch, dlarnd
484  EXTERNAL dlamch, dlarnd
485 * ..
486 * .. External Subroutines ..
487  EXTERNAL alasum, dbdsdc, dbdsqr, dbdt01, dbdt02, dbdt03,
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( 'DCHKBD', -info )
560  RETURN
561  END IF
562 *
563 * Initialize constants
564 *
565  path( 1: 1 ) = 'Double precision'
566  path( 2: 3 ) = 'BD'
567  nfail = 0
568  ntest = 0
569  unfl = dlamch( 'Safe minimum' )
570  ovfl = dlamch( 'Overflow' )
571  CALL dlabad( unfl, ovfl )
572  ulp = dlamch( '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 dlaset( '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 dlatms( 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 dlatms( 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 dlatms( 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 dlatmr( 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 dlatmr( 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 dlatmr( 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*dlarnd( 2, iseed ) )
724  IF( j.LT.mnmin )
725  $ be( j ) = exp( temp1*dlarnd( 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 dlatmr( 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 dlatmr( 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 DGEBRD and DORGBR 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 dlacpy( ' ', m, n, a, lda, q, ldq )
777  CALL dgebrd( 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 DGEBRD.
781 *
782  IF( iinfo.NE.0 ) THEN
783  WRITE( nout, fmt = 9998 )'DGEBRD', iinfo, m, n,
784  $ jtype, ioldsd
785  info = abs( iinfo )
786  RETURN
787  END IF
788 *
789  CALL dlacpy( ' ', 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 dorgbr( 'Q', m, mq, n, q, ldq, work,
802  $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
803 *
804 * Check error code from DORGBR.
805 *
806  IF( iinfo.NE.0 ) THEN
807  WRITE( nout, fmt = 9998 )'DORGBR(Q)', iinfo, m, n,
808  $ jtype, ioldsd
809  info = abs( iinfo )
810  RETURN
811  END IF
812 *
813 * Generate P'
814 *
815  CALL dorgbr( 'P', mnmin, n, m, pt, ldpt, work( mnmin+1 ),
816  $ work( 2*mnmin+1 ), lwork-2*mnmin, iinfo )
817 *
818 * Check error code from DORGBR.
819 *
820  IF( iinfo.NE.0 ) THEN
821  WRITE( nout, fmt = 9998 )'DORGBR(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 dgemm( '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 dbdt01( m, n, 1, a, lda, q, ldq, bd, be, pt, ldpt,
837  $ work, result( 1 ) )
838  CALL dort01( 'Columns', m, mq, q, ldq, work, lwork,
839  $ result( 2 ) )
840  CALL dort01( 'Rows', mnmin, n, pt, ldpt, work, lwork,
841  $ result( 3 ) )
842  END IF
843 *
844 * Use DBDSQR to form the SVD of the bidiagonal matrix B:
845 * B := U * S1 * VT, and compute Z = U' * Y.
846 *
847  CALL dcopy( mnmin, bd, 1, s1, 1 )
848  IF( mnmin.GT.0 )
849  $ CALL dcopy( mnmin-1, be, 1, work, 1 )
850  CALL dlacpy( ' ', m, nrhs, y, ldx, z, ldx )
851  CALL dlaset( 'Full', mnmin, mnmin, zero, one, u, ldpt )
852  CALL dlaset( 'Full', mnmin, mnmin, zero, one, vt, ldpt )
853 *
854  CALL dbdsqr( uplo, mnmin, mnmin, mnmin, nrhs, s1, work, vt,
855  $ ldpt, u, ldpt, z, ldx, work( mnmin+1 ), iinfo )
856 *
857 * Check error code from DBDSQR.
858 *
859  IF( iinfo.NE.0 ) THEN
860  WRITE( nout, fmt = 9998 )'DBDSQR(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 DBDSQR to compute only the singular values of the
872 * bidiagonal matrix B; U, VT, and Z should not be modified.
873 *
874  CALL dcopy( mnmin, bd, 1, s2, 1 )
875  IF( mnmin.GT.0 )
876  $ CALL dcopy( mnmin-1, be, 1, work, 1 )
877 *
878  CALL dbdsqr( uplo, mnmin, 0, 0, 0, s2, work, vt, ldpt, u,
879  $ ldpt, z, ldx, work( mnmin+1 ), iinfo )
880 *
881 * Check error code from DBDSQR.
882 *
883  IF( iinfo.NE.0 ) THEN
884  WRITE( nout, fmt = 9998 )'DBDSQR(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 dbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
901  $ work, result( 4 ) )
902  CALL dbdt02( mnmin, nrhs, y, ldx, z, ldx, u, ldpt, work,
903  $ result( 5 ) )
904  CALL dort01( 'Columns', mnmin, mnmin, u, ldpt, work, lwork,
905  $ result( 6 ) )
906  CALL dort01( '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 DBDSQR 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 DSVDCH( 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 DBDSQR 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 dcopy( mnmin, bd, 1, s2, 1 )
957  IF( mnmin.GT.0 )
958  $ CALL dcopy( mnmin-1, be, 1, work, 1 )
959 *
960  CALL dbdsqr( 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 dbdt01( m, n, 0, a, lda, q, ldq, s2, dumma, pt,
969  $ ldpt, work, result( 11 ) )
970  CALL dbdt02( m, nrhs, x, ldx, y, ldx, q, ldq, work,
971  $ result( 12 ) )
972  CALL dort01( 'Columns', m, mq, q, ldq, work, lwork,
973  $ result( 13 ) )
974  CALL dort01( 'Rows', mnmin, n, pt, ldpt, work, lwork,
975  $ result( 14 ) )
976  END IF
977 *
978 * Use DBDSDC to form the SVD of the bidiagonal matrix B:
979 * B := U * S1 * VT
980 *
981  CALL dcopy( mnmin, bd, 1, s1, 1 )
982  IF( mnmin.GT.0 )
983  $ CALL dcopy( mnmin-1, be, 1, work, 1 )
984  CALL dlaset( 'Full', mnmin, mnmin, zero, one, u, ldpt )
985  CALL dlaset( 'Full', mnmin, mnmin, zero, one, vt, ldpt )
986 *
987  CALL dbdsdc( uplo, 'I', mnmin, s1, work, u, ldpt, vt, ldpt,
988  $ dum, idum, work( mnmin+1 ), iwork, iinfo )
989 *
990 * Check error code from DBDSDC.
991 *
992  IF( iinfo.NE.0 ) THEN
993  WRITE( nout, fmt = 9998 )'DBDSDC(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 DBDSDC to compute only the singular values of the
1005 * bidiagonal matrix B; U and VT should not be modified.
1006 *
1007  CALL dcopy( mnmin, bd, 1, s2, 1 )
1008  IF( mnmin.GT.0 )
1009  $ CALL dcopy( mnmin-1, be, 1, work, 1 )
1010 *
1011  CALL dbdsdc( uplo, 'N', mnmin, s2, work, dum, 1, dum, 1,
1012  $ dum, idum, work( mnmin+1 ), iwork, iinfo )
1013 *
1014 * Check error code from DBDSDC.
1015 *
1016  IF( iinfo.NE.0 ) THEN
1017  WRITE( nout, fmt = 9998 )'DBDSDC(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 dbdt03( uplo, mnmin, 1, bd, be, u, ldpt, s1, vt, ldpt,
1033  $ work, result( 15 ) )
1034  CALL dort01( 'Columns', mnmin, mnmin, u, ldpt, work, lwork,
1035  $ result( 16 ) )
1036  CALL dort01( '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 DBDSQR 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 dlahd2( 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 DCHKBD
1095 *
1096  9999 FORMAT( ' M=', i5, ', N=', i5, ', type ', i2, ', seed=',
1097  $ 4( i4, ',' ), ' test(', i2, ')=', g11.4 )
1098  9998 FORMAT( ' DCHKBD: ', a, ' returned INFO=', i6, '.', / 9x, 'M=',
1099  $ i6, ', N=', i6, ', JTYPE=', i6, ', ISEED=(', 3( i5, ',' ),
1100  $ i5, ')' )
1101 *
subroutine dbdsdc(UPLO, COMPQ, N, D, E, U, LDU, VT, LDVT, Q, IQ, WORK, IWORK, INFO)
DBDSDC
Definition: dbdsdc.f:207
subroutine dlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
DLATMS
Definition: dlatms.f:323
double precision function dlarnd(IDIST, ISEED)
DLARND
Definition: dlarnd.f:75
subroutine dbdt01(M, N, KD, A, LDA, Q, LDQ, D, E, PT, LDPT, WORK, RESID)
DBDT01
Definition: dbdt01.f:142
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine dlatmr(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, RSIGN, GRADE, DL, MODEL, CONDL, DR, MODER, CONDR, PIVTNG, IPIVOT, KL, KU, SPARSE, ANORM, PACK, A, LDA, IWORK, INFO)
DLATMR
Definition: dlatmr.f:473
subroutine dgebrd(M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK, INFO)
DGEBRD
Definition: dgebrd.f:207
subroutine alasum(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASUM
Definition: alasum.f:75
subroutine dbdt02(M, N, B, LDB, C, LDC, U, LDU, WORK, RESID)
DBDT02
Definition: dbdt02.f:113
subroutine dorgbr(VECT, M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGBR
Definition: dorgbr.f:159
subroutine dlahd2(IOUNIT, PATH)
DLAHD2
Definition: dlahd2.f:67
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:112
subroutine dbdt03(UPLO, N, KD, D, E, U, LDU, S, VT, LDVT, WORK, RESID)
DBDT03
Definition: dbdt03.f:137
subroutine dbdsqr(UPLO, N, NCVT, NRU, NCC, D, E, VT, LDVT, U, LDU, C, LDC, WORK, INFO)
DBDSQR
Definition: dbdsqr.f:232
subroutine dort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
DORT01
Definition: dort01.f:118

Here is the call graph for this function:

Here is the caller graph for this function: