#include "f2c.h"
#include "blaswrap.h"
/* Table of constant values */
static doublecomplex c_b1 = {0.,0.};
static doublecomplex c_b2 = {1.,0.};
static integer c__4 = 4;
static integer c__1 = 1;
static integer c__0 = 0;
/* Subroutine */ int zdrvbd_(integer *nsizes, integer *mm, integer *nn,
integer *ntypes, logical *dotype, integer *iseed, doublereal *thresh,
doublecomplex *a, integer *lda, doublecomplex *u, integer *ldu,
doublecomplex *vt, integer *ldvt, doublecomplex *asav, doublecomplex *
usav, doublecomplex *vtsav, doublereal *s, doublereal *ssav,
doublereal *e, doublecomplex *work, integer *lwork, doublereal *rwork,
integer *iwork, integer *nounit, integer *info)
{
/* Initialized data */
static char cjob[1*4] = "N" "O" "S" "A";
/* Format strings */
static char fmt_9996[] = "(\002 ZDRVBD: \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 ZDRVBD: \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 -- Complex Singular Value Decomposit"
"ion Driver \002,/\002 Matrix types (see ZDRVBD 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 sp"
"aced singular values near underflow\002,/\002 5 = Evenly spaced "
"singular values near overflow\002,//\002 Tests performed: ( A is"
" dense, U and V are unitary,\002,/19x,\002 S is an array, and Up"
"artial, VTpartial, and\002,/19x,\002 Spartial are partially comp"
"uted U, VT and S),\002,/)";
static char fmt_9998[] = "(\002 Tests performed with Test Threshold ="
" \002,f8.2,/\002 ZGESVD: \002,/\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 o"
"rder, else 1/ulp\002,/\002 5 = | U - Upartial | / ( M ulp )\002,/"
"\002 6 = | VT - VTpartial | / ( N ulp )\002,/\002 7 = | S - Spar"
"tial | / ( min(M,N) ulp |S| )\002,/\002 ZGESDD: \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 min(M,N) nonnegative values in"
"\002,\002 decreasing order, else 1/ulp\002,/\00212 = | U - Upart"
"ial | / ( M ulp )\002,/\00213 = | VT - VTpartial | / ( N ulp "
")\002,/\00214 = | 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,i1,\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;
doublereal d__1, d__2, d__3;
/* Builtin functions */
integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
/* Local variables */
integer i__, j, m, n;
doublereal dif, div;
integer ijq, iju;
doublereal ulp;
char jobq[1], jobu[1];
integer mmax, nmax;
doublereal unfl, ovfl;
integer ijvt;
logical badmm, badnn;
integer nfail, iinfo;
extern /* Subroutine */ int zbdt01_(integer *, integer *, integer *,
doublecomplex *, integer *, doublecomplex *, integer *,
doublereal *, doublereal *, doublecomplex *, integer *,
doublecomplex *, doublereal *, doublereal *);
doublereal anorm;
integer mnmin, mnmax;
char jobvt[1];
integer iwspc, jsize, nerrs, jtype, ntest, iwtmp;
extern /* Subroutine */ int zunt01_(char *, integer *, integer *,
doublecomplex *, integer *, doublecomplex *, integer *,
doublereal *, doublereal *), zunt03_(char *, integer *,
integer *, integer *, integer *, doublecomplex *, integer *,
doublecomplex *, integer *, doublecomplex *, integer *,
doublereal *, doublereal *, integer *);
extern doublereal dlamch_(char *);
extern /* Subroutine */ int xerbla_(char *, integer *);
integer ioldsd[4];
extern /* Subroutine */ int zgesdd_(char *, integer *, integer *,
doublecomplex *, integer *, doublereal *, doublecomplex *,
integer *, doublecomplex *, integer *, doublecomplex *, integer *,
doublereal *, integer *, integer *), alasvm_(char *,
integer *, integer *, integer *, integer *), zgesvd_(char
*, char *, integer *, integer *, doublecomplex *, integer *,
doublereal *, doublecomplex *, integer *, doublecomplex *,
integer *, doublecomplex *, integer *, doublereal *, integer *);
integer ntestf;
extern /* Subroutine */ int zlacpy_(char *, integer *, integer *,
doublecomplex *, integer *, doublecomplex *, integer *),
zlaset_(char *, integer *, integer *, doublecomplex *,
doublecomplex *, doublecomplex *, integer *);
integer minwrk;
extern /* Subroutine */ int zlatms_(integer *, integer *, char *, integer
*, char *, doublereal *, integer *, doublereal *, doublereal *,
integer *, integer *, char *, doublecomplex *, integer *,
doublecomplex *, integer *);
doublereal ulpinv, result[14];
integer lswork, mtypes, ntestt;
/* Fortran I/O blocks */
static cilist io___27 = { 0, 0, 0, fmt_9996, 0 };
static cilist io___32 = { 0, 0, 0, fmt_9995, 0 };
static cilist io___39 = { 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 */
/* ======= */
/* ZDRVBD checks the singular value decomposition (SVD) driver ZGESVD */
/* and ZGESDD. */
/* ZGESVD and CGESDD factors A = U diag(S) VT, where U and VT are */
/* unitary 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 ZDRVBD 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 ZGESVD: */
/* (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 ZGESDD: */
/* (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 */
/* 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 unitary 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 sizes of matrices to use. If it is zero, */
/* ZDRVBD does nothing. It must be at least zero. */
/* MM (input) INTEGER array, dimension (NSIZES) */
/* An array containing the matrix "heights" to be used. For */
/* each j=1,...,NSIZES, if MM(j) is zero, then MM(j) and NN(j) */
/* will be ignored. The MM(j) values must be at least zero. */
/* NN (input) INTEGER array, dimension (NSIZES) */
/* An array containing the matrix "widths" to be used. For */
/* each j=1,...,NSIZES, if NN(j) is zero, then MM(j) and NN(j) */
/* will be ignored. The NN(j) values must be at least zero. */
/* NTYPES (input) INTEGER */
/* The number of elements in DOTYPE. If it is zero, ZDRVBD */
/* 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 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 random number generator uses a linear */
/* congruential sequence limited to small integers, and so */
/* should produce machine independent random numbers. The */
/* values of ISEED are changed on exit, and can be used in the */
/* next call to ZDRVBD to continue the same random number */
/* sequence. */
/* THRESH (input) DOUBLE PRECISION */
/* A test will count as "failed" if the "error", computed as */
/* described above, exceeds THRESH. Note that the error */
/* is scaled to be O(1), so THRESH should be a reasonably */
/* small multiple of 1, e.g., 10 or 100. In particular, */
/* it should not depend on the precision (single vs. double) */
/* or the size of the matrix. It must be at least zero. */
/* NOUNIT (input) INTEGER */
/* The FORTRAN unit number for printing out error messages */
/* (e.g., if a routine returns IINFO not equal to 0.) */
/* A (output) COMPLEX*16 array, dimension (LDA,max(NN)) */
/* Used to hold the matrix whose singular values are to be */
/* computed. On exit, A contains the last matrix actually */
/* used. */
/* LDA (input) INTEGER */
/* The leading dimension of A. It must be at */
/* least 1 and at least max( MM ). */
/* U (output) COMPLEX*16 array, dimension (LDU,max(MM)) */
/* Used to hold the computed matrix of right singular vectors. */
/* On exit, U contains the last such vectors actually computed. */
/* LDU (input) INTEGER */
/* The leading dimension of U. It must be at */
/* least 1 and at least max( MM ). */
/* VT (output) COMPLEX*16 array, dimension (LDVT,max(NN)) */
/* Used to hold the computed matrix of left singular vectors. */
/* On exit, VT contains the last such vectors actually computed. */
/* LDVT (input) INTEGER */
/* The leading dimension of VT. It must be at */
/* least 1 and at least max( NN ). */
/* ASAV (output) COMPLEX*16 array, dimension (LDA,max(NN)) */
/* Used to hold a different copy of the matrix whose singular */
/* values are to be computed. On exit, A contains the last */
/* matrix actually used. */
/* USAV (output) COMPLEX*16 array, dimension (LDU,max(MM)) */
/* Used to hold a different copy of the computed matrix of */
/* right singular vectors. On exit, USAV contains the last such */
/* vectors actually computed. */
/* VTSAV (output) COMPLEX*16 array, dimension (LDVT,max(NN)) */
/* Used to hold a different copy of the computed matrix of */
/* left singular vectors. On exit, VTSAV contains the last such */
/* vectors actually computed. */
/* S (output) DOUBLE PRECISION array, dimension (max(min(MM,NN))) */
/* Contains the computed singular values. */
/* SSAV (output) DOUBLE PRECISION array, dimension (max(min(MM,NN))) */
/* Contains another copy of the computed singular values. */
/* E (output) DOUBLE PRECISION array, dimension (max(min(MM,NN))) */
/* Workspace for ZGESVD. */
/* WORK (workspace) COMPLEX*16 array, dimension (LWORK) */
/* LWORK (input) INTEGER */
/* The number of entries in WORK. This must be at least */
/* MAX(3*MIN(M,N)+MAX(M,N)**2,5*MIN(M,N),3*MAX(M,N)) for all */
/* pairs (M,N)=(MM(j),NN(j)) */
/* RWORK (workspace) DOUBLE PRECISION array, */
/* dimension ( 5*max(max(MM,NN)) ) */
/* IWORK (workspace) INTEGER array, dimension at least 8*min(M,N) */
/* RESULT (output) DOUBLE PRECISION array, dimension (7) */
/* The values computed by the 7 tests described above. */
/* The values are currently limited to 1/ULP, to avoid */
/* overflow. */
/* 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 ZLATMS, or ZGESVD returns an error code, the */
/* absolute value of it is returned. */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. Local Arrays .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. 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;
--rwork;
--iwork;
/* Function Body */
/* .. */
/* .. Executable Statements .. */
/* Check for errors */
*info = 0;
/* Important constants */
nerrs = 0;
ntestt = 0;
ntestf = 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__9 = mm[j], i__10 = nn[j];
/* Computing 2nd power */
i__8 = max(i__9,i__10);
/* Computing MIN */
i__11 = mm[j], i__12 = nn[j];
/* Computing MAX */
i__13 = mm[j], i__14 = nn[j];
i__4 = min(i__6,i__7) * 3 + i__8 * i__8, i__5 = min(i__11,i__12) * 5,
i__4 = max(i__4,i__5), i__5 = max(i__13,i__14) * 3;
i__2 = minwrk, i__3 = max(i__4,i__5);
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_("ZDRVBD", &i__1);
return 0;
}
/* Quick return if nothing to do */
if (*nsizes == 0 || *ntypes == 0) {
return 0;
}
/* More Important constants */
unfl = dlamch_("S");
ovfl = 1. / unfl;
ulp = dlamch_("E");
ulpinv = 1. / ulp;
/* Loop over sizes, types */
nerrs = 0;
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 L170;
}
ntest = 0;
for (j = 1; j <= 4; ++j) {
ioldsd[j - 1] = iseed[j];
/* L20: */
}
/* Compute "A" */
if (mtypes > 5) {
goto L50;
}
if (jtype == 1) {
/* Zero matrix */
zlaset_("Full", &m, &n, &c_b1, &c_b1, &a[a_offset], lda);
i__3 = min(m,n);
for (i__ = 1; i__ <= i__3; ++i__) {
s[i__] = 0.;
/* L30: */
}
} else if (jtype == 2) {
/* Identity matrix */
zlaset_("Full", &m, &n, &c_b1, &c_b2, &a[a_offset], lda);
i__3 = min(m,n);
for (i__ = 1; i__ <= i__3; ++i__) {
s[i__] = 1.;
/* L40: */
}
} else {
/* (Scaled) random matrix */
if (jtype == 3) {
anorm = 1.;
}
if (jtype == 4) {
anorm = unfl / ulp;
}
if (jtype == 5) {
anorm = ovfl * ulp;
}
d__1 = (doublereal) mnmin;
i__3 = m - 1;
i__4 = n - 1;
zlatms_(&m, &n, "U", &iseed[1], "N", &s[1], &c__4, &d__1, &
anorm, &i__3, &i__4, "N", &a[a_offset], lda, &work[1],
&iinfo);
if (iinfo != 0) {
io___27.ciunit = *nounit;
s_wsfe(&io___27);
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;
}
}
L50:
zlacpy_("F", &m, &n, &a[a_offset], lda, &asav[asav_offset], lda);
/* Do for minimal and adequate (for blocking) workspace */
for (iwspc = 1; iwspc <= 4; ++iwspc) {
/* Test for ZGESVD */
iwtmp = (min(m,n) << 1) + max(m,n);
lswork = iwtmp + (iwspc - 1) * (*lwork - iwtmp) / 3;
lswork = min(lswork,*lwork);
lswork = max(lswork,1);
if (iwspc == 4) {
lswork = *lwork;
}
for (j = 1; j <= 14; ++j) {
result[j - 1] = -1.;
/* L60: */
}
/* Factorize A */
if (iwspc > 1) {
zlacpy_("F", &m, &n, &asav[asav_offset], lda, &a[a_offset]
, lda);
}
zgesvd_("A", "A", &m, &n, &a[a_offset], lda, &ssav[1], &usav[
usav_offset], ldu, &vtsav[vtsav_offset], ldvt, &work[
1], &lswork, &rwork[1], &iinfo);
if (iinfo != 0) {
io___32.ciunit = *nounit;
s_wsfe(&io___32);
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 */
zbdt01_(&m, &n, &c__0, &asav[asav_offset], lda, &usav[
usav_offset], ldu, &ssav[1], &e[1], &vtsav[
vtsav_offset], ldvt, &work[1], &rwork[1], result);
if (m != 0 && n != 0) {
zunt01_("Columns", &mnmin, &m, &usav[usav_offset], ldu, &
work[1], lwork, &rwork[1], &result[1]);
zunt01_("Rows", &mnmin, &n, &vtsav[vtsav_offset], ldvt, &
work[1], lwork, &rwork[1], &result[2]);
}
result[3] = 0.;
i__3 = mnmin - 1;
for (i__ = 1; i__ <= i__3; ++i__) {
if (ssav[i__] < ssav[i__ + 1]) {
result[3] = ulpinv;
}
if (ssav[i__] < 0.) {
result[3] = ulpinv;
}
/* L70: */
}
if (mnmin >= 1) {
if (ssav[mnmin] < 0.) {
result[3] = ulpinv;
}
}
/* Do partial SVDs, comparing to SSAV, USAV, and VTSAV */
result[4] = 0.;
result[5] = 0.;
result[6] = 0.;
for (iju = 0; iju <= 3; ++iju) {
for (ijvt = 0; ijvt <= 3; ++ijvt) {
if (iju == 3 && ijvt == 3 || iju == 1 && ijvt == 1) {
goto L90;
}
*(unsigned char *)jobu = *(unsigned char *)&cjob[iju];
*(unsigned char *)jobvt = *(unsigned char *)&cjob[
ijvt];
zlacpy_("F", &m, &n, &asav[asav_offset], lda, &a[
a_offset], lda);
zgesvd_(jobu, jobvt, &m, &n, &a[a_offset], lda, &s[1],
&u[u_offset], ldu, &vt[vt_offset], ldvt, &
work[1], &lswork, &rwork[1], &iinfo);
/* Compare U */
dif = 0.;
if (m > 0 && n > 0) {
if (iju == 1) {
zunt03_("C", &m, &mnmin, &m, &mnmin, &usav[
usav_offset], ldu, &a[a_offset], lda,
&work[1], lwork, &rwork[1], &dif, &
iinfo);
} else if (iju == 2) {
zunt03_("C", &m, &mnmin, &m, &mnmin, &usav[
usav_offset], ldu, &u[u_offset], ldu,
&work[1], lwork, &rwork[1], &dif, &
iinfo);
} else if (iju == 3) {
zunt03_("C", &m, &m, &m, &mnmin, &usav[
usav_offset], ldu, &u[u_offset], ldu,
&work[1], lwork, &rwork[1], &dif, &
iinfo);
}
}
result[4] = max(result[4],dif);
/* Compare VT */
dif = 0.;
if (m > 0 && n > 0) {
if (ijvt == 1) {
zunt03_("R", &n, &mnmin, &n, &mnmin, &vtsav[
vtsav_offset], ldvt, &a[a_offset],
lda, &work[1], lwork, &rwork[1], &dif,
&iinfo);
} else if (ijvt == 2) {
zunt03_("R", &n, &mnmin, &n, &mnmin, &vtsav[
vtsav_offset], ldvt, &vt[vt_offset],
ldvt, &work[1], lwork, &rwork[1], &
dif, &iinfo);
} else if (ijvt == 3) {
zunt03_("R", &n, &n, &n, &mnmin, &vtsav[
vtsav_offset], ldvt, &vt[vt_offset],
ldvt, &work[1], lwork, &rwork[1], &
dif, &iinfo);
}
}
result[5] = max(result[5],dif);
/* Compare S */
dif = 0.;
/* Computing MAX */
d__1 = (doublereal) mnmin * ulp * s[1], d__2 =
dlamch_("Safe minimum");
div = max(d__1,d__2);
i__3 = mnmin - 1;
for (i__ = 1; i__ <= i__3; ++i__) {
if (ssav[i__] < ssav[i__ + 1]) {
dif = ulpinv;
}
if (ssav[i__] < 0.) {
dif = ulpinv;
}
/* Computing MAX */
d__2 = dif, d__3 = (d__1 = ssav[i__] - s[i__],
abs(d__1)) / div;
dif = max(d__2,d__3);
/* L80: */
}
result[6] = max(result[6],dif);
L90:
;
}
/* L100: */
}
/* Test for ZGESDD */
iwtmp = (mnmin << 1) * mnmin + (mnmin << 1) + max(m,n);
lswork = iwtmp + (iwspc - 1) * (*lwork - iwtmp) / 3;
lswork = min(lswork,*lwork);
lswork = max(lswork,1);
if (iwspc == 4) {
lswork = *lwork;
}
/* Factorize A */
zlacpy_("F", &m, &n, &asav[asav_offset], lda, &a[a_offset],
lda);
zgesdd_("A", &m, &n, &a[a_offset], lda, &ssav[1], &usav[
usav_offset], ldu, &vtsav[vtsav_offset], ldvt, &work[
1], &lswork, &rwork[1], &iwork[1], &iinfo);
if (iinfo != 0) {
io___39.ciunit = *nounit;
s_wsfe(&io___39);
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 1--4 */
zbdt01_(&m, &n, &c__0, &asav[asav_offset], lda, &usav[
usav_offset], ldu, &ssav[1], &e[1], &vtsav[
vtsav_offset], ldvt, &work[1], &rwork[1], &result[7]);
if (m != 0 && n != 0) {
zunt01_("Columns", &mnmin, &m, &usav[usav_offset], ldu, &
work[1], lwork, &rwork[1], &result[8]);
zunt01_("Rows", &mnmin, &n, &vtsav[vtsav_offset], ldvt, &
work[1], lwork, &rwork[1], &result[9]);
}
result[10] = 0.;
i__3 = mnmin - 1;
for (i__ = 1; i__ <= i__3; ++i__) {
if (ssav[i__] < ssav[i__ + 1]) {
result[10] = ulpinv;
}
if (ssav[i__] < 0.) {
result[10] = ulpinv;
}
/* L110: */
}
if (mnmin >= 1) {
if (ssav[mnmin] < 0.) {
result[10] = ulpinv;
}
}
/* Do partial SVDs, comparing to SSAV, USAV, and VTSAV */
result[11] = 0.;
result[12] = 0.;
result[13] = 0.;
for (ijq = 0; ijq <= 2; ++ijq) {
*(unsigned char *)jobq = *(unsigned char *)&cjob[ijq];
zlacpy_("F", &m, &n, &asav[asav_offset], lda, &a[a_offset]
, lda);
zgesdd_(jobq, &m, &n, &a[a_offset], lda, &s[1], &u[
u_offset], ldu, &vt[vt_offset], ldvt, &work[1], &
lswork, &rwork[1], &iwork[1], &iinfo);
/* Compare U */
dif = 0.;
if (m > 0 && n > 0) {
if (ijq == 1) {
if (m >= n) {
zunt03_("C", &m, &mnmin, &m, &mnmin, &usav[
usav_offset], ldu, &a[a_offset], lda,
&work[1], lwork, &rwork[1], &dif, &
iinfo);
} else {
zunt03_("C", &m, &mnmin, &m, &mnmin, &usav[
usav_offset], ldu, &u[u_offset], ldu,
&work[1], lwork, &rwork[1], &dif, &
iinfo);
}
} else if (ijq == 2) {
zunt03_("C", &m, &mnmin, &m, &mnmin, &usav[
usav_offset], ldu, &u[u_offset], ldu, &
work[1], lwork, &rwork[1], &dif, &iinfo);
}
}
result[11] = max(result[11],dif);
/* Compare VT */
dif = 0.;
if (m > 0 && n > 0) {
if (ijq == 1) {
if (m >= n) {
zunt03_("R", &n, &mnmin, &n, &mnmin, &vtsav[
vtsav_offset], ldvt, &vt[vt_offset],
ldvt, &work[1], lwork, &rwork[1], &
dif, &iinfo);
} else {
zunt03_("R", &n, &mnmin, &n, &mnmin, &vtsav[
vtsav_offset], ldvt, &a[a_offset],
lda, &work[1], lwork, &rwork[1], &dif,
&iinfo);
}
} else if (ijq == 2) {
zunt03_("R", &n, &mnmin, &n, &mnmin, &vtsav[
vtsav_offset], ldvt, &vt[vt_offset], ldvt,
&work[1], lwork, &rwork[1], &dif, &iinfo);
}
}
result[12] = max(result[12],dif);
/* Compare S */
dif = 0.;
/* Computing MAX */
d__1 = (doublereal) mnmin * ulp * s[1], d__2 = dlamch_(
"Safe minimum");
div = max(d__1,d__2);
i__3 = mnmin - 1;
for (i__ = 1; i__ <= i__3; ++i__) {
if (ssav[i__] < ssav[i__ + 1]) {
dif = ulpinv;
}
if (ssav[i__] < 0.) {
dif = ulpinv;
}
/* Computing MAX */
d__2 = dif, d__3 = (d__1 = ssav[i__] - s[i__], abs(
d__1)) / div;
dif = max(d__2,d__3);
/* L120: */
}
result[13] = max(result[13],dif);
/* L130: */
}
/* End of Loop -- Check for RESULT(j) > THRESH */
ntest = 0;
nfail = 0;
for (j = 1; j <= 14; ++j) {
if (result[j - 1] >= 0.) {
++ntest;
}
if (result[j - 1] >= *thresh) {
++nfail;
}
/* L140: */
}
if (nfail > 0) {
++ntestf;
}
if (ntestf == 1) {
io___43.ciunit = *nounit;
s_wsfe(&io___43);
e_wsfe();
io___44.ciunit = *nounit;
s_wsfe(&io___44);
do_fio(&c__1, (char *)&(*thresh), (ftnlen)sizeof(
doublereal));
e_wsfe();
ntestf = 2;
}
for (j = 1; j <= 14; ++j) {
if (result[j - 1] >= *thresh) {
io___45.ciunit = *nounit;
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 *)&iwspc, (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(
doublereal));
e_wsfe();
}
/* L150: */
}
nerrs += nfail;
ntestt += ntest;
/* L160: */
}
L170:
;
}
/* L180: */
}
/* Summary */
alasvm_("ZBD", nounit, &nerrs, &ntestt, &c__0);
return 0;
/* End of ZDRVBD */
} /* zdrvbd_ */