sdrvbd.c

Go to the documentation of this file.

/* sdrvbd.f -- translated by f2c (version 20061008).
   You must link the resulting object file with libf2c:
        on Microsoft Windows system, link with libf2c.lib;
        on Linux or Unix systems, link with .../path/to/libf2c.a -lm
        or, if you install libf2c.a in a standard place, with -lf2c -lm
        -- in that order, at the end of the command line, as in
                cc *.o -lf2c -lm
        Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

                http://www.netlib.org/f2c/libf2c.zip
*/

#include "f2c.h"
#include "blaswrap.h"

/* Common Block Declarations */

struct {
    integer infot, nunit;
    logical ok, lerr;
} infoc_;

#define infoc_1 infoc_

struct {
    char srnamt[32];
} srnamc_;

#define srnamc_1 srnamc_

/* Table of constant values */

static real c_b13 = 0.f;
static real c_b17 = 1.f;
static integer c__4 = 4;
static integer c__1 = 1;
static integer c__0 = 0;

/* Subroutine */ int sdrvbd_(integer *nsizes, integer *mm, integer *nn, 
        integer *ntypes, logical *dotype, integer *iseed, real *thresh, real *
        a, integer *lda, real *u, integer *ldu, real *vt, integer *ldvt, real 
        *asav, real *usav, real *vtsav, real *s, real *ssav, real *e, real *
        work, integer *lwork, integer *iwork, integer *nout, integer *info)
{
    /* Initialized data */

    static char cjob[1*4] = "N" "O" "S" "A";

    /* Format strings */
    static char fmt_9996[] = "(\002 SDRVBD: \002,a,\002 returned INFO=\002,i"
            "6,\002.\002,/9x,\002M=\002,i6,\002, N=\002,i6,\002, JTYPE=\002,i"
            "6,\002, ISEED=(\002,3(i5,\002,\002),i5,\002)\002)";
    static char fmt_9995[] = "(\002 SDRVBD: \002,a,\002 returned INFO=\002,i"
            "6,\002.\002,/9x,\002M=\002,i6,\002, N=\002,i6,\002, JTYPE=\002,i"
            "6,\002, LSWORK=\002,i6,/9x,\002ISEED=(\002,3(i5,\002,\002),i5"
            ",\002)\002)";
    static char fmt_9999[] = "(\002 SVD -- Real Singular Value Decomposition"
            " Driver \002,/\002 Matrix types (see SDRVBD for details):\002,/"
            "/\002 1 = Zero matrix\002,/\002 2 = Identity matrix\002,/\002 3 "
            "= Evenly spaced singular values near 1\002,/\002 4 = Evenly spac"
            "ed singular values near underflow\002,/\002 5 = Evenly spaced si"
            "ngular values near overflow\002,//\002 Tests performed: ( A is d"
            "ense, U and V are orthogonal,\002,/19x,\002 S is an array, and U"
            "partial, VTpartial, and\002,/19x,\002 Spartial are partially com"
            "puted U, VT and S),\002,/)";
    static char fmt_9998[] = "(\002 1 = | A - U diag(S) VT | / ( |A| max(M,N"
            ") ulp ) \002,/\002 2 = | I - U**T U | / ( M ulp ) \002,/\002 3 ="
            " | I - VT VT**T | / ( N ulp ) \002,/\002 4 = 0 if S contains min"
            "(M,N) nonnegative values in\002,\002 decreasing order, else 1/ulp"
            "\002,/\002 5 = | U - Upartial | / ( M ulp )\002,/\002 6 = | VT -"
            " VTpartial | / ( N ulp )\002,/\002 7 = | S - Spartial | / ( min("
            "M,N) ulp |S| )\002,/\002 8 = | A - U diag(S) VT | / ( |A| max(M,"
            "N) ulp ) \002,/\002 9 = | I - U**T U | / ( M ulp ) \002,/\00210 "
            "= | I - VT VT**T | / ( N ulp ) \002,/\00211 = 0 if S contains mi"
            "n(M,N) nonnegative values in\002,\002 decreasing order, else 1/u"
            "lp\002,/\00212 = | U - Upartial | / ( M ulp )\002,/\00213 = | VT"
            " - VTpartial | / ( N ulp )\002,/\00214 = | S - Spartial | / ( mi"
            "n(M,N) ulp |S| )\002,/\00215 = | A - U diag(S) VT | / ( |A| max("
            "M,N) ulp ) \002,/\00216 = | I - U**T U | / ( M ulp ) \002,/\0021"
            "7 = | I - VT VT**T | / ( N ulp ) \002,/\00218 = 0 if S contains "
            "min(M,N) nonnegative values in\002,\002 decreasing order, else 1"
            "/ulp\002,/\00219 = | U - Upartial | / ( M ulp )\002,/\00220 = | "
            "VT - VTpartial | / ( N ulp )\002,/\00221 = | S - Spartial | / ( "
            "min(M,N) ulp |S| )\002,//)";
    static char fmt_9997[] = "(\002 M=\002,i5,\002, N=\002,i5,\002, type "
            "\002,i1,\002, IWS=\002,i1,\002, seed=\002,4(i4,\002,\002),\002 t"
            "est(\002,i2,\002)=\002,g11.4)";

    /* System generated locals */
    integer a_dim1, a_offset, asav_dim1, asav_offset, u_dim1, u_offset, 
            usav_dim1, usav_offset, vt_dim1, vt_offset, vtsav_dim1, 
            vtsav_offset, i__1, i__2, i__3, i__4, i__5, i__6, i__7, i__8, 
            i__9, i__10, i__11, i__12, i__13, i__14;
    real r__1, r__2, r__3;

    /* Builtin functions */
    /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
    integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);

    /* Local variables */
    integer i__, j, m, n;
    real dif, div;
    integer ijq, iju;
    real ulp;
    integer iws;
    char jobq[1], path[3], jobu[1];
    integer mmax, nmax;
    real unfl, ovfl;
    integer ijvt;
    logical badmm, badnn;
    integer nfail;
    extern /* Subroutine */ int sbdt01_(integer *, integer *, integer *, real 
            *, integer *, real *, integer *, real *, real *, real *, integer *
, real *, real *);
    integer iinfo;
    real anorm;
    integer mnmin, mnmax;
    char jobvt[1];
    integer jsize;
    extern /* Subroutine */ int sort01_(char *, integer *, integer *, real *, 
            integer *, real *, integer *, real *), sort03_(char *, 
            integer *, integer *, integer *, integer *, real *, integer *, 
            real *, integer *, real *, integer *, real *, integer *);
    integer jtype, ntest, iwtmp;
    extern /* Subroutine */ int slabad_(real *, real *);
    extern doublereal slamch_(char *);
    extern /* Subroutine */ int sgesdd_(char *, integer *, integer *, real *, 
            integer *, real *, real *, integer *, real *, integer *, real *, 
            integer *, integer *, integer *), xerbla_(char *, integer 
            *);
    integer ioldsd[4];
    extern /* Subroutine */ int alasvm_(char *, integer *, integer *, integer 
            *, integer *), sgesvd_(char *, char *, integer *, integer 
            *, real *, integer *, real *, real *, integer *, real *, integer *
, real *, integer *, integer *), slacpy_(char *, 
            integer *, integer *, real *, integer *, real *, integer *), slaset_(char *, integer *, integer *, real *, real *, 
            real *, integer *), sgesvj_(char *, char *, char *, 
            integer *, integer *, real *, integer *, real *, integer *, real *
, integer *, real *, integer *, integer *)
            , sgejsv_(char *, char *, char *, char *, char *, char *, integer 
            *, integer *, real *, integer *, real *, real *, integer *, real *
, integer *, real *, integer *, integer *, integer *), slatms_(integer *, 
            integer *, char *, integer *, char *, real *, integer *, real *, 
            real *, integer *, integer *, char *, real *, integer *, real *, 
            integer *);
    integer minwrk;
    real ulpinv, result[22];
    integer lswork, mtypes;

    /* Fortran I/O blocks */
    static cilist io___25 = { 0, 0, 0, fmt_9996, 0 };
    static cilist io___30 = { 0, 0, 0, fmt_9995, 0 };
    static cilist io___38 = { 0, 0, 0, fmt_9995, 0 };
    static cilist io___41 = { 0, 0, 0, fmt_9995, 0 };
    static cilist io___42 = { 0, 0, 0, fmt_9995, 0 };
    static cilist io___43 = { 0, 0, 0, fmt_9999, 0 };
    static cilist io___44 = { 0, 0, 0, fmt_9998, 0 };
    static cilist io___45 = { 0, 0, 0, fmt_9997, 0 };



