cdrvls.c
Go to the documentation of this file.
00001 /* cdrvls.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Common Block Declarations */
00017 
00018 struct {
00019     integer infot, iounit;
00020     logical ok, lerr;
00021 } infoc_;
00022 
00023 #define infoc_1 infoc_
00024 
00025 struct {
00026     char srnamt[32];
00027 } srnamc_;
00028 
00029 #define srnamc_1 srnamc_
00030 
00031 /* Table of constant values */
00032 
00033 static complex c_b1 = {1.f,0.f};
00034 static complex c_b2 = {0.f,0.f};
00035 static integer c__9 = 9;
00036 static integer c__25 = 25;
00037 static integer c__1 = 1;
00038 static integer c__3 = 3;
00039 static integer c__2 = 2;
00040 static integer c__0 = 0;
00041 static integer c_n1 = -1;
00042 static real c_b91 = -1.f;
00043 
00044 /* Subroutine */ int cdrvls_(logical *dotype, integer *nm, integer *mval, 
00045         integer *nn, integer *nval, integer *nns, integer *nsval, integer *
00046         nnb, integer *nbval, integer *nxval, real *thresh, logical *tsterr, 
00047         complex *a, complex *copya, complex *b, complex *copyb, complex *c__, 
00048         real *s, real *copys, complex *work, real *rwork, integer *iwork, 
00049         integer *nout)
00050 {
00051     /* Initialized data */
00052 
00053     static integer iseedy[4] = { 1988,1989,1990,1991 };
00054 
00055     /* Format strings */
00056     static char fmt_9999[] = "(\002 TRANS='\002,a1,\002', M=\002,i5,\002, N"
00057             "=\002,i5,\002, NRHS=\002,i4,\002, NB=\002,i4,\002, type\002,i2"
00058             ",\002, test(\002,i2,\002)=\002,g12.5)";
00059     static char fmt_9998[] = "(\002 M=\002,i5,\002, N=\002,i5,\002, NRHS="
00060             "\002,i4,\002, NB=\002,i4,\002, type\002,i2,\002, test(\002,i2"
00061             ",\002)=\002,g12.5)";
00062 
00063     /* System generated locals */
00064     integer i__1, i__2, i__3, i__4, i__5, i__6;
00065     real r__1;
00066 
00067     /* Builtin functions */
00068     /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
00069     double sqrt(doublereal);
00070     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
00071 
00072     /* Local variables */
00073     integer i__, j, k, m, n, nb, im, in, lda, ldb, inb;
00074     real eps;
00075     integer ins, info;
00076     char path[3];
00077     integer rank, nrhs, nrun;
00078     extern /* Subroutine */ int alahd_(integer *, char *), cgemm_(
00079             char *, char *, integer *, integer *, integer *, complex *, 
00080             complex *, integer *, complex *, integer *, complex *, complex *, 
00081             integer *);
00082     integer nfail, iseed[4];
00083     extern /* Subroutine */ int cgels_(char *, integer *, integer *, integer *
00084 , complex *, integer *, complex *, integer *, complex *, integer *
00085 , integer *);
00086     integer crank, irank;
00087     real rcond;
00088     integer itran, mnmin, ncols;
00089     real norma, normb;
00090     extern doublereal cqrt12_(integer *, integer *, complex *, integer *, 
00091             real *, complex *, integer *, real *), cqrt14_(char *, integer *, 
00092             integer *, integer *, complex *, integer *, complex *, integer *, 
00093             complex *, integer *), cqrt17_(char *, integer *, integer 
00094             *, integer *, integer *, complex *, integer *, complex *, integer 
00095             *, complex *, integer *, complex *, complex *, integer *);
00096     char trans[1];
00097     integer nerrs, itype;
00098     extern doublereal sasum_(integer *, real *, integer *);
00099     integer lwork;
00100     extern /* Subroutine */ int cqrt13_(integer *, integer *, integer *, 
00101             complex *, integer *, real *, integer *), cqrt15_(integer *, 
00102             integer *, integer *, integer *, integer *, complex *, integer *, 
00103             complex *, integer *, real *, integer *, real *, real *, integer *
00104 , complex *, integer *), cqrt16_(char *, integer *, integer *, 
00105             integer *, complex *, integer *, complex *, integer *, complex *, 
00106             integer *, real *, real *), saxpy_(integer *, real *, 
00107             real *, integer *, real *, integer *);
00108     integer nrows, lwlsy;
00109     extern /* Subroutine */ int alaerh_(char *, char *, integer *, integer *, 
00110             char *, integer *, integer *, integer *, integer *, integer *, 
00111             integer *, integer *, integer *, integer *);
00112     integer iscale;
00113     extern /* Subroutine */ int cgelsd_(integer *, integer *, integer *, 
00114             complex *, integer *, complex *, integer *, real *, real *, 
00115             integer *, complex *, integer *, real *, integer *, integer *);
00116     extern doublereal slamch_(char *);
00117     extern /* Subroutine */ int csscal_(integer *, real *, complex *, integer 
00118             *), clacpy_(char *, integer *, integer *, complex *, integer *, 
00119             complex *, integer *), cgelss_(integer *, integer *, 
00120             integer *, complex *, integer *, complex *, integer *, real *, 
00121             real *, integer *, complex *, integer *, real *, integer *), 
00122             alasvm_(char *, integer *, integer *, integer *, integer *), cgelsx_(integer *, integer *, integer *, complex *, 
00123             integer *, complex *, integer *, integer *, real *, integer *, 
00124             complex *, real *, integer *), cgelsy_(integer *, integer *, 
00125             integer *, complex *, integer *, complex *, integer *, integer *, 
00126             real *, integer *, complex *, integer *, real *, integer *), 
00127             clarnv_(integer *, integer *, integer *, complex *), cerrls_(char 
00128             *, integer *), xlaenv_(integer *, integer *);
00129     integer ldwork;
00130     real result[18];
00131 
00132     /* Fortran I/O blocks */
00133     static cilist io___34 = { 0, 0, 0, fmt_9999, 0 };
00134     static cilist io___39 = { 0, 0, 0, fmt_9998, 0 };
00135     static cilist io___41 = { 0, 0, 0, fmt_9998, 0 };
00136 
00137 
00138 
00139 /*  -- LAPACK test routine (version 3.1.1) -- */
00140 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00141 /*     January 2007 */
00142 
00143 /*     .. Scalar Arguments .. */
00144 /*     .. */
00145 /*     .. Array Arguments .. */
00146 /*     .. */
00147 
00148 /*  Purpose */
00149 /*  ======= */
00150 
00151 /*  CDRVLS tests the least squares driver routines CGELS, CGELSX, CGELSS, */
00152 /*  CGELSY and CGELSD. */
00153 
00154 /*  Arguments */
00155 /*  ========= */
00156 
00157 /*  DOTYPE  (input) LOGICAL array, dimension (NTYPES) */
00158 /*          The matrix types to be used for testing.  Matrices of type j */
00159 /*          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = */
00160 /*          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. */
00161 /*          The matrix of type j is generated as follows: */
00162 /*          j=1: A = U*D*V where U and V are random unitary matrices */
00163 /*               and D has random entries (> 0.1) taken from a uniform */
00164 /*               distribution (0,1). A is full rank. */
00165 /*          j=2: The same of 1, but A is scaled up. */
00166 /*          j=3: The same of 1, but A is scaled down. */
00167 /*          j=4: A = U*D*V where U and V are random unitary matrices */
00168 /*               and D has 3*min(M,N)/4 random entries (> 0.1) taken */
00169 /*               from a uniform distribution (0,1) and the remaining */
00170 /*               entries set to 0. A is rank-deficient. */
00171 /*          j=5: The same of 4, but A is scaled up. */
00172 /*          j=6: The same of 5, but A is scaled down. */
00173 
00174 /*  NM      (input) INTEGER */
00175 /*          The number of values of M contained in the vector MVAL. */
00176 
00177 /*  MVAL    (input) INTEGER array, dimension (NM) */
00178 /*          The values of the matrix row dimension M. */
00179 
00180 /*  NN      (input) INTEGER */
00181 /*          The number of values of N contained in the vector NVAL. */
00182 
00183 /*  NVAL    (input) INTEGER array, dimension (NN) */
00184 /*          The values of the matrix column dimension N. */
00185 
00186 /*  NNB     (input) INTEGER */
00187 /*          The number of values of NB and NX contained in the */
00188 /*          vectors NBVAL and NXVAL.  The blocking parameters are used */
00189 /*          in pairs (NB,NX). */
00190 
00191 /*  NBVAL   (input) INTEGER array, dimension (NNB) */
00192 /*          The values of the blocksize NB. */
00193 
00194 /*  NXVAL   (input) INTEGER array, dimension (NNB) */
00195 /*          The values of the crossover point NX. */
00196 
00197 /*  NNS     (input) INTEGER */
00198 /*          The number of values of NRHS contained in the vector NSVAL. */
00199 
00200 /*  NSVAL   (input) INTEGER array, dimension (NNS) */
00201 /*          The values of the number of right hand sides NRHS. */
00202 
00203 /*  THRESH  (input) REAL */
00204 /*          The threshold value for the test ratios.  A result is */
00205 /*          included in the output file if RESULT >= THRESH.  To have */
00206 /*          every test ratio printed, use THRESH = 0. */
00207 
00208 /*  TSTERR  (input) LOGICAL */
00209 /*          Flag that indicates whether error exits are to be tested. */
00210 
00211 /*  A       (workspace) COMPLEX array, dimension (MMAX*NMAX) */
00212 /*          where MMAX is the maximum value of M in MVAL and NMAX is the */
00213 /*          maximum value of N in NVAL. */
00214 
00215 /*  COPYA   (workspace) COMPLEX array, dimension (MMAX*NMAX) */
00216 
00217 /*  B       (workspace) COMPLEX array, dimension (MMAX*NSMAX) */
00218 /*          where MMAX is the maximum value of M in MVAL and NSMAX is the */
00219 /*          maximum value of NRHS in NSVAL. */
00220 
00221 /*  COPYB   (workspace) COMPLEX array, dimension (MMAX*NSMAX) */
00222 
00223 /*  C       (workspace) COMPLEX array, dimension (MMAX*NSMAX) */
00224 
00225 /*  S       (workspace) REAL array, dimension */
00226 /*                      (min(MMAX,NMAX)) */
00227 
00228 /*  COPYS   (workspace) REAL array, dimension */
00229 /*                      (min(MMAX,NMAX)) */
00230 
00231 /*  WORK    (workspace) COMPLEX array, dimension */
00232 /*                      (MMAX*NMAX + 4*NMAX + MMAX). */
00233 
00234 /*  RWORK   (workspace) REAL array, dimension (5*NMAX-1) */
00235 
00236 /*  IWORK   (workspace) INTEGER array, dimension (15*NMAX) */
00237 
00238 /*  NOUT    (input) INTEGER */
00239 /*          The unit number for output. */
00240 
00241 /*  ===================================================================== */
00242 
00243 /*     .. Parameters .. */
00244 /*     .. */
00245 /*     .. Local Scalars .. */
00246 /*     .. */
00247 /*     .. Local Arrays .. */
00248 /*     .. */
00249 /*     .. External Functions .. */
00250 /*     .. */
00251 /*     .. External Subroutines .. */
00252 /*     .. */
00253 /*     .. Intrinsic Functions .. */
00254 /*     .. */
00255 /*     .. Scalars in Common .. */
00256 /*     .. */
00257 /*     .. Common blocks .. */
00258 /*     .. */
00259 /*     .. Data statements .. */
00260     /* Parameter adjustments */
00261     --iwork;
00262     --rwork;
00263     --work;
00264     --copys;
00265     --s;
00266     --c__;
00267     --copyb;
00268     --b;
00269     --copya;
00270     --a;
00271     --nxval;
00272     --nbval;
00273     --nsval;
00274     --nval;
00275     --mval;
00276     --dotype;
00277 
00278     /* Function Body */
00279 /*     .. */
00280 /*     .. Executable Statements .. */
00281 
00282 /*     Initialize constants and the random number seed. */
00283 
00284     s_copy(path, "Complex precision", (ftnlen)1, (ftnlen)17);
00285     s_copy(path + 1, "LS", (ftnlen)2, (ftnlen)2);
00286     nrun = 0;
00287     nfail = 0;
00288     nerrs = 0;
00289     for (i__ = 1; i__ <= 4; ++i__) {
00290         iseed[i__ - 1] = iseedy[i__ - 1];
00291 /* L10: */
00292     }
00293     eps = slamch_("Epsilon");
00294 
00295 /*     Threshold for rank estimation */
00296 
00297     rcond = sqrt(eps) - (sqrt(eps) - eps) / 2;
00298 
00299 /*     Test the error exits */
00300 
00301     xlaenv_(&c__9, &c__25);
00302     if (*tsterr) {
00303         cerrls_(path, nout);
00304     }
00305 
00306 /*     Print the header if NM = 0 or NN = 0 and THRESH = 0. */
00307 
00308     if ((*nm == 0 || *nn == 0) && *thresh == 0.f) {
00309         alahd_(nout, path);
00310     }
00311     infoc_1.infot = 0;
00312 
00313     i__1 = *nm;
00314     for (im = 1; im <= i__1; ++im) {
00315         m = mval[im];
00316         lda = max(1,m);
00317 
00318         i__2 = *nn;
00319         for (in = 1; in <= i__2; ++in) {
00320             n = nval[in];
00321             mnmin = min(m,n);
00322 /* Computing MAX */
00323             i__3 = max(1,m);
00324             ldb = max(i__3,n);
00325 
00326             i__3 = *nns;
00327             for (ins = 1; ins <= i__3; ++ins) {
00328                 nrhs = nsval[ins];
00329 /* Computing MAX */
00330                 i__4 = 1, i__5 = (m + nrhs) * (n + 2), i__4 = max(i__4,i__5), 
00331                         i__5 = (n + nrhs) * (m + 2), i__4 = max(i__4,i__5), 
00332                         i__5 = m * n + (mnmin << 2) + max(m,n), i__4 = max(
00333                         i__4,i__5), i__5 = (n << 1) + m;
00334                 lwork = max(i__4,i__5);
00335 
00336                 for (irank = 1; irank <= 2; ++irank) {
00337                     for (iscale = 1; iscale <= 3; ++iscale) {
00338                         itype = (irank - 1) * 3 + iscale;
00339                         if (! dotype[itype]) {
00340                             goto L100;
00341                         }
00342 
00343                         if (irank == 1) {
00344 
00345 /*                       Test CGELS */
00346 
00347 /*                       Generate a matrix of scaling type ISCALE */
00348 
00349                             cqrt13_(&iscale, &m, &n, &copya[1], &lda, &norma, 
00350                                     iseed);
00351                             i__4 = *nnb;
00352                             for (inb = 1; inb <= i__4; ++inb) {
00353                                 nb = nbval[inb];
00354                                 xlaenv_(&c__1, &nb);
00355                                 xlaenv_(&c__3, &nxval[inb]);
00356 
00357                                 for (itran = 1; itran <= 2; ++itran) {
00358                                     if (itran == 1) {
00359                                         *(unsigned char *)trans = 'N';
00360                                         nrows = m;
00361                                         ncols = n;
00362                                     } else {
00363                                         *(unsigned char *)trans = 'C';
00364                                         nrows = n;
00365                                         ncols = m;
00366                                     }
00367                                     ldwork = max(1,ncols);
00368 
00369 /*                             Set up a consistent rhs */
00370 
00371                                     if (ncols > 0) {
00372                                         i__5 = ncols * nrhs;
00373                                         clarnv_(&c__2, iseed, &i__5, &work[1])
00374                                                 ;
00375                                         i__5 = ncols * nrhs;
00376                                         r__1 = 1.f / (real) ncols;
00377                                         csscal_(&i__5, &r__1, &work[1], &c__1)
00378                                                 ;
00379                                     }
00380                                     cgemm_(trans, "No transpose", &nrows, &
00381                                             nrhs, &ncols, &c_b1, &copya[1], &
00382                                             lda, &work[1], &ldwork, &c_b2, &b[
00383                                             1], &ldb);
00384                                     clacpy_("Full", &nrows, &nrhs, &b[1], &
00385                                             ldb, &copyb[1], &ldb);
00386 
00387 /*                             Solve LS or overdetermined system */
00388 
00389                                     if (m > 0 && n > 0) {
00390                                         clacpy_("Full", &m, &n, &copya[1], &
00391                                                 lda, &a[1], &lda);
00392                                         clacpy_("Full", &nrows, &nrhs, &copyb[
00393                                                 1], &ldb, &b[1], &ldb);
00394                                     }
00395                                     s_copy(srnamc_1.srnamt, "CGELS ", (ftnlen)
00396                                             32, (ftnlen)6);
00397                                     cgels_(trans, &m, &n, &nrhs, &a[1], &lda, 
00398                                             &b[1], &ldb, &work[1], &lwork, &
00399                                             info);
00400 
00401                                     if (info != 0) {
00402                                         alaerh_(path, "CGELS ", &info, &c__0, 
00403                                                 trans, &m, &n, &nrhs, &c_n1, &
00404                                                 nb, &itype, &nfail, &nerrs, 
00405                                                 nout);
00406                                     }
00407 
00408 /*                             Check correctness of results */
00409 
00410                                     ldwork = max(1,nrows);
00411                                     if (nrows > 0 && nrhs > 0) {
00412                                         clacpy_("Full", &nrows, &nrhs, &copyb[
00413                                                 1], &ldb, &c__[1], &ldb);
00414                                     }
00415                                     cqrt16_(trans, &m, &n, &nrhs, &copya[1], &
00416                                             lda, &b[1], &ldb, &c__[1], &ldb, &
00417                                             rwork[1], result);
00418 
00419                                     if (itran == 1 && m >= n || itran == 2 && 
00420                                             m < n) {
00421 
00422 /*                                Solving LS system */
00423 
00424                                         result[1] = cqrt17_(trans, &c__1, &m, 
00425                                                 &n, &nrhs, &copya[1], &lda, &
00426                                                 b[1], &ldb, &copyb[1], &ldb, &
00427                                                 c__[1], &work[1], &lwork);
00428                                     } else {
00429 
00430 /*                                Solving overdetermined system */
00431 
00432                                         result[1] = cqrt14_(trans, &m, &n, &
00433                                                 nrhs, &copya[1], &lda, &b[1], 
00434                                                 &ldb, &work[1], &lwork);
00435                                     }
00436 
00437 /*                             Print information about the tests that */
00438 /*                             did not pass the threshold. */
00439 
00440                                     for (k = 1; k <= 2; ++k) {
00441                                         if (result[k - 1] >= *thresh) {
00442                                             if (nfail == 0 && nerrs == 0) {
00443                           alahd_(nout, path);
00444                                             }
00445                                             io___34.ciunit = *nout;
00446                                             s_wsfe(&io___34);
00447                                             do_fio(&c__1, trans, (ftnlen)1);
00448                                             do_fio(&c__1, (char *)&m, (ftnlen)
00449                                                     sizeof(integer));
00450                                             do_fio(&c__1, (char *)&n, (ftnlen)
00451                                                     sizeof(integer));
00452                                             do_fio(&c__1, (char *)&nrhs, (
00453                                                     ftnlen)sizeof(integer));
00454                                             do_fio(&c__1, (char *)&nb, (
00455                                                     ftnlen)sizeof(integer));
00456                                             do_fio(&c__1, (char *)&itype, (
00457                                                     ftnlen)sizeof(integer));
00458                                             do_fio(&c__1, (char *)&k, (ftnlen)
00459                                                     sizeof(integer));
00460                                             do_fio(&c__1, (char *)&result[k - 
00461                                                     1], (ftnlen)sizeof(real));
00462                                             e_wsfe();
00463                                             ++nfail;
00464                                         }
00465 /* L20: */
00466                                     }
00467                                     nrun += 2;
00468 /* L30: */
00469                                 }
00470 /* L40: */
00471                             }
00472                         }
00473 
00474 /*                    Generate a matrix of scaling type ISCALE and rank */
00475 /*                    type IRANK. */
00476 
00477                         cqrt15_(&iscale, &irank, &m, &n, &nrhs, &copya[1], &
00478                                 lda, &copyb[1], &ldb, &copys[1], &rank, &
00479                                 norma, &normb, iseed, &work[1], &lwork);
00480 
00481 /*                    workspace used: MAX(M+MIN(M,N),NRHS*MIN(M,N),2*N+M) */
00482 
00483                         i__4 = n;
00484                         for (j = 1; j <= i__4; ++j) {
00485                             iwork[j] = 0;
00486 /* L50: */
00487                         }
00488                         ldwork = max(1,m);
00489 
00490 /*                    Test CGELSX */
00491 
00492 /*                    CGELSX:  Compute the minimum-norm solution X */
00493 /*                    to min( norm( A * X - B ) ) */
00494 /*                    using a complete orthogonal factorization. */
00495 
00496                         clacpy_("Full", &m, &n, &copya[1], &lda, &a[1], &lda);
00497                         clacpy_("Full", &m, &nrhs, &copyb[1], &ldb, &b[1], &
00498                                 ldb);
00499 
00500                         s_copy(srnamc_1.srnamt, "CGELSX", (ftnlen)32, (ftnlen)
00501                                 6);
00502                         cgelsx_(&m, &n, &nrhs, &a[1], &lda, &b[1], &ldb, &
00503                                 iwork[1], &rcond, &crank, &work[1], &rwork[1], 
00504                                  &info);
00505 
00506                         if (info != 0) {
00507                             alaerh_(path, "CGELSX", &info, &c__0, " ", &m, &n, 
00508                                      &nrhs, &c_n1, &nb, &itype, &nfail, &
00509                                     nerrs, nout);
00510                         }
00511 
00512 /*                    workspace used: MAX( MNMIN+3*N, 2*MNMIN+NRHS ) */
00513 
00514 /*                    Test 3:  Compute relative error in svd */
00515 /*                             workspace: M*N + 4*MIN(M,N) + MAX(M,N) */
00516 
00517                         result[2] = cqrt12_(&crank, &crank, &a[1], &lda, &
00518                                 copys[1], &work[1], &lwork, &rwork[1]);
00519 
00520 /*                    Test 4:  Compute error in solution */
00521 /*                             workspace:  M*NRHS + M */
00522 
00523                         clacpy_("Full", &m, &nrhs, &copyb[1], &ldb, &work[1], 
00524                                 &ldwork);
00525                         cqrt16_("No transpose", &m, &n, &nrhs, &copya[1], &
00526                                 lda, &b[1], &ldb, &work[1], &ldwork, &rwork[1]
00527 , &result[3]);
00528 
00529 /*                    Test 5:  Check norm of r'*A */
00530 /*                             workspace: NRHS*(M+N) */
00531 
00532                         result[4] = 0.f;
00533                         if (m > crank) {
00534                             result[4] = cqrt17_("No transpose", &c__1, &m, &n, 
00535                                      &nrhs, &copya[1], &lda, &b[1], &ldb, &
00536                                     copyb[1], &ldb, &c__[1], &work[1], &lwork);
00537                         }
00538 
00539 /*                    Test 6:  Check if x is in the rowspace of A */
00540 /*                             workspace: (M+NRHS)*(N+2) */
00541 
00542                         result[5] = 0.f;
00543 
00544                         if (n > crank) {
00545                             result[5] = cqrt14_("No transpose", &m, &n, &nrhs, 
00546                                      &copya[1], &lda, &b[1], &ldb, &work[1], &
00547                                     lwork);
00548                         }
00549 
00550 /*                    Print information about the tests that did not */
00551 /*                    pass the threshold. */
00552 
00553                         for (k = 3; k <= 6; ++k) {
00554                             if (result[k - 1] >= *thresh) {
00555                                 if (nfail == 0 && nerrs == 0) {
00556                                     alahd_(nout, path);
00557                                 }
00558                                 io___39.ciunit = *nout;
00559                                 s_wsfe(&io___39);
00560                                 do_fio(&c__1, (char *)&m, (ftnlen)sizeof(
00561                                         integer));
00562                                 do_fio(&c__1, (char *)&n, (ftnlen)sizeof(
00563                                         integer));
00564                                 do_fio(&c__1, (char *)&nrhs, (ftnlen)sizeof(
00565                                         integer));
00566                                 do_fio(&c__1, (char *)&c__0, (ftnlen)sizeof(
00567                                         integer));
00568                                 do_fio(&c__1, (char *)&itype, (ftnlen)sizeof(
00569                                         integer));
00570                                 do_fio(&c__1, (char *)&k, (ftnlen)sizeof(
00571                                         integer));
00572                                 do_fio(&c__1, (char *)&result[k - 1], (ftnlen)
00573                                         sizeof(real));
00574                                 e_wsfe();
00575                                 ++nfail;
00576                             }
00577 /* L60: */
00578                         }
00579                         nrun += 4;
00580 
00581 /*                    Loop for testing different block sizes. */
00582 
00583                         i__4 = *nnb;
00584                         for (inb = 1; inb <= i__4; ++inb) {
00585                             nb = nbval[inb];
00586                             xlaenv_(&c__1, &nb);
00587                             xlaenv_(&c__3, &nxval[inb]);
00588 
00589 /*                       Test CGELSY */
00590 
00591 /*                       CGELSY:  Compute the minimum-norm solution */
00592 /*                       X to min( norm( A * X - B ) ) */
00593 /*                       using the rank-revealing orthogonal */
00594 /*                       factorization. */
00595 
00596                             clacpy_("Full", &m, &n, &copya[1], &lda, &a[1], &
00597                                     lda);
00598                             clacpy_("Full", &m, &nrhs, &copyb[1], &ldb, &b[1], 
00599                                      &ldb);
00600 
00601 /*                       Initialize vector IWORK. */
00602 
00603                             i__5 = n;
00604                             for (j = 1; j <= i__5; ++j) {
00605                                 iwork[j] = 0;
00606 /* L70: */
00607                             }
00608 
00609 /*                       Set LWLSY to the adequate value. */
00610 
00611 /* Computing MAX */
00612                             i__5 = mnmin << 1, i__6 = nb * (n + 1), i__5 = 
00613                                     max(i__5,i__6), i__6 = mnmin + nb * nrhs;
00614                             lwlsy = mnmin + max(i__5,i__6);
00615                             lwlsy = max(1,lwlsy);
00616 
00617                             s_copy(srnamc_1.srnamt, "CGELSY", (ftnlen)32, (
00618                                     ftnlen)6);
00619                             cgelsy_(&m, &n, &nrhs, &a[1], &lda, &b[1], &ldb, &
00620                                     iwork[1], &rcond, &crank, &work[1], &
00621                                     lwlsy, &rwork[1], &info);
00622                             if (info != 0) {
00623                                 alaerh_(path, "CGELSY", &info, &c__0, " ", &m, 
00624                                          &n, &nrhs, &c_n1, &nb, &itype, &
00625                                         nfail, &nerrs, nout);
00626                             }
00627 
00628 /*                       workspace used: 2*MNMIN+NB*NB+NB*MAX(N,NRHS) */
00629 
00630 /*                       Test 7:  Compute relative error in svd */
00631 /*                                workspace: M*N + 4*MIN(M,N) + MAX(M,N) */
00632 
00633                             result[6] = cqrt12_(&crank, &crank, &a[1], &lda, &
00634                                     copys[1], &work[1], &lwork, &rwork[1]);
00635 
00636 /*                       Test 8:  Compute error in solution */
00637 /*                                workspace:  M*NRHS + M */
00638 
00639                             clacpy_("Full", &m, &nrhs, &copyb[1], &ldb, &work[
00640                                     1], &ldwork);
00641                             cqrt16_("No transpose", &m, &n, &nrhs, &copya[1], 
00642                                     &lda, &b[1], &ldb, &work[1], &ldwork, &
00643                                     rwork[1], &result[7]);
00644 
00645 /*                       Test 9:  Check norm of r'*A */
00646 /*                                workspace: NRHS*(M+N) */
00647 
00648                             result[8] = 0.f;
00649                             if (m > crank) {
00650                                 result[8] = cqrt17_("No transpose", &c__1, &m, 
00651                                          &n, &nrhs, &copya[1], &lda, &b[1], &
00652                                         ldb, &copyb[1], &ldb, &c__[1], &work[
00653                                         1], &lwork);
00654                             }
00655 
00656 /*                       Test 10:  Check if x is in the rowspace of A */
00657 /*                                workspace: (M+NRHS)*(N+2) */
00658 
00659                             result[9] = 0.f;
00660 
00661                             if (n > crank) {
00662                                 result[9] = cqrt14_("No transpose", &m, &n, &
00663                                         nrhs, &copya[1], &lda, &b[1], &ldb, &
00664                                         work[1], &lwork);
00665                             }
00666 
00667 /*                       Test CGELSS */
00668 
00669 /*                       CGELSS:  Compute the minimum-norm solution */
00670 /*                       X to min( norm( A * X - B ) ) */
00671 /*                       using the SVD. */
00672 
00673                             clacpy_("Full", &m, &n, &copya[1], &lda, &a[1], &
00674                                     lda);
00675                             clacpy_("Full", &m, &nrhs, &copyb[1], &ldb, &b[1], 
00676                                      &ldb);
00677                             s_copy(srnamc_1.srnamt, "CGELSS", (ftnlen)32, (
00678                                     ftnlen)6);
00679                             cgelss_(&m, &n, &nrhs, &a[1], &lda, &b[1], &ldb, &
00680                                     s[1], &rcond, &crank, &work[1], &lwork, &
00681                                     rwork[1], &info);
00682 
00683                             if (info != 0) {
00684                                 alaerh_(path, "CGELSS", &info, &c__0, " ", &m, 
00685                                          &n, &nrhs, &c_n1, &nb, &itype, &
00686                                         nfail, &nerrs, nout);
00687                             }
00688 
00689 /*                       workspace used: 3*min(m,n) + */
00690 /*                                       max(2*min(m,n),nrhs,max(m,n)) */
00691 
00692 /*                       Test 11:  Compute relative error in svd */
00693 
00694                             if (rank > 0) {
00695                                 saxpy_(&mnmin, &c_b91, &copys[1], &c__1, &s[1]
00696 , &c__1);
00697                                 result[10] = sasum_(&mnmin, &s[1], &c__1) / 
00698                                         sasum_(&mnmin, &copys[1], &c__1) / (
00699                                         eps * (real) mnmin);
00700                             } else {
00701                                 result[10] = 0.f;
00702                             }
00703 
00704 /*                       Test 12:  Compute error in solution */
00705 
00706                             clacpy_("Full", &m, &nrhs, &copyb[1], &ldb, &work[
00707                                     1], &ldwork);
00708                             cqrt16_("No transpose", &m, &n, &nrhs, &copya[1], 
00709                                     &lda, &b[1], &ldb, &work[1], &ldwork, &
00710                                     rwork[1], &result[11]);
00711 
00712 /*                       Test 13:  Check norm of r'*A */
00713 
00714                             result[12] = 0.f;
00715                             if (m > crank) {
00716                                 result[12] = cqrt17_("No transpose", &c__1, &
00717                                         m, &n, &nrhs, &copya[1], &lda, &b[1], 
00718                                         &ldb, &copyb[1], &ldb, &c__[1], &work[
00719                                         1], &lwork);
00720                             }
00721 
00722 /*                       Test 14:  Check if x is in the rowspace of A */
00723 
00724                             result[13] = 0.f;
00725                             if (n > crank) {
00726                                 result[13] = cqrt14_("No transpose", &m, &n, &
00727                                         nrhs, &copya[1], &lda, &b[1], &ldb, &
00728                                         work[1], &lwork);
00729                             }
00730 
00731 /*                       Test CGELSD */
00732 
00733 /*                       CGELSD:  Compute the minimum-norm solution X */
00734 /*                       to min( norm( A * X - B ) ) using a */
00735 /*                       divide and conquer SVD. */
00736 
00737                             xlaenv_(&c__9, &c__25);
00738 
00739                             clacpy_("Full", &m, &n, &copya[1], &lda, &a[1], &
00740                                     lda);
00741                             clacpy_("Full", &m, &nrhs, &copyb[1], &ldb, &b[1], 
00742                                      &ldb);
00743 
00744                             s_copy(srnamc_1.srnamt, "CGELSD", (ftnlen)32, (
00745                                     ftnlen)6);
00746                             cgelsd_(&m, &n, &nrhs, &a[1], &lda, &b[1], &ldb, &
00747                                     s[1], &rcond, &crank, &work[1], &lwork, &
00748                                     rwork[1], &iwork[1], &info);
00749                             if (info != 0) {
00750                                 alaerh_(path, "CGELSD", &info, &c__0, " ", &m, 
00751                                          &n, &nrhs, &c_n1, &nb, &itype, &
00752                                         nfail, &nerrs, nout);
00753                             }
00754 
00755 /*                       Test 15:  Compute relative error in svd */
00756 
00757                             if (rank > 0) {
00758                                 saxpy_(&mnmin, &c_b91, &copys[1], &c__1, &s[1]
00759 , &c__1);
00760                                 result[14] = sasum_(&mnmin, &s[1], &c__1) / 
00761                                         sasum_(&mnmin, &copys[1], &c__1) / (
00762                                         eps * (real) mnmin);
00763                             } else {
00764                                 result[14] = 0.f;
00765                             }
00766 
00767 /*                       Test 16:  Compute error in solution */
00768 
00769                             clacpy_("Full", &m, &nrhs, &copyb[1], &ldb, &work[
00770                                     1], &ldwork);
00771                             cqrt16_("No transpose", &m, &n, &nrhs, &copya[1], 
00772                                     &lda, &b[1], &ldb, &work[1], &ldwork, &
00773                                     rwork[1], &result[15]);
00774 
00775 /*                       Test 17:  Check norm of r'*A */
00776 
00777                             result[16] = 0.f;
00778                             if (m > crank) {
00779                                 result[16] = cqrt17_("No transpose", &c__1, &
00780                                         m, &n, &nrhs, &copya[1], &lda, &b[1], 
00781                                         &ldb, &copyb[1], &ldb, &c__[1], &work[
00782                                         1], &lwork);
00783                             }
00784 
00785 /*                       Test 18:  Check if x is in the rowspace of A */
00786 
00787                             result[17] = 0.f;
00788                             if (n > crank) {
00789                                 result[17] = cqrt14_("No transpose", &m, &n, &
00790                                         nrhs, &copya[1], &lda, &b[1], &ldb, &
00791                                         work[1], &lwork);
00792                             }
00793 
00794 /*                       Print information about the tests that did not */
00795 /*                       pass the threshold. */
00796 
00797                             for (k = 7; k <= 18; ++k) {
00798                                 if (result[k - 1] >= *thresh) {
00799                                     if (nfail == 0 && nerrs == 0) {
00800                                         alahd_(nout, path);
00801                                     }
00802                                     io___41.ciunit = *nout;
00803                                     s_wsfe(&io___41);
00804                                     do_fio(&c__1, (char *)&m, (ftnlen)sizeof(
00805                                             integer));
00806                                     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(
00807                                             integer));
00808                                     do_fio(&c__1, (char *)&nrhs, (ftnlen)
00809                                             sizeof(integer));
00810                                     do_fio(&c__1, (char *)&nb, (ftnlen)sizeof(
00811                                             integer));
00812                                     do_fio(&c__1, (char *)&itype, (ftnlen)
00813                                             sizeof(integer));
00814                                     do_fio(&c__1, (char *)&k, (ftnlen)sizeof(
00815                                             integer));
00816                                     do_fio(&c__1, (char *)&result[k - 1], (
00817                                             ftnlen)sizeof(real));
00818                                     e_wsfe();
00819                                     ++nfail;
00820                                 }
00821 /* L80: */
00822                             }
00823                             nrun += 12;
00824 
00825 /* L90: */
00826                         }
00827 L100:
00828                         ;
00829                     }
00830 /* L110: */
00831                 }
00832 /* L120: */
00833             }
00834 /* L130: */
00835         }
00836 /* L140: */
00837     }
00838 
00839 /*     Print a summary of the results. */
00840 
00841     alasvm_(path, nout, &nfail, &nrun, &nerrs);
00842 
00843     return 0;
00844 
00845 /*     End of CDRVLS */
00846 
00847 } /* cdrvls_ */

swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:21