/* cchkbd.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 complex c_b1 = {0.f,0.f}; static complex c_b2 = {1.f,0.f}; static integer c__0 = 0; static integer c__6 = 6; static real c_b37 = 1.f; static integer c__1 = 1; static real c_b47 = 0.f; static integer c__2 = 2; static integer c__4 = 4; /* Subroutine */ int cchkbd_(integer *nsizes, integer *mval, integer *nval, integer *ntypes, logical *dotype, integer *nrhs, integer *iseed, real *thresh, complex *a, integer *lda, real *bd, real *be, real *s1, real *s2, complex *x, integer *ldx, complex *y, complex *z__, complex *q, integer *ldq, complex *pt, integer *ldpt, complex *u, complex *vt, complex *work, integer *lwork, real *rwork, integer *nout, integer * info) { /* Initialized data */ static integer ktype[16] = { 1,2,4,4,4,4,4,6,6,6,6,6,9,9,9,10 }; static integer kmagn[16] = { 1,1,1,1,1,2,3,1,1,1,2,3,1,2,3,0 }; static integer kmode[16] = { 0,0,4,3,1,4,4,4,3,1,4,4,0,0,0,0 }; /* Format strings */ static char fmt_9998[] = "(\002 CCHKBD: \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_9999[] = "(\002 M=\002,i5,\002, N=\002,i5,\002, type " "\002,i2,\002, seed=\002,4(i4,\002,\002),\002 test(\002,i2,\002)" "=\002,g11.4)"; /* System generated locals */ integer a_dim1, a_offset, pt_dim1, pt_offset, q_dim1, q_offset, u_dim1, u_offset, vt_dim1, vt_offset, x_dim1, x_offset, y_dim1, y_offset, z_dim1, z_offset, i__1, i__2, i__3, i__4, i__5, i__6, i__7; real r__1, r__2, r__3, r__4, r__5, r__6, r__7; /* Builtin functions */ /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen); double log(doublereal), sqrt(doublereal), exp(doublereal); integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void); /* Local variables */ integer i__, j, m, n, mq; real ulp, cond; integer jcol; char path[3]; integer mmax, nmax; real unfl, ovfl; char uplo[1]; real temp1, temp2; extern /* Subroutine */ int cbdt01_(integer *, integer *, integer *, complex *, integer *, complex *, integer *, real *, real *, complex *, integer *, complex *, real *, real *), cbdt02_(integer *, integer *, complex *, integer *, complex *, integer *, complex *, integer *, complex *, real *, real *), cbdt03_(char *, integer *, integer *, real *, real *, complex *, integer *, real *, complex *, integer *, complex *, real *); logical badmm, badnn; extern /* Subroutine */ int cgemm_(char *, char *, integer *, integer *, integer *, complex *, complex *, integer *, complex *, integer *, complex *, complex *, integer *); integer nfail, imode; real dumma[1]; integer iinfo; extern /* Subroutine */ int cunt01_(char *, integer *, integer *, complex *, integer *, complex *, integer *, real *, real *); real anorm; integer mnmin, mnmax, jsize, itype, jtype, iwork[1], ntest; extern /* Subroutine */ int scopy_(integer *, real *, integer *, real *, integer *), slahd2_(integer *, char *); integer log2ui; logical bidiag; extern /* Subroutine */ int cgebrd_(integer *, integer *, complex *, integer *, real *, real *, complex *, complex *, complex *, integer *, integer *), slabad_(real *, real *); extern doublereal slamch_(char *); extern /* Subroutine */ int clacpy_(char *, integer *, integer *, complex *, integer *, complex *, integer *), claset_(char *, integer *, integer *, complex *, complex *, complex *, integer *), xerbla_(char *, integer *); integer ioldsd[4]; extern /* Subroutine */ int cbdsqr_(char *, integer *, integer *, integer *, integer *, real *, real *, complex *, integer *, complex *, integer *, complex *, integer *, real *, integer *), cungbr_(char *, integer *, integer *, integer *, complex *, integer *, complex *, complex *, integer *, integer *), alasum_(char *, integer *, integer *, integer *, integer *); extern doublereal slarnd_(integer *, integer *); extern /* Subroutine */ int clatmr_(integer *, integer *, char *, integer *, char *, complex *, integer *, real *, complex *, char *, char * , complex *, integer *, real *, complex *, integer *, real *, char *, integer *, integer *, integer *, real *, real *, char *, complex *, integer *, integer *, integer *), clatms_(integer *, integer *, char *, integer *, char *, real *, integer *, real *, real *, integer *, integer *, char *, complex *, integer *, complex *, integer *); real amninv; extern /* Subroutine */ int ssvdch_(integer *, real *, real *, real *, real *, integer *); integer minwrk; real rtunfl, rtovfl, ulpinv, result[14]; integer mtypes; /* Fortran I/O blocks */ static cilist io___40 = { 0, 0, 0, fmt_9998, 0 }; static cilist io___41 = { 0, 0, 0, fmt_9998, 0 }; static cilist io___43 = { 0, 0, 0, fmt_9998, 0 }; static cilist io___44 = { 0, 0, 0, fmt_9998, 0 }; static cilist io___45 = { 0, 0, 0, fmt_9998, 0 }; static cilist io___46 = { 0, 0, 0, fmt_9998, 0 }; static cilist io___50 = { 0, 0, 0, fmt_9999, 0 }; /* -- LAPACK test routine (version 3.1) -- */ /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ /* November 2006 */ /* .. Scalar Arguments .. */ /* .. */ /* .. Array Arguments .. */ /* .. */ /* Purpose */ /* ======= */ /* CCHKBD checks the singular value decomposition (SVD) routines. */ /* CGEBRD reduces a complex general m by n matrix A to real upper or */ /* lower bidiagonal form 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. */ /* CUNGBR generates the orthogonal matrices Q and P' from CGEBRD. */ /* Note that Q and P are not necessarily square. */ /* CBDSQR 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, CBDSQR has an option to apply the left orthogonal matrix */ /* U to a matrix X, useful in least squares applications. */ /* 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 CGEBRD and CUNGBR */ /* (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 CBDSQR 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) 0 if the true singular values of B are within THRESH of */ /* those in S1. 2*THRESH if they are not. (Tested using */ /* SSVDCH) */ /* (10) | S1 - S2 | / ( |S1| ulp ), where S2 is computed without */ /* computing U and V. */ /* Test CBDSQR 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 ) */ /* 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) CGEBRD is not called to reduce it to bidiagonal form. */ /* (b) the bidiagonal is min(M,N) x min(M,N); if M= 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. */ /* A (workspace) COMPLEX array, dimension (LDA,NMAX) */ /* where NMAX is the maximum value of N in NVAL. */ /* LDA (input) INTEGER */ /* The leading dimension of the array A. LDA >= max(1,MMAX), */ /* where MMAX is the maximum value of M in MVAL. */ /* BD (workspace) REAL array, dimension */ /* (max(min(MVAL(j),NVAL(j)))) */ /* BE (workspace) REAL array, dimension */ /* (max(min(MVAL(j),NVAL(j)))) */ /* S1 (workspace) REAL array, dimension */ /* (max(min(MVAL(j),NVAL(j)))) */ /* S2 (workspace) REAL array, dimension */ /* (max(min(MVAL(j),NVAL(j)))) */ /* X (workspace) COMPLEX array, dimension (LDX,NRHS) */ /* LDX (input) INTEGER */ /* The leading dimension of the arrays X, Y, and Z. */ /* LDX >= max(1,MMAX). */ /* Y (workspace) COMPLEX array, dimension (LDX,NRHS) */ /* Z (workspace) COMPLEX array, dimension (LDX,NRHS) */ /* Q (workspace) COMPLEX array, dimension (LDQ,MMAX) */ /* LDQ (input) INTEGER */ /* The leading dimension of the array Q. LDQ >= max(1,MMAX). */ /* PT (workspace) COMPLEX array, dimension (LDPT,NMAX) */ /* LDPT (input) INTEGER */ /* The leading dimension of the arrays PT, U, and V. */ /* LDPT >= max(1, max(min(MVAL(j),NVAL(j)))). */ /* U (workspace) COMPLEX array, dimension */ /* (LDPT,max(min(MVAL(j),NVAL(j)))) */ /* V (workspace) COMPLEX array, dimension */ /* (LDPT,max(min(MVAL(j),NVAL(j)))) */ /* WORK (workspace) COMPLEX array, dimension (LWORK) */ /* LWORK (input) 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)) */ /* RWORK (workspace) REAL array, dimension */ /* (5*max(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 */ /* -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: LDP < 1 or LDP < MNMAX. */ /* -27: LWORK too small. */ /* If CLATMR, CLATMS, CGEBRD, CUNGBR, or CBDSQR, */ /* 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) ) */ /* ====================================================================== */ /* .. Parameters .. */ /* .. */ /* .. Local Scalars .. */ /* .. */ /* .. Local Arrays .. */ /* .. */ /* .. External Functions .. */ /* .. */ /* .. External Subroutines .. */ /* .. */ /* .. Intrinsic Functions .. */ /* .. */ /* .. Scalars in Common .. */ /* .. */ /* .. Common blocks .. */ /* .. */ /* .. Data statements .. */ /* Parameter adjustments */ --mval; --nval; --dotype; --iseed; a_dim1 = *lda; a_offset = 1 + a_dim1; a -= a_offset; --bd; --be; --s1; --s2; z_dim1 = *ldx; z_offset = 1 + z_dim1; z__ -= z_offset; y_dim1 = *ldx; y_offset = 1 + y_dim1; y -= y_offset; x_dim1 = *ldx; x_offset = 1 + x_dim1; x -= x_offset; q_dim1 = *ldq; q_offset = 1 + q_dim1; q -= q_offset; vt_dim1 = *ldpt; vt_offset = 1 + vt_dim1; vt -= vt_offset; u_dim1 = *ldpt; u_offset = 1 + u_dim1; u -= u_offset; pt_dim1 = *ldpt; pt_offset = 1 + pt_dim1; pt -= pt_offset; --work; --rwork; /* 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 = mval[j]; mmax = max(i__2,i__3); if (mval[j] < 0) { badmm = TRUE_; } /* Computing MAX */ i__2 = nmax, i__3 = nval[j]; nmax = max(i__2,i__3); if (nval[j] < 0) { badnn = TRUE_; } /* Computing MAX */ /* Computing MIN */ i__4 = mval[j], i__5 = nval[j]; i__2 = mnmax, i__3 = min(i__4,i__5); mnmax = max(i__2,i__3); /* Computing MAX */ /* Computing MAX */ i__4 = mval[j], i__5 = nval[j], i__4 = max(i__4,i__5); /* Computing MIN */ i__6 = nval[j], i__7 = mval[j]; i__2 = minwrk, i__3 = (mval[j] + nval[j]) * 3, i__2 = max(i__2,i__3), i__3 = mval[j] * (mval[j] + max(i__4,*nrhs) + 1) + nval[j] * min(i__6,i__7); 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 (*nrhs < 0) { *info = -6; } else if (*lda < mmax) { *info = -11; } else if (*ldx < mmax) { *info = -17; } else if (*ldq < mmax) { *info = -21; } else if (*ldpt < mnmax) { *info = -23; } else if (minwrk > *lwork) { *info = -27; } if (*info != 0) { i__1 = -(*info); xerbla_("CCHKBD", &i__1); return 0; } /* Initialize constants */ s_copy(path, "Complex precision", (ftnlen)1, (ftnlen)17); s_copy(path + 1, "BD", (ftnlen)2, (ftnlen)2); nfail = 0; ntest = 0; unfl = slamch_("Safe minimum"); ovfl = slamch_("Overflow"); slabad_(&unfl, &ovfl); ulp = slamch_("Precision"); ulpinv = 1.f / ulp; log2ui = (integer) (log(ulpinv) / log(2.f)); rtunfl = sqrt(unfl); rtovfl = sqrt(ovfl); infoc_1.infot = 0; /* Loop over sizes, types */ i__1 = *nsizes; for (jsize = 1; jsize <= i__1; ++jsize) { m = mval[jsize]; n = nval[jsize]; mnmin = min(m,n); /* Computing MAX */ i__2 = max(m,n); amninv = 1.f / max(i__2,1); if (*nsizes != 1) { mtypes = min(16,*ntypes); } else { mtypes = min(17,*ntypes); } i__2 = mtypes; for (jtype = 1; jtype <= i__2; ++jtype) { if (! dotype[jtype]) { goto L170; } for (j = 1; j <= 4; ++j) { ioldsd[j - 1] = iseed[j]; /* L20: */ } for (j = 1; j <= 14; ++j) { result[j - 1] = -1.f; /* L30: */ } *(unsigned char *)uplo = ' '; /* Compute "A" */ /* Control parameters: */ /* KMAGN KMODE KTYPE */ /* =1 O(1) clustered 1 zero */ /* =2 large clustered 2 identity */ /* =3 small exponential (none) */ /* =4 arithmetic diagonal, (w/ eigenvalues) */ /* =5 random symmetric, w/ eigenvalues */ /* =6 nonsymmetric, w/ singular values */ /* =7 random diagonal */ /* =8 random symmetric */ /* =9 random nonsymmetric */ /* =10 random bidiagonal (log. distrib.) */ if (mtypes > 16) { goto L100; } itype = ktype[jtype - 1]; imode = kmode[jtype - 1]; /* Compute norm */ switch (kmagn[jtype - 1]) { case 1: goto L40; case 2: goto L50; case 3: goto L60; } L40: anorm = 1.f; goto L70; L50: anorm = rtovfl * ulp * amninv; goto L70; L60: anorm = rtunfl * max(m,n) * ulpinv; goto L70; L70: claset_("Full", lda, &n, &c_b1, &c_b1, &a[a_offset], lda); iinfo = 0; cond = ulpinv; bidiag = FALSE_; if (itype == 1) { /* Zero matrix */ iinfo = 0; } else if (itype == 2) { /* Identity */ i__3 = mnmin; for (jcol = 1; jcol <= i__3; ++jcol) { i__4 = jcol + jcol * a_dim1; a[i__4].r = anorm, a[i__4].i = 0.f; /* L80: */ } } else if (itype == 4) { /* Diagonal Matrix, [Eigen]values Specified */ clatms_(&mnmin, &mnmin, "S", &iseed[1], "N", &rwork[1], & imode, &cond, &anorm, &c__0, &c__0, "N", &a[a_offset], lda, &work[1], &iinfo); } else if (itype == 5) { /* Symmetric, eigenvalues specified */ clatms_(&mnmin, &mnmin, "S", &iseed[1], "S", &rwork[1], & imode, &cond, &anorm, &m, &n, "N", &a[a_offset], lda, &work[1], &iinfo); } else if (itype == 6) { /* Nonsymmetric, singular values specified */ clatms_(&m, &n, "S", &iseed[1], "N", &rwork[1], &imode, &cond, &anorm, &m, &n, "N", &a[a_offset], lda, &work[1], & iinfo); } else if (itype == 7) { /* Diagonal, random entries */ clatmr_(&mnmin, &mnmin, "S", &iseed[1], "N", &work[1], &c__6, &c_b37, &c_b2, "T", "N", &work[mnmin + 1], &c__1, & c_b37, &work[(mnmin << 1) + 1], &c__1, &c_b37, "N", iwork, &c__0, &c__0, &c_b47, &anorm, "NO", &a[ a_offset], lda, iwork, &iinfo); } else if (itype == 8) { /* Symmetric, random entries */ clatmr_(&mnmin, &mnmin, "S", &iseed[1], "S", &work[1], &c__6, &c_b37, &c_b2, "T", "N", &work[mnmin + 1], &c__1, & c_b37, &work[m + mnmin + 1], &c__1, &c_b37, "N", iwork, &m, &n, &c_b47, &anorm, "NO", &a[a_offset], lda, iwork, &iinfo); } else if (itype == 9) { /* Nonsymmetric, random entries */ clatmr_(&m, &n, "S", &iseed[1], "N", &work[1], &c__6, &c_b37, &c_b2, "T", "N", &work[mnmin + 1], &c__1, &c_b37, & work[m + mnmin + 1], &c__1, &c_b37, "N", iwork, &m, & n, &c_b47, &anorm, "NO", &a[a_offset], lda, iwork, & iinfo); } else if (itype == 10) { /* Bidiagonal, random entries */ temp1 = log(ulp) * -2.f; i__3 = mnmin; for (j = 1; j <= i__3; ++j) { bd[j] = exp(temp1 * slarnd_(&c__2, &iseed[1])); if (j < mnmin) { be[j] = exp(temp1 * slarnd_(&c__2, &iseed[1])); } /* L90: */ } iinfo = 0; bidiag = TRUE_; if (m >= n) { *(unsigned char *)uplo = 'U'; } else { *(unsigned char *)uplo = 'L'; } } else { iinfo = 1; } if (iinfo == 0) { /* Generate Right-Hand Side */ if (bidiag) { clatmr_(&mnmin, nrhs, "S", &iseed[1], "N", &work[1], & c__6, &c_b37, &c_b2, "T", "N", &work[mnmin + 1], & c__1, &c_b37, &work[(mnmin << 1) + 1], &c__1, & c_b37, "N", iwork, &mnmin, nrhs, &c_b47, &c_b37, "NO", &y[y_offset], ldx, iwork, &iinfo); } else { clatmr_(&m, nrhs, "S", &iseed[1], "N", &work[1], &c__6, & c_b37, &c_b2, "T", "N", &work[m + 1], &c__1, & c_b37, &work[(m << 1) + 1], &c__1, &c_b37, "N", iwork, &m, nrhs, &c_b47, &c_b37, "NO", &x[ x_offset], ldx, iwork, &iinfo); } } /* Error Exit */ if (iinfo != 0) { io___40.ciunit = *nout; s_wsfe(&io___40); 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; } L100: /* Call CGEBRD and CUNGBR to compute B, Q, and P, do tests. */ if (! bidiag) { /* Compute transformations to reduce A to bidiagonal form: */ /* B := Q' * A * P. */ clacpy_(" ", &m, &n, &a[a_offset], lda, &q[q_offset], ldq); i__3 = *lwork - (mnmin << 1); cgebrd_(&m, &n, &q[q_offset], ldq, &bd[1], &be[1], &work[1], & work[mnmin + 1], &work[(mnmin << 1) + 1], &i__3, & iinfo); /* Check error code from CGEBRD. */ if (iinfo != 0) { io___41.ciunit = *nout; s_wsfe(&io___41); do_fio(&c__1, "CGEBRD", (ftnlen)6); 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; } clacpy_(" ", &m, &n, &q[q_offset], ldq, &pt[pt_offset], ldpt); if (m >= n) { *(unsigned char *)uplo = 'U'; } else { *(unsigned char *)uplo = 'L'; } /* Generate Q */ mq = m; if (*nrhs <= 0) { mq = mnmin; } i__3 = *lwork - (mnmin << 1); cungbr_("Q", &m, &mq, &n, &q[q_offset], ldq, &work[1], &work[( mnmin << 1) + 1], &i__3, &iinfo); /* Check error code from CUNGBR. */ if (iinfo != 0) { io___43.ciunit = *nout; s_wsfe(&io___43); do_fio(&c__1, "CUNGBR(Q)", (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; } /* Generate P' */ i__3 = *lwork - (mnmin << 1); cungbr_("P", &mnmin, &n, &m, &pt[pt_offset], ldpt, &work[ mnmin + 1], &work[(mnmin << 1) + 1], &i__3, &iinfo); /* Check error code from CUNGBR. */ if (iinfo != 0) { io___44.ciunit = *nout; s_wsfe(&io___44); do_fio(&c__1, "CUNGBR(P)", (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; } /* Apply Q' to an M by NRHS matrix X: Y := Q' * X. */ cgemm_("Conjugate transpose", "No transpose", &m, nrhs, &m, & c_b2, &q[q_offset], ldq, &x[x_offset], ldx, &c_b1, &y[ y_offset], ldx); /* Test 1: Check the decomposition A := Q * B * PT */ /* 2: Check the orthogonality of Q */ /* 3: Check the orthogonality of PT */ cbdt01_(&m, &n, &c__1, &a[a_offset], lda, &q[q_offset], ldq, & bd[1], &be[1], &pt[pt_offset], ldpt, &work[1], &rwork[ 1], result); cunt01_("Columns", &m, &mq, &q[q_offset], ldq, &work[1], lwork, &rwork[1], &result[1]); cunt01_("Rows", &mnmin, &n, &pt[pt_offset], ldpt, &work[1], lwork, &rwork[1], &result[2]); } /* Use CBDSQR to form the SVD of the bidiagonal matrix B: */ /* B := U * S1 * VT, and compute Z = U' * Y. */ scopy_(&mnmin, &bd[1], &c__1, &s1[1], &c__1); if (mnmin > 0) { i__3 = mnmin - 1; scopy_(&i__3, &be[1], &c__1, &rwork[1], &c__1); } clacpy_(" ", &m, nrhs, &y[y_offset], ldx, &z__[z_offset], ldx); claset_("Full", &mnmin, &mnmin, &c_b1, &c_b2, &u[u_offset], ldpt); claset_("Full", &mnmin, &mnmin, &c_b1, &c_b2, &vt[vt_offset], ldpt); cbdsqr_(uplo, &mnmin, &mnmin, &mnmin, nrhs, &s1[1], &rwork[1], & vt[vt_offset], ldpt, &u[u_offset], ldpt, &z__[z_offset], ldx, &rwork[mnmin + 1], &iinfo); /* Check error code from CBDSQR. */ if (iinfo != 0) { io___45.ciunit = *nout; s_wsfe(&io___45); do_fio(&c__1, "CBDSQR(vects)", (ftnlen)13); 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); if (iinfo < 0) { return 0; } else { result[3] = ulpinv; goto L150; } } /* Use CBDSQR to compute only the singular values of the */ /* bidiagonal matrix B; U, VT, and Z should not be modified. */ scopy_(&mnmin, &bd[1], &c__1, &s2[1], &c__1); if (mnmin > 0) { i__3 = mnmin - 1; scopy_(&i__3, &be[1], &c__1, &rwork[1], &c__1); } cbdsqr_(uplo, &mnmin, &c__0, &c__0, &c__0, &s2[1], &rwork[1], &vt[ vt_offset], ldpt, &u[u_offset], ldpt, &z__[z_offset], ldx, &rwork[mnmin + 1], &iinfo); /* Check error code from CBDSQR. */ if (iinfo != 0) { io___46.ciunit = *nout; s_wsfe(&io___46); do_fio(&c__1, "CBDSQR(values)", (ftnlen)14); 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); if (iinfo < 0) { return 0; } else { result[8] = ulpinv; goto L150; } } /* Test 4: Check the decomposition B := U * S1 * VT */ /* 5: Check the computation Z := U' * Y */ /* 6: Check the orthogonality of U */ /* 7: Check the orthogonality of VT */ cbdt03_(uplo, &mnmin, &c__1, &bd[1], &be[1], &u[u_offset], ldpt, & s1[1], &vt[vt_offset], ldpt, &work[1], &result[3]); cbdt02_(&mnmin, nrhs, &y[y_offset], ldx, &z__[z_offset], ldx, &u[ u_offset], ldpt, &work[1], &rwork[1], &result[4]); cunt01_("Columns", &mnmin, &mnmin, &u[u_offset], ldpt, &work[1], lwork, &rwork[1], &result[5]); cunt01_("Rows", &mnmin, &mnmin, &vt[vt_offset], ldpt, &work[1], lwork, &rwork[1], &result[6]); /* Test 8: Check that the singular values are sorted in */ /* non-increasing order and are non-negative */ result[7] = 0.f; i__3 = mnmin - 1; for (i__ = 1; i__ <= i__3; ++i__) { if (s1[i__] < s1[i__ + 1]) { result[7] = ulpinv; } if (s1[i__] < 0.f) { result[7] = ulpinv; } /* L110: */ } if (mnmin >= 1) { if (s1[mnmin] < 0.f) { result[7] = ulpinv; } } /* Test 9: Compare CBDSQR with and without singular vectors */ temp2 = 0.f; i__3 = mnmin; for (j = 1; j <= i__3; ++j) { /* Computing MAX */ /* Computing MAX */ r__6 = (r__1 = s1[j], dabs(r__1)), r__7 = (r__2 = s2[j], dabs( r__2)); r__4 = sqrt(unfl) * dmax(s1[1],1.f), r__5 = ulp * dmax(r__6, r__7); temp1 = (r__3 = s1[j] - s2[j], dabs(r__3)) / dmax(r__4,r__5); temp2 = dmax(temp1,temp2); /* L120: */ } result[8] = temp2; /* Test 10: Sturm sequence test of singular values */ /* Go up by factors of two until it succeeds */ temp1 = *thresh * (.5f - ulp); i__3 = log2ui; for (j = 0; j <= i__3; ++j) { ssvdch_(&mnmin, &bd[1], &be[1], &s1[1], &temp1, &iinfo); if (iinfo == 0) { goto L140; } temp1 *= 2.f; /* L130: */ } L140: result[9] = temp1; /* Use CBDSQR to form the decomposition A := (QU) S (VT PT) */ /* from the bidiagonal form A := Q B PT. */ if (! bidiag) { scopy_(&mnmin, &bd[1], &c__1, &s2[1], &c__1); if (mnmin > 0) { i__3 = mnmin - 1; scopy_(&i__3, &be[1], &c__1, &rwork[1], &c__1); } cbdsqr_(uplo, &mnmin, &n, &m, nrhs, &s2[1], &rwork[1], &pt[ pt_offset], ldpt, &q[q_offset], ldq, &y[y_offset], ldx, &rwork[mnmin + 1], &iinfo); /* Test 11: Check the decomposition A := Q*U * S2 * VT*PT */ /* 12: Check the computation Z := U' * Q' * X */ /* 13: Check the orthogonality of Q*U */ /* 14: Check the orthogonality of VT*PT */ cbdt01_(&m, &n, &c__0, &a[a_offset], lda, &q[q_offset], ldq, & s2[1], dumma, &pt[pt_offset], ldpt, &work[1], &rwork[ 1], &result[10]); cbdt02_(&m, nrhs, &x[x_offset], ldx, &y[y_offset], ldx, &q[ q_offset], ldq, &work[1], &rwork[1], &result[11]); cunt01_("Columns", &m, &mq, &q[q_offset], ldq, &work[1], lwork, &rwork[1], &result[12]); cunt01_("Rows", &mnmin, &n, &pt[pt_offset], ldpt, &work[1], lwork, &rwork[1], &result[13]); } /* End of Loop -- Check for RESULT(j) > THRESH */ L150: for (j = 1; j <= 14; ++j) { if (result[j - 1] >= *thresh) { if (nfail == 0) { slahd2_(nout, path); } io___50.ciunit = *nout; s_wsfe(&io___50); 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)) ; do_fio(&c__1, (char *)&j, (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&result[j - 1], (ftnlen)sizeof(real) ); e_wsfe(); ++nfail; } /* L160: */ } if (! bidiag) { ntest += 14; } else { ntest += 5; } L170: ; } /* L180: */ } /* Summary */ alasum_(path, nout, &nfail, &ntest, &c__0); return 0; /* End of CCHKBD */ } /* cchkbd_ */