/*  -- LAPACK test routine (version 3.1) -- */
/*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
/*     November 2006 */

/*     .. Scalar Arguments .. */
/*     .. */
/*     .. Array Arguments .. */
/*     .. */

/*  Purpose */
/*  ======= */

/*  SDRVBD checks the singular value decomposition (SVD) drivers */
/*  SGESVD, SGESDD, SGESVJ, and SGEJSV. */

/*  Both SGESVD and SGESDD factor A = U diag(S) VT, where U and VT are */
/*  orthogonal and diag(S) is diagonal with the entries of the array S */
/*  on its diagonal. The entries of S are the singular values, */
/*  nonnegative and stored in decreasing order.  U and VT can be */
/*  optionally not computed, overwritten on A, or computed partially. */

/*  A is M by N. Let MNMIN = min( M, N ). S has dimension MNMIN. */
/*  U can be M by M or M by MNMIN. VT can be N by N or MNMIN by N. */

/*  When SDRVBD is called, a number of matrix "sizes" (M's and N's) */
/*  and a number of matrix "types" are specified.  For each size (M,N) */
/*  and each type of matrix, and for the minimal workspace as well as */
/*  workspace adequate to permit blocking, an  M x N  matrix "A" will be */
/*  generated and used to test the SVD routines.  For each matrix, A will */
/*  be factored as A = U diag(S) VT and the following 12 tests computed: */

/*  Test for SGESVD: */

/*  (1)    | A - U diag(S) VT | / ( |A| max(M,N) ulp ) */

/*  (2)    | I - U'U | / ( M ulp ) */

/*  (3)    | I - VT VT' | / ( N ulp ) */

/*  (4)    S contains MNMIN nonnegative values in decreasing order. */
/*         (Return 0 if true, 1/ULP if false.) */

/*  (5)    | U - Upartial | / ( M ulp ) where Upartial is a partially */
/*         computed U. */

/*  (6)    | VT - VTpartial | / ( N ulp ) where VTpartial is a partially */
/*         computed VT. */

/*  (7)    | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the */
/*         vector of singular values from the partial SVD */

/*  Test for SGESDD: */

/*  (8)    | A - U diag(S) VT | / ( |A| max(M,N) ulp ) */

/*  (9)    | I - U'U | / ( M ulp ) */

/*  (10)   | I - VT VT' | / ( N ulp ) */

/*  (11)   S contains MNMIN nonnegative values in decreasing order. */
/*         (Return 0 if true, 1/ULP if false.) */

/*  (12)   | U - Upartial | / ( M ulp ) where Upartial is a partially */
/*         computed U. */

/*  (13)   | VT - VTpartial | / ( N ulp ) where VTpartial is a partially */
/*         computed VT. */

/*  (14)   | S - Spartial | / ( MNMIN ulp |S| ) where Spartial is the */
/*         vector of singular values from the partial SVD */

/*  Test for SGESVJ: */

/*  (15)    | A - U diag(S) VT | / ( |A| max(M,N) ulp ) */

/*  (16)    | I - U'U | / ( M ulp ) */

/*  (17)   | I - VT VT' | / ( N ulp ) */

/*  (18)   S contains MNMIN nonnegative values in decreasing order. */
/*         (Return 0 if true, 1/ULP if false.) */

/*  Test for SGEJSV: */

/*  (19)    | A - U diag(S) VT | / ( |A| max(M,N) ulp ) */

/*  (20)    | I - U'U | / ( M ulp ) */

/*  (21)   | I - VT VT' | / ( N ulp ) */

/*  (22)   S contains MNMIN nonnegative values in decreasing order. */
/*         (Return 0 if true, 1/ULP if false.) */

/*  The "sizes" are specified by the arrays MM(1:NSIZES) and */
/*  NN(1:NSIZES); the value of each element pair (MM(j),NN(j)) */
/*  specifies one size.  The "types" are specified by a logical array */
/*  DOTYPE( 1:NTYPES ); if DOTYPE(j) is .TRUE., then matrix type "j" */
/*  will be generated. */
/*  Currently, the list of possible types is: */

/*  (1)  The zero matrix. */
/*  (2)  The identity matrix. */
/*  (3)  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. */
/*  (4)  Same as (3), but multiplied by the underflow-threshold / ULP. */
/*  (5)  Same as (3), but multiplied by the overflow-threshold * ULP. */

/*  Arguments */
/*  ========== */

/*  NSIZES  (input) INTEGER */
/*          The number of matrix sizes (M,N) contained in the vectors */
/*          MM and NN. */

/*  MM      (input) INTEGER array, dimension (NSIZES) */
/*          The values of the matrix row dimension M. */

/*  NN      (input) INTEGER array, dimension (NSIZES) */
/*          The values of the matrix column dimension N. */

/*  NTYPES  (input) INTEGER */
/*          The number of elements in DOTYPE.   If it is zero, SDRVBD */
/*          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. . */

/*  DOTYPE  (input) 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. */

/*  ISEED   (input/output) INTEGER array, dimension (4) */
/*          On entry, 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. */
/*          On exit, ISEED is changed and can be used in the next call to */
/*          SDRVBD to continue the same random number sequence. */

/*  THRESH  (input) REAL */
/*          The threshold value for the test ratios.  A result is */
/*          included in the output file if RESULT >= THRESH.  The test */
/*          ratios are scaled to be O(1), so THRESH should be a small */
/*          multiple of 1, e.g., 10 or 100.  To have every test ratio */
/*          printed, use THRESH = 0. */

/*  A       (workspace) REAL array, dimension (LDA,NMAX) */
/*          where NMAX is the maximum value of N in NN. */

/*  LDA     (input) INTEGER */
/*          The leading dimension of the array A.  LDA >= max(1,MMAX), */
/*          where MMAX is the maximum value of M in MM. */

/*  U       (workspace) REAL array, dimension (LDU,MMAX) */

/*  LDU     (input) INTEGER */
/*          The leading dimension of the array U.  LDU >= max(1,MMAX). */

/*  VT      (workspace) REAL array, dimension (LDVT,NMAX) */

/*  LDVT    (input) INTEGER */
/*          The leading dimension of the array VT.  LDVT >= max(1,NMAX). */

/*  ASAV    (workspace) REAL array, dimension (LDA,NMAX) */

/*  USAV    (workspace) REAL array, dimension (LDU,MMAX) */

/*  VTSAV   (workspace) REAL array, dimension (LDVT,NMAX) */

/*  S       (workspace) REAL array, dimension */
/*                      (max(min(MM,NN))) */

/*  SSAV    (workspace) REAL array, dimension */
/*                      (max(min(MM,NN))) */

/*  E       (workspace) REAL array, dimension */
/*                      (max(min(MM,NN))) */

/*  WORK    (workspace) REAL array, dimension (LWORK) */

/*  LWORK   (input) INTEGER */
/*          The number of entries in WORK.  This must be at least */
/*          max(3*MN+MX,5*MN-4)+2*MN**2 for all pairs */
/*          pairs  (MN,MX)=( min(MM(j),NN(j), max(MM(j),NN(j)) ) */

/*  IWORK   (workspace) INTEGER array, dimension at least 8*min(M,N) */

/*  NOUT    (input) INTEGER */
/*          The FORTRAN unit number for printing out error messages */
/*          (e.g., if a routine returns IINFO not equal to 0.) */

/*  INFO    (output) INTEGER */
/*          If 0, then everything ran OK. */
/*           -1: NSIZES < 0 */
/*           -2: Some MM(j) < 0 */
/*           -3: Some NN(j) < 0 */
/*           -4: NTYPES < 0 */
/*           -7: THRESH < 0 */
/*          -10: LDA < 1 or LDA < MMAX, where MMAX is max( MM(j) ). */
/*          -12: LDU < 1 or LDU < MMAX. */
/*          -14: LDVT < 1 or LDVT < NMAX, where NMAX is max( NN(j) ). */
/*          -21: LWORK too small. */
/*          If  SLATMS, or SGESVD returns an error code, the */
/*              absolute value of it is returned. */

/*  ===================================================================== */

/*     .. Parameters .. */
/*     .. */
/*     .. Local Scalars .. */
/*     .. */
/*     .. Local Arrays .. */
/*     .. */
/*     .. External Functions .. */
/*     .. */
/*     .. External Subroutines .. */
/*     .. */
/*     .. Intrinsic Functions .. */
/*     .. */
/*     .. Scalars in Common .. */
/*     .. */
/*     .. Common blocks .. */
/*     .. */
/*     .. Data statements .. */
    /* Parameter adjustments */
    --mm;
    --nn;
    --dotype;
    --iseed;
    asav_dim1 = *lda;
    asav_offset = 1 + asav_dim1;
    asav -= asav_offset;
    a_dim1 = *lda;
    a_offset = 1 + a_dim1;
    a -= a_offset;
    usav_dim1 = *ldu;
    usav_offset = 1 + usav_dim1;
    usav -= usav_offset;
    u_dim1 = *ldu;
    u_offset = 1 + u_dim1;
    u -= u_offset;
    vtsav_dim1 = *ldvt;
    vtsav_offset = 1 + vtsav_dim1;
    vtsav -= vtsav_offset;
    vt_dim1 = *ldvt;
    vt_offset = 1 + vt_dim1;
    vt -= vt_offset;
    --s;
    --ssav;
    --e;
    --work;
    --iwork;

    /* Function Body */
/*     .. */
/*     .. Executable Statements .. */

/*     Check for errors */

    *info = 0;
    badmm = FALSE_;
    badnn = FALSE_;
    mmax = 1;
    nmax = 1;
    mnmax = 1;
    minwrk = 1;
    i__1 = *nsizes;
    for (j = 1; j <= i__1; ++j) {
/* Computing MAX */
        i__2 = mmax, i__3 = mm[j];
        mmax = max(i__2,i__3);
        if (mm[j] < 0) {
            badmm = TRUE_;
        }
/* Computing MAX */
        i__2 = nmax, i__3 = nn[j];
        nmax = max(i__2,i__3);
        if (nn[j] < 0) {
            badnn = TRUE_;
        }
/* Computing MAX */
/* Computing MIN */
        i__4 = mm[j], i__5 = nn[j];
        i__2 = mnmax, i__3 = min(i__4,i__5);
        mnmax = max(i__2,i__3);
/* Computing MAX */
/* Computing MAX */
/* Computing MIN */
        i__6 = mm[j], i__7 = nn[j];
/* Computing MAX */
        i__8 = mm[j], i__9 = nn[j];
/* Computing MIN */
        i__10 = mm[j], i__11 = nn[j] - 4;
        i__4 = min(i__6,i__7) * 3 + max(i__8,i__9), i__5 = min(i__10,i__11) * 
                5;
/* Computing MIN */
        i__13 = mm[j], i__14 = nn[j];
/* Computing 2nd power */
        i__12 = min(i__13,i__14);
        i__2 = minwrk, i__3 = max(i__4,i__5) + (i__12 * i__12 << 1);
        minwrk = max(i__2,i__3);
/* L10: */
    }

/*     Check for errors */

    if (*nsizes < 0) {
        *info = -1;
    } else if (badmm) {
        *info = -2;
    } else if (badnn) {
        *info = -3;
    } else if (*ntypes < 0) {
        *info = -4;
    } else if (*lda < max(1,mmax)) {
        *info = -10;
    } else if (*ldu < max(1,mmax)) {
        *info = -12;
    } else if (*ldvt < max(1,nmax)) {
        *info = -14;
    } else if (minwrk > *lwork) {
        *info = -21;
    }

    if (*info != 0) {
        i__1 = -(*info);
        xerbla_("SDRVBD", &i__1);
        return 0;
    }

/*     Initialize constants */

    s_copy(path, "Single precision", (ftnlen)1, (ftnlen)16);
    s_copy(path + 1, "BD", (ftnlen)2, (ftnlen)2);
    nfail = 0;
    ntest = 0;
    unfl = slamch_("Safe minimum");
    ovfl = 1.f / unfl;
    slabad_(&unfl, &ovfl);
    ulp = slamch_("Precision");
    ulpinv = 1.f / ulp;
    infoc_1.infot = 0;

/*     Loop over sizes, types */

    i__1 = *nsizes;
    for (jsize = 1; jsize <= i__1; ++jsize) {
        m = mm[jsize];
        n = nn[jsize];
        mnmin = min(m,n);

        if (*nsizes != 1) {
            mtypes = min(5,*ntypes);
        } else {
            mtypes = min(6,*ntypes);
        }

        i__2 = mtypes;
        for (jtype = 1; jtype <= i__2; ++jtype) {
            if (! dotype[jtype]) {
                goto L140;
            }

            for (j = 1; j <= 4; ++j) {
                ioldsd[j - 1] = iseed[j];
/* L20: */
            }

/*           Compute "A" */

            if (mtypes > 5) {
                goto L30;
            }

            if (jtype == 1) {

/*              Zero matrix */

                slaset_("Full", &m, &n, &c_b13, &c_b13, &a[a_offset], lda);

            } else if (jtype == 2) {

/*              Identity matrix */

                slaset_("Full", &m, &n, &c_b13, &c_b17, &a[a_offset], lda);

            } else {

/*              (Scaled) random matrix */

                if (jtype == 3) {
                    anorm = 1.f;
                }
                if (jtype == 4) {
                    anorm = unfl / ulp;
                }
                if (jtype == 5) {
                    anorm = ovfl * ulp;
                }
                r__1 = (real) mnmin;
                i__3 = m - 1;
                i__4 = n - 1;
                slatms_(&m, &n, "U", &iseed[1], "N", &s[1], &c__4, &r__1, &
                        anorm, &i__3, &i__4, "N", &a[a_offset], lda, &work[1], 
                         &iinfo);
                if (iinfo != 0) {
                    io___25.ciunit = *nout;
                    s_wsfe(&io___25);
                    do_fio(&c__1, "Generator", (ftnlen)9);
                    do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
                    do_fio(&c__1, (char *)&m, (ftnlen)sizeof(integer));
                    do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
                    do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
                    do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer))
                            ;
                    e_wsfe();
                    *info = abs(iinfo);
                    return 0;
                }
            }

L30:
            slacpy_("F", &m, &n, &a[a_offset], lda, &asav[asav_offset], lda);

/*           Do for minimal and adequate (for blocking) workspace */

            for (iws = 1; iws <= 4; ++iws) {

                for (j = 1; j <= 14; ++j) {
                    result[j - 1] = -1.f;
/* L40: */
                }

/*              Test SGESVD: Factorize A */

/* Computing MAX */
                i__3 = min(m,n) * 3 + max(m,n), i__4 = min(m,n) * 5;
                iwtmp = max(i__3,i__4);
                lswork = iwtmp + (iws - 1) * (*lwork - iwtmp) / 3;
                lswork = min(lswork,*lwork);
                lswork = max(lswork,1);
                if (iws == 4) {
                    lswork = *lwork;
                }

                if (iws > 1) {
                    slacpy_("F", &m, &n, &asav[asav_offset], lda, &a[a_offset]
, lda);
                }
                s_copy(srnamc_1.srnamt, "SGESVD", (ftnlen)32, (ftnlen)6);
                sgesvd_("A", "A", &m, &n, &a[a_offset], lda, &ssav[1], &usav[
                        usav_offset], ldu, &vtsav[vtsav_offset], ldvt, &work[
                        1], &lswork, &iinfo);
                if (iinfo != 0) {
                    io___30.ciunit = *nout;
                    s_wsfe(&io___30);
                    do_fio(&c__1, "GESVD", (ftnlen)5);
                    do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
                    do_fio(&c__1, (char *)&m, (ftnlen)sizeof(integer));
                    do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
                    do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
                    do_fio(&c__1, (char *)&lswork, (ftnlen)sizeof(integer));
                    do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer))
                            ;
                    e_wsfe();
                    *info = abs(iinfo);
                    return 0;
                }

/*              Do tests 1--4 */

                sbdt01_(&m, &n, &c__0, &asav[asav_offset], lda, &usav[
                        usav_offset], ldu, &ssav[1], &e[1], &vtsav[
                        vtsav_offset], ldvt, &work[1], result);
                if (m != 0 && n != 0) {
                    sort01_("Columns", &m, &m, &usav[usav_offset], ldu, &work[
                            1], lwork, &result[1]);
                    sort01_("Rows", &n, &n, &vtsav[vtsav_offset], ldvt, &work[
                            1], lwork, &result[2]);
                }
                result[3] = 0.f;
                i__3 = mnmin - 1;
                for (i__ = 1; i__ <= i__3; ++i__) {
                    if (ssav[i__] < ssav[i__ + 1]) {
                        result[3] = ulpinv;
                    }
                    if (ssav[i__] < 0.f) {
                        result[3] = ulpinv;
                    }
/* L50: */
                }
                if (mnmin >= 1) {
                    if (ssav[mnmin] < 0.f) {
                        result[3] = ulpinv;
                    }
                }

/*              Do partial SVDs, comparing to SSAV, USAV, and VTSAV */

                result[4] = 0.f;
                result[5] = 0.f;
                result[6] = 0.f;
                for (iju = 0; iju <= 3; ++iju) {
                    for (ijvt = 0; ijvt <= 3; ++ijvt) {
                        if (iju == 3 && ijvt == 3 || iju == 1 && ijvt == 1) {
                            goto L70;
                        }
                        *(unsigned char *)jobu = *(unsigned char *)&cjob[iju];
                        *(unsigned char *)jobvt = *(unsigned char *)&cjob[
                                ijvt];
                        slacpy_("F", &m, &n, &asav[asav_offset], lda, &a[
                                a_offset], lda);
                        s_copy(srnamc_1.srnamt, "SGESVD", (ftnlen)32, (ftnlen)
                                6);
                        sgesvd_(jobu, jobvt, &m, &n, &a[a_offset], lda, &s[1], 
                                 &u[u_offset], ldu, &vt[vt_offset], ldvt, &
                                work[1], &lswork, &iinfo);

/*                    Compare U */

                        dif = 0.f;
                        if (m > 0 && n > 0) {
                            if (iju == 1) {
                                sort03_("C", &m, &mnmin, &m, &mnmin, &usav[
                                        usav_offset], ldu, &a[a_offset], lda, 
                                        &work[1], lwork, &dif, &iinfo);
                            } else if (iju == 2) {
                                sort03_("C", &m, &mnmin, &m, &mnmin, &usav[
                                        usav_offset], ldu, &u[u_offset], ldu, 
                                        &work[1], lwork, &dif, &iinfo);
                            } else if (iju == 3) {
                                sort03_("C", &m, &m, &m, &mnmin, &usav[
                                        usav_offset], ldu, &u[u_offset], ldu, 
                                        &work[1], lwork, &dif, &iinfo);
                            }
                        }
                        result[4] = dmax(result[4],dif);

/*                    Compare VT */

                        dif = 0.f;
                        if (m > 0 && n > 0) {
                            if (ijvt == 1) {
                                sort03_("R", &n, &mnmin, &n, &mnmin, &vtsav[
                                        vtsav_offset], ldvt, &a[a_offset], 
                                        lda, &work[1], lwork, &dif, &iinfo);
                            } else if (ijvt == 2) {
                                sort03_("R", &n, &mnmin, &n, &mnmin, &vtsav[
                                        vtsav_offset], ldvt, &vt[vt_offset], 
                                        ldvt, &work[1], lwork, &dif, &iinfo);
                            } else if (ijvt == 3) {
                                sort03_("R", &n, &n, &n, &mnmin, &vtsav[
                                        vtsav_offset], ldvt, &vt[vt_offset], 
                                        ldvt, &work[1], lwork, &dif, &iinfo);
                            }
                        }
                        result[5] = dmax(result[5],dif);

/*                    Compare S */

                        dif = 0.f;
/* Computing MAX */
                        r__1 = (real) mnmin * ulp * s[1];
                        div = dmax(r__1,unfl);
                        i__3 = mnmin - 1;
                        for (i__ = 1; i__ <= i__3; ++i__) {
                            if (ssav[i__] < ssav[i__ + 1]) {
                                dif = ulpinv;
                            }
                            if (ssav[i__] < 0.f) {
                                dif = ulpinv;
                            }
/* Computing MAX */
                            r__2 = dif, r__3 = (r__1 = ssav[i__] - s[i__], 
                                    dabs(r__1)) / div;
                            dif = dmax(r__2,r__3);
/* L60: */
                        }
                        result[6] = dmax(result[6],dif);
L70:
                        ;
                    }
/* L80: */
                }

/*              Test SGESDD: Factorize A */

                iwtmp = mnmin * 5 * mnmin + mnmin * 9 + max(m,n);
                lswork = iwtmp + (iws - 1) * (*lwork - iwtmp) / 3;
                lswork = min(lswork,*lwork);
                lswork = max(lswork,1);
                if (iws == 4) {
                    lswork = *lwork;
                }

                slacpy_("F", &m, &n, &asav[asav_offset], lda, &a[a_offset], 
                        lda);
                s_copy(srnamc_1.srnamt, "SGESDD", (ftnlen)32, (ftnlen)6);
                sgesdd_("A", &m, &n, &a[a_offset], lda, &ssav[1], &usav[
                        usav_offset], ldu, &vtsav[vtsav_offset], ldvt, &work[
                        1], &lswork, &iwork[1], &iinfo);
                if (iinfo != 0) {
                    io___38.ciunit = *nout;
                    s_wsfe(&io___38);
                    do_fio(&c__1, "GESDD", (ftnlen)5);
                    do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer));
                    do_fio(&c__1, (char *)&m, (ftnlen)sizeof(integer));
                    do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
                    do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer));
                    do_fio(&c__1, (char *)&lswork, (ftnlen)sizeof(integer));
                    do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(integer))
                            ;
                    e_wsfe();
                    *info = abs(iinfo);
                    return 0;
                }

/*              Do tests 8--11 */

                sbdt01_(&m, &n, &c__0, &asav[asav_offset], lda, &usav[
                        usav_offset], ldu, &ssav[1], &e[1], &vtsav[
                        vtsav_offset], ldvt, &work[1], &result[7]);
                if (m != 0 && n != 0) {
                    sort01_("Columns", &m, &m, &usav[usav_offset], ldu, &work[
                            1], lwork, &result[8]);
                    sort01_("Rows", &n, &n, &vtsav[vtsav_offset], ldvt, &work[
                            1], lwork, &result[9]);
                }
                result[10] = 0.f;
                i__3 = mnmin - 1;
                for (i__ = 1; i__ <= i__3; ++i__) {
                    if (ssav[i__] < ssav[i__ + 1]) {
                        result[10] = ulpinv;
                    }
                    if (ssav[i__] < 0.f) {
                        result[10] = ulpinv;
                    }
/* L90: */
                }
                if (mnmin >= 1) {
                    if (ssav[mnmin] < 0.f) {
                        result[10] = ulpinv;
                    }
                }

/*              Do partial SVDs, comparing to SSAV, USAV, and VTSAV */

                result[11] = 0.f;
                result[12] = 0.f;
                result[13] = 0.f;
                for (ijq = 0; ijq <= 2; ++ijq) {
                    *(unsigned char *)jobq = *(unsigned char *)&cjob[ijq];
                    slacpy_("F", &m, &n, &asav[asav_offset], lda, &a[a_offset]
, lda);
                    s_copy(srnamc_1.srnamt, "SGESDD", (ftnlen)32, (ftnlen)6);
                    sgesdd_(jobq, &m, &n, &a[a_offset], lda, &s[1], &u[
                            u_offset], ldu, &vt[vt_offset], ldvt, &work[1], &
                            lswork, &iwork[1], &iinfo);

/*                 Compare U */

                    dif = 0.f;
                    if (m > 0 && n > 0) {
                        if (ijq == 1) {
                            if (m >= n) {
                                sort03_("C", &m, &mnmin, &m, &mnmin, &usav[
                                        usav_offset], ldu, &a[a_offset], lda, 
                                        &work[1], lwork, &dif, info);
                            } else {
                                sort03_("C", &m, &mnmin, &m, &mnmin, &usav[
                                        usav_offset], ldu, &u[u_offset], ldu, 
                                        &work[1], lwork, &dif, info);
                            }
                        } else if (ijq == 2) {
                            sort03_("C", &m, &mnmin, &m, &mnmin, &usav[
                                    usav_offset], ldu, &u[u_offset], ldu, &
                                    work[1], lwork, &dif, info);
                        }
                    }
                    result[11] = dmax(result[11],dif);

/*                 Compare VT */

                    dif = 0.f;
                    if (m > 0 && n > 0) {
                        if (ijq == 1) {
                            if (m >= n) {
                                sort03_("R", &n, &mnmin, &n, &mnmin, &vtsav[
                                        vtsav_offset], ldvt, &vt[vt_offset], 
                                        ldvt, &work[1], lwork, &dif, info);
                            } else {
                                sort03_("R", &n, &mnmin, &n, &mnmin, &vtsav[
                                        vtsav_offset], ldvt, &a[a_offset], 
                                        lda, &work[1], lwork, &dif, info);
                            }
                        } else if (ijq == 2) {
                            sort03_("R", &n, &mnmin, &n, &mnmin, &vtsav[
                                    vtsav_offset], ldvt, &vt[vt_offset], ldvt, 
                                     &work[1], lwork, &dif, info);
                        }
                    }
                    result[12] = dmax(result[12],dif);

/*                 Compare S */

                    dif = 0.f;
/* Computing MAX */
                    r__1 = (real) mnmin * ulp * s[1];
                    div = dmax(r__1,unfl);
                    i__3 = mnmin - 1;
                    for (i__ = 1; i__ <= i__3; ++i__) {
                        if (ssav[i__] < ssav[i__ + 1]) {
                            dif = ulpinv;
                        }
                        if (ssav[i__] < 0.f) {
                            dif = ulpinv;
                        }
/* Computing MAX */
                        r__2 = dif, r__3 = (r__1 = ssav[i__] - s[i__], dabs(
                                r__1)) / div;
                        dif = dmax(r__2,r__3);
/* L100: */
                    }
                    result[13] = dmax(result[13],dif);
/* L110: */
                }

/*              Test SGESVJ: Factorize A */
/*              Note: SGESVJ does not work for M < N */

                result[14] = 0.f;
                result[15] = 0.f;
                result[16] = 0.f;
                result[17] = 0.f;

                if (m >= n) {
                    iwtmp = mnmin * 5 * mnmin + mnmin * 9 + max(m,n);
                    lswork = iwtmp + (iws - 1) * (*lwork - iwtmp) / 3;
                    lswork = min(lswork,*lwork);
                    lswork = max(lswork,1);
                    if (iws == 4) {
                        lswork = *lwork;
                    }

                    slacpy_("F", &m, &n, &asav[asav_offset], lda, &usav[
                            usav_offset], lda);
                    s_copy(srnamc_1.srnamt, "SGESVJ", (ftnlen)32, (ftnlen)6);
                    sgesvj_("G", "U", "V", &m, &n, &usav[usav_offset], lda, &
                            ssav[1], &c__0, &a[a_offset], ldvt, &work[1], 
                            lwork, info);

/*                 SGESVJ retuns V not VT, so we transpose to use the same */
/*                 test suite. */

                    i__3 = n;
                    for (j = 1; j <= i__3; ++j) {
                        i__4 = n;
                        for (i__ = 1; i__ <= i__4; ++i__) {
                            vtsav[j + i__ * vtsav_dim1] = a[i__ + j * a_dim1];
                        }
                    }

                    if (iinfo != 0) {
                        io___41.ciunit = *nout;
                        s_wsfe(&io___41);
                        do_fio(&c__1, "GESVJ", (ftnlen)5);
                        do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer))
                                ;
                        do_fio(&c__1, (char *)&m, (ftnlen)sizeof(integer));
                        do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
                        do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer))
                                ;
                        do_fio(&c__1, (char *)&lswork, (ftnlen)sizeof(integer)
                                );
                        do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(
                                integer));
                        e_wsfe();
                        *info = abs(iinfo);
                        return 0;
                    }

/*                 Do tests 15--18 */

                    sbdt01_(&m, &n, &c__0, &asav[asav_offset], lda, &usav[
                            usav_offset], ldu, &ssav[1], &e[1], &vtsav[
                            vtsav_offset], ldvt, &work[1], &result[14]);
                    if (m != 0 && n != 0) {
                        sort01_("Columns", &m, &m, &usav[usav_offset], ldu, &
                                work[1], lwork, &result[15]);
                        sort01_("Rows", &n, &n, &vtsav[vtsav_offset], ldvt, &
                                work[1], lwork, &result[16]);
                    }
                    result[17] = 0.f;
                    i__3 = mnmin - 1;
                    for (i__ = 1; i__ <= i__3; ++i__) {
                        if (ssav[i__] < ssav[i__ + 1]) {
                            result[17] = ulpinv;
                        }
                        if (ssav[i__] < 0.f) {
                            result[17] = ulpinv;
                        }
/* L200: */
                    }
                    if (mnmin >= 1) {
                        if (ssav[mnmin] < 0.f) {
                            result[17] = ulpinv;
                        }
                    }
                }

/*              Test SGEJSV: Factorize A */
/*              Note: SGEJSV does not work for M < N */

                result[18] = 0.f;
                result[19] = 0.f;
                result[20] = 0.f;
                result[21] = 0.f;
                if (m >= n) {
                    iwtmp = mnmin * 5 * mnmin + mnmin * 9 + max(m,n);
                    lswork = iwtmp + (iws - 1) * (*lwork - iwtmp) / 3;
                    lswork = min(lswork,*lwork);
                    lswork = max(lswork,1);
                    if (iws == 4) {
                        lswork = *lwork;
                    }

                    slacpy_("F", &m, &n, &asav[asav_offset], lda, &vtsav[
                            vtsav_offset], lda);
                    s_copy(srnamc_1.srnamt, "SGEJSV", (ftnlen)32, (ftnlen)6);
                    sgejsv_("G", "U", "V", "R", "N", "N", &m, &n, &vtsav[
                            vtsav_offset], lda, &ssav[1], &usav[usav_offset], 
                            ldu, &a[a_offset], ldvt, &work[1], lwork, &iwork[
                            1], info);

/*                 SGEJSV retuns V not VT, so we transpose to use the same */
/*                 test suite. */

                    i__3 = n;
                    for (j = 1; j <= i__3; ++j) {
                        i__4 = n;
                        for (i__ = 1; i__ <= i__4; ++i__) {
                            vtsav[j + i__ * vtsav_dim1] = a[i__ + j * a_dim1];
                        }
                    }

                    if (iinfo != 0) {
                        io___42.ciunit = *nout;
                        s_wsfe(&io___42);
                        do_fio(&c__1, "GESVJ", (ftnlen)5);
                        do_fio(&c__1, (char *)&iinfo, (ftnlen)sizeof(integer))
                                ;
                        do_fio(&c__1, (char *)&m, (ftnlen)sizeof(integer));
                        do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
                        do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer))
                                ;
                        do_fio(&c__1, (char *)&lswork, (ftnlen)sizeof(integer)
                                );
                        do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(
                                integer));
                        e_wsfe();
                        *info = abs(iinfo);
                        return 0;
                    }

/*                 Do tests 19--22 */

                    sbdt01_(&m, &n, &c__0, &asav[asav_offset], lda, &usav[
                            usav_offset], ldu, &ssav[1], &e[1], &vtsav[
                            vtsav_offset], ldvt, &work[1], &result[18]);
                    if (m != 0 && n != 0) {
                        sort01_("Columns", &m, &m, &usav[usav_offset], ldu, &
                                work[1], lwork, &result[19]);
                        sort01_("Rows", &n, &n, &vtsav[vtsav_offset], ldvt, &
                                work[1], lwork, &result[20]);
                    }
                    result[21] = 0.f;
                    i__3 = mnmin - 1;
                    for (i__ = 1; i__ <= i__3; ++i__) {
                        if (ssav[i__] < ssav[i__ + 1]) {
                            result[21] = ulpinv;
                        }
                        if (ssav[i__] < 0.f) {
                            result[21] = ulpinv;
                        }
/* L300: */
                    }
                    if (mnmin >= 1) {
                        if (ssav[mnmin] < 0.f) {
                            result[21] = ulpinv;
                        }
                    }
                }

/*              End of Loop -- Check for RESULT(j) > THRESH */

                for (j = 1; j <= 22; ++j) {
                    if (result[j - 1] >= *thresh) {
                        if (nfail == 0) {
                            io___43.ciunit = *nout;
                            s_wsfe(&io___43);
                            e_wsfe();
                            io___44.ciunit = *nout;
                            s_wsfe(&io___44);
                            e_wsfe();
                        }
                        io___45.ciunit = *nout;
                        s_wsfe(&io___45);
                        do_fio(&c__1, (char *)&m, (ftnlen)sizeof(integer));
                        do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
                        do_fio(&c__1, (char *)&jtype, (ftnlen)sizeof(integer))
                                ;
                        do_fio(&c__1, (char *)&iws, (ftnlen)sizeof(integer));
                        do_fio(&c__4, (char *)&ioldsd[0], (ftnlen)sizeof(
                                integer));
                        do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer));
                        do_fio(&c__1, (char *)&result[j - 1], (ftnlen)sizeof(
                                real));
                        e_wsfe();
                        ++nfail;
                    }
/* L120: */
                }
                ntest += 22;

/* L130: */
            }
L140:
            ;
        }
/* L150: */
    }

/*     Summary */

    alasvm_(path, nout, &nfail, &ntest, &c__0);


    return 0;

/*     End of SDRVBD */

} /* sdrvbd_ */