sdrvls.c
Go to the documentation of this file.
00001 /* sdrvls.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 integer c__2 = 2;
00034 static integer c__9 = 9;
00035 static integer c__25 = 25;
00036 static integer c__1 = 1;
00037 static integer c__3 = 3;
00038 static real c_b24 = 1.f;
00039 static real c_b25 = 0.f;
00040 static integer c__0 = 0;
00041 static integer c_n1 = -1;
00042 static real c_b92 = -1.f;
00043 
00044 /* Subroutine */ int sdrvls_(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         real *a, real *copya, real *b, real *copyb, real *c__, real *s, real *
00048         copys, real *work, integer *iwork, integer *nout)
00049 {
00050     /* Initialized data */
00051 
00052     static integer iseedy[4] = { 1988,1989,1990,1991 };
00053 
00054     /* Format strings */
00055     static char fmt_9999[] = "(\002 TRANS='\002,a1,\002', M=\002,i5,\002, N"
00056             "=\002,i5,\002, NRHS=\002,i4,\002, NB=\002,i4,\002, type\002,i2"
00057             ",\002, test(\002,i2,\002)=\002,g12.5)";
00058     static char fmt_9998[] = "(\002 M=\002,i5,\002, N=\002,i5,\002, NRHS="
00059             "\002,i4,\002, NB=\002,i4,\002, type\002,i2,\002, test(\002,i2"
00060             ",\002)=\002,g12.5)";
00061 
00062     /* System generated locals */
00063     integer i__1, i__2, i__3, i__4, i__5, i__6;
00064     real r__1, r__2;
00065 
00066     /* Builtin functions */
00067     /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
00068     double sqrt(doublereal), log(doublereal);
00069     integer s_wsfe(cilist *), do_fio(integer *, char *, ftnlen), e_wsfe(void);
00070 
00071     /* Local variables */
00072     integer i__, j, k, m, n, nb, im, in, lda, ldb, inb;
00073     real eps;
00074     integer ins, info;
00075     char path[3];
00076     integer rank, nrhs, nlvl, nrun;
00077     extern /* Subroutine */ int alahd_(integer *, char *);
00078     integer nfail, iseed[4], crank, irank;
00079     real rcond;
00080     extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *), 
00081             sgemm_(char *, char *, integer *, integer *, integer *, real *, 
00082             real *, integer *, real *, integer *, real *, real *, integer *);
00083     integer itran, mnmin, ncols;
00084     real norma, normb;
00085     extern /* Subroutine */ int sgels_(char *, integer *, integer *, integer *
00086 , real *, integer *, real *, integer *, real *, integer *, 
00087             integer *);
00088     char trans[1];
00089     integer nerrs, itype;
00090     extern doublereal sasum_(integer *, real *, integer *);
00091     integer lwork;
00092     extern doublereal sqrt12_(integer *, integer *, real *, integer *, real *, 
00093              real *, integer *), sqrt14_(char *, integer *, integer *, 
00094             integer *, real *, integer *, real *, integer *, real *, integer *
00095 ), sqrt17_(char *, integer *, integer *, integer *, 
00096             integer *, real *, integer *, real *, integer *, real *, integer *
00097 , real *, real *, integer *);
00098     extern /* Subroutine */ int sqrt13_(integer *, integer *, integer *, real 
00099             *, integer *, real *, integer *), sqrt15_(integer *, integer *, 
00100             integer *, integer *, integer *, real *, integer *, real *, 
00101             integer *, real *, integer *, real *, real *, integer *, real *, 
00102             integer *), saxpy_(integer *, real *, real *, integer *, real *, 
00103             integer *), sqrt16_(char *, integer *, integer *, integer *, real 
00104             *, integer *, real *, integer *, real *, integer *, real *, real *
00105 );
00106     integer nrows, lwlsy;
00107     extern /* Subroutine */ int alaerh_(char *, char *, integer *, integer *, 
00108             char *, integer *, integer *, integer *, integer *, integer *, 
00109             integer *, integer *, integer *, integer *);
00110     integer iscale;
00111     extern doublereal slamch_(char *);
00112     extern /* Subroutine */ int sgelsd_(integer *, integer *, integer *, real 
00113             *, integer *, real *, integer *, real *, real *, integer *, real *
00114 , integer *, integer *, integer *), alasvm_(char *, integer *, 
00115             integer *, integer *, integer *), slacpy_(char *, integer 
00116             *, integer *, real *, integer *, real *, integer *), 
00117             xlaenv_(integer *, integer *), sgelss_(integer *, integer *, 
00118             integer *, real *, integer *, real *, integer *, real *, real *, 
00119             integer *, real *, integer *, integer *);
00120     integer ldwork;
00121     extern /* Subroutine */ int sgelsx_(integer *, integer *, integer *, real 
00122             *, integer *, real *, integer *, integer *, real *, integer *, 
00123             real *, integer *), sgelsy_(integer *, integer *, integer *, real 
00124             *, integer *, real *, integer *, integer *, real *, integer *, 
00125             real *, integer *, integer *), slarnv_(integer *, integer *, 
00126             integer *, real *), serrls_(char *, integer *);
00127     real result[18];
00128 
00129     /* Fortran I/O blocks */
00130     static cilist io___35 = { 0, 0, 0, fmt_9999, 0 };
00131     static cilist io___40 = { 0, 0, 0, fmt_9998, 0 };
00132     static cilist io___42 = { 0, 0, 0, fmt_9998, 0 };
00133 
00134 
00135 
00136 /*  -- LAPACK test routine (version 3.1.1) -- */
00137 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00138 /*     January 2007 */
00139 
00140 /*     .. Scalar Arguments .. */
00141 /*     .. */
00142 /*     .. Array Arguments .. */
00143 /*     .. */
00144 
00145 /*  Purpose */
00146 /*  ======= */
00147 
00148 /*  SDRVLS tests the least squares driver routines SGELS, SGELSS, SGELSX, */
00149 /*  SGELSY and SGELSD. */
00150 
00151 /*  Arguments */
00152 /*  ========= */
00153 
00154 /*  DOTYPE  (input) LOGICAL array, dimension (NTYPES) */
00155 /*          The matrix types to be used for testing.  Matrices of type j */
00156 /*          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) = */
00157 /*          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used. */
00158 /*          The matrix of type j is generated as follows: */
00159 /*          j=1: A = U*D*V where U and V are random orthogonal matrices */
00160 /*               and D has random entries (> 0.1) taken from a uniform */
00161 /*               distribution (0,1). A is full rank. */
00162 /*          j=2: The same of 1, but A is scaled up. */
00163 /*          j=3: The same of 1, but A is scaled down. */
00164 /*          j=4: A = U*D*V where U and V are random orthogonal matrices */
00165 /*               and D has 3*min(M,N)/4 random entries (> 0.1) taken */
00166 /*               from a uniform distribution (0,1) and the remaining */
00167 /*               entries set to 0. A is rank-deficient. */
00168 /*          j=5: The same of 4, but A is scaled up. */
00169 /*          j=6: The same of 5, but A is scaled down. */
00170 
00171 /*  NM      (input) INTEGER */
00172 /*          The number of values of M contained in the vector MVAL. */
00173 
00174 /*  MVAL    (input) INTEGER array, dimension (NM) */
00175 /*          The values of the matrix row dimension M. */
00176 
00177 /*  NN      (input) INTEGER */
00178 /*          The number of values of N contained in the vector NVAL. */
00179 
00180 /*  NVAL    (input) INTEGER array, dimension (NN) */
00181 /*          The values of the matrix column dimension N. */
00182 
00183 /*  NNS     (input) INTEGER */
00184 /*          The number of values of NRHS contained in the vector NSVAL. */
00185 
00186 /*  NSVAL   (input) INTEGER array, dimension (NNS) */
00187 /*          The values of the number of right hand sides NRHS. */
00188 
00189 /*  NNB     (input) INTEGER */
00190 /*          The number of values of NB and NX contained in the */
00191 /*          vectors NBVAL and NXVAL.  The blocking parameters are used */
00192 /*          in pairs (NB,NX). */
00193 
00194 /*  NBVAL   (input) INTEGER array, dimension (NNB) */
00195 /*          The values of the blocksize NB. */
00196 
00197 /*  NXVAL   (input) INTEGER array, dimension (NNB) */
00198 /*          The values of the crossover point NX. */
00199 
00200 /*  THRESH  (input) REAL */
00201 /*          The threshold value for the test ratios.  A result is */
00202 /*          included in the output file if RESULT >= THRESH.  To have */
00203 /*          every test ratio printed, use THRESH = 0. */
00204 
00205 /*  TSTERR  (input) LOGICAL */
00206 /*          Flag that indicates whether error exits are to be tested. */
00207 
00208 /*  A       (workspace) REAL array, dimension (MMAX*NMAX) */
00209 /*          where MMAX is the maximum value of M in MVAL and NMAX is the */
00210 /*          maximum value of N in NVAL. */
00211 
00212 /*  COPYA   (workspace) REAL array, dimension (MMAX*NMAX) */
00213 
00214 /*  B       (workspace) REAL array, dimension (MMAX*NSMAX) */
00215 /*          where MMAX is the maximum value of M in MVAL and NSMAX is the */
00216 /*          maximum value of NRHS in NSVAL. */
00217 
00218 /*  COPYB   (workspace) REAL array, dimension (MMAX*NSMAX) */
00219 
00220 /*  C       (workspace) REAL array, dimension (MMAX*NSMAX) */
00221 
00222 /*  S       (workspace) REAL array, dimension */
00223 /*                      (min(MMAX,NMAX)) */
00224 
00225 /*  COPYS   (workspace) REAL array, dimension */
00226 /*                      (min(MMAX,NMAX)) */
00227 
00228 /*  WORK    (workspace) REAL array, */
00229 /*                      dimension (MMAX*NMAX + 4*NMAX + MMAX). */
00230 
00231 /*  IWORK   (workspace) INTEGER array, dimension (15*NMAX) */
00232 
00233 /*  NOUT    (input) INTEGER */
00234 /*          The unit number for output. */
00235 
00236 /*  ===================================================================== */
00237 
00238 /*     .. Parameters .. */
00239 /*     .. */
00240 /*     .. Local Scalars .. */
00241 /*     .. */
00242 /*     .. Local Arrays .. */
00243 /*     .. */
00244 /*     .. External Functions .. */
00245 /*     .. */
00246 /*     .. External Subroutines .. */
00247 /*     .. */
00248 /*     .. Intrinsic Functions .. */
00249 /*     .. */
00250 /*     .. Scalars in Common .. */
00251 /*     .. */
00252 /*     .. Common blocks .. */
00253 /*     .. */
00254 /*     .. Data statements .. */
00255     /* Parameter adjustments */
00256     --iwork;
00257     --work;
00258     --copys;
00259     --s;
00260     --c__;
00261     --copyb;
00262     --b;
00263     --copya;
00264     --a;
00265     --nxval;
00266     --nbval;
00267     --nsval;
00268     --nval;
00269     --mval;
00270     --dotype;
00271 
00272     /* Function Body */
00273 /*     .. */
00274 /*     .. Executable Statements .. */
00275 
00276 /*     Initialize constants and the random number seed. */
00277 
00278     s_copy(path, "Single precision", (ftnlen)1, (ftnlen)16);
00279     s_copy(path + 1, "LS", (ftnlen)2, (ftnlen)2);
00280     nrun = 0;
00281     nfail = 0;
00282     nerrs = 0;
00283     for (i__ = 1; i__ <= 4; ++i__) {
00284         iseed[i__ - 1] = iseedy[i__ - 1];
00285 /* L10: */
00286     }
00287     eps = slamch_("Epsilon");
00288 
00289 /*     Threshold for rank estimation */
00290 
00291     rcond = sqrt(eps) - (sqrt(eps) - eps) / 2;
00292 
00293 /*     Test the error exits */
00294 
00295     xlaenv_(&c__2, &c__2);
00296     xlaenv_(&c__9, &c__25);
00297     if (*tsterr) {
00298         serrls_(path, nout);
00299     }
00300 
00301 /*     Print the header if NM = 0 or NN = 0 and THRESH = 0. */
00302 
00303     if ((*nm == 0 || *nn == 0) && *thresh == 0.f) {
00304         alahd_(nout, path);
00305     }
00306     infoc_1.infot = 0;
00307 
00308     i__1 = *nm;
00309     for (im = 1; im <= i__1; ++im) {
00310         m = mval[im];
00311         lda = max(1,m);
00312 
00313         i__2 = *nn;
00314         for (in = 1; in <= i__2; ++in) {
00315             n = nval[in];
00316             mnmin = min(m,n);
00317 /* Computing MAX */
00318             i__3 = max(1,m);
00319             ldb = max(i__3,n);
00320 
00321             i__3 = *nns;
00322             for (ins = 1; ins <= i__3; ++ins) {
00323                 nrhs = nsval[ins];
00324 /* Computing MAX */
00325 /* Computing MAX */
00326                 r__1 = 1.f, r__2 = (real) mnmin;
00327                 i__4 = (integer) (log(dmax(r__1,r__2) / 26.f) / log(2.f)) + 1;
00328                 nlvl = max(i__4,0);
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 = mnmin * 12 + mnmin * 50 + (mnmin <<
00334                          3) * nlvl + mnmin * nrhs + 676;
00335                 lwork = max(i__4,i__5);
00336 
00337                 for (irank = 1; irank <= 2; ++irank) {
00338                     for (iscale = 1; iscale <= 3; ++iscale) {
00339                         itype = (irank - 1) * 3 + iscale;
00340                         if (! dotype[itype]) {
00341                             goto L110;
00342                         }
00343 
00344                         if (irank == 1) {
00345 
00346 /*                       Test SGELS */
00347 
00348 /*                       Generate a matrix of scaling type ISCALE */
00349 
00350                             sqrt13_(&iscale, &m, &n, &copya[1], &lda, &norma, 
00351                                     iseed);
00352                             i__4 = *nnb;
00353                             for (inb = 1; inb <= i__4; ++inb) {
00354                                 nb = nbval[inb];
00355                                 xlaenv_(&c__1, &nb);
00356                                 xlaenv_(&c__3, &nxval[inb]);
00357 
00358                                 for (itran = 1; itran <= 2; ++itran) {
00359                                     if (itran == 1) {
00360                                         *(unsigned char *)trans = 'N';
00361                                         nrows = m;
00362                                         ncols = n;
00363                                     } else {
00364                                         *(unsigned char *)trans = 'T';
00365                                         nrows = n;
00366                                         ncols = m;
00367                                     }
00368                                     ldwork = max(1,ncols);
00369 
00370 /*                             Set up a consistent rhs */
00371 
00372                                     if (ncols > 0) {
00373                                         i__5 = ncols * nrhs;
00374                                         slarnv_(&c__2, iseed, &i__5, &work[1])
00375                                                 ;
00376                                         i__5 = ncols * nrhs;
00377                                         r__1 = 1.f / (real) ncols;
00378                                         sscal_(&i__5, &r__1, &work[1], &c__1);
00379                                     }
00380                                     sgemm_(trans, "No transpose", &nrows, &
00381                                             nrhs, &ncols, &c_b24, &copya[1], &
00382                                             lda, &work[1], &ldwork, &c_b25, &
00383                                             b[1], &ldb)
00384                                             ;
00385                                     slacpy_("Full", &nrows, &nrhs, &b[1], &
00386                                             ldb, &copyb[1], &ldb);
00387 
00388 /*                             Solve LS or overdetermined system */
00389 
00390                                     if (m > 0 && n > 0) {
00391                                         slacpy_("Full", &m, &n, &copya[1], &
00392                                                 lda, &a[1], &lda);
00393                                         slacpy_("Full", &nrows, &nrhs, &copyb[
00394                                                 1], &ldb, &b[1], &ldb);
00395                                     }
00396                                     s_copy(srnamc_1.srnamt, "SGELS ", (ftnlen)
00397                                             32, (ftnlen)6);
00398                                     sgels_(trans, &m, &n, &nrhs, &a[1], &lda, 
00399                                             &b[1], &ldb, &work[1], &lwork, &
00400                                             info);
00401                                     if (info != 0) {
00402                                         alaerh_(path, "SGELS ", &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                                         slacpy_("Full", &nrows, &nrhs, &copyb[
00413                                                 1], &ldb, &c__[1], &ldb);
00414                                     }
00415                                     sqrt16_(trans, &m, &n, &nrhs, &copya[1], &
00416                                             lda, &b[1], &ldb, &c__[1], &ldb, &
00417                                             work[1], result);
00418 
00419                                     if (itran == 1 && m >= n || itran == 2 && 
00420                                             m < n) {
00421 
00422 /*                                Solving LS system */
00423 
00424                                         result[1] = sqrt17_(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] = sqrt14_(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___35.ciunit = *nout;
00446                                             s_wsfe(&io___35);
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                         sqrt15_(&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 /*                    Initialize vector IWORK. */
00484 
00485                         i__4 = n;
00486                         for (j = 1; j <= i__4; ++j) {
00487                             iwork[j] = 0;
00488 /* L50: */
00489                         }
00490                         ldwork = max(1,m);
00491 
00492 /*                    Test SGELSX */
00493 
00494 /*                    SGELSX:  Compute the minimum-norm solution X */
00495 /*                    to min( norm( A * X - B ) ) using a complete */
00496 /*                    orthogonal factorization. */
00497 
00498                         slacpy_("Full", &m, &n, &copya[1], &lda, &a[1], &lda);
00499                         slacpy_("Full", &m, &nrhs, &copyb[1], &ldb, &b[1], &
00500                                 ldb);
00501 
00502                         s_copy(srnamc_1.srnamt, "SGELSX", (ftnlen)32, (ftnlen)
00503                                 6);
00504                         sgelsx_(&m, &n, &nrhs, &a[1], &lda, &b[1], &ldb, &
00505                                 iwork[1], &rcond, &crank, &work[1], &info);
00506                         if (info != 0) {
00507                             alaerh_(path, "SGELSX", &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] = sqrt12_(&crank, &crank, &a[1], &lda, &
00518                                 copys[1], &work[1], &lwork);
00519 
00520 /*                    Test 4:  Compute error in solution */
00521 /*                             workspace:  M*NRHS + M */
00522 
00523                         slacpy_("Full", &m, &nrhs, &copyb[1], &ldb, &work[1], 
00524                                 &ldwork);
00525                         sqrt16_("No transpose", &m, &n, &nrhs, &copya[1], &
00526                                 lda, &b[1], &ldb, &work[1], &ldwork, &work[m *
00527                                  nrhs + 1], &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] = sqrt17_("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] = sqrt14_("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___40.ciunit = *nout;
00559                                 s_wsfe(&io___40);
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 *)&nb, (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 SGELSY */
00590 
00591 /*                       SGELSY:  Compute the minimum-norm solution X */
00592 /*                       to min( norm( A * X - B ) ) */
00593 /*                       using the rank-revealing orthogonal */
00594 /*                       factorization. */
00595 
00596 /*                       Initialize vector IWORK. */
00597 
00598                             i__5 = n;
00599                             for (j = 1; j <= i__5; ++j) {
00600                                 iwork[j] = 0;
00601 /* L70: */
00602                             }
00603 
00604 /*                       Set LWLSY to the adequate value. */
00605 
00606 /* Computing MAX */
00607                             i__5 = 1, i__6 = mnmin + (n << 1) + nb * (n + 1), 
00608                                     i__5 = max(i__5,i__6), i__6 = (mnmin << 1)
00609                                      + nb * nrhs;
00610                             lwlsy = max(i__5,i__6);
00611 
00612                             slacpy_("Full", &m, &n, &copya[1], &lda, &a[1], &
00613                                     lda);
00614                             slacpy_("Full", &m, &nrhs, &copyb[1], &ldb, &b[1], 
00615                                      &ldb);
00616 
00617                             s_copy(srnamc_1.srnamt, "SGELSY", (ftnlen)32, (
00618                                     ftnlen)6);
00619                             sgelsy_(&m, &n, &nrhs, &a[1], &lda, &b[1], &ldb, &
00620                                     iwork[1], &rcond, &crank, &work[1], &
00621                                     lwlsy, &info);
00622                             if (info != 0) {
00623                                 alaerh_(path, "SGELSY", &info, &c__0, " ", &m, 
00624                                          &n, &nrhs, &c_n1, &nb, &itype, &
00625                                         nfail, &nerrs, nout);
00626                             }
00627 
00628 /*                       Test 7:  Compute relative error in svd */
00629 /*                                workspace: M*N + 4*MIN(M,N) + MAX(M,N) */
00630 
00631                             result[6] = sqrt12_(&crank, &crank, &a[1], &lda, &
00632                                     copys[1], &work[1], &lwork);
00633 
00634 /*                       Test 8:  Compute error in solution */
00635 /*                                workspace:  M*NRHS + M */
00636 
00637                             slacpy_("Full", &m, &nrhs, &copyb[1], &ldb, &work[
00638                                     1], &ldwork);
00639                             sqrt16_("No transpose", &m, &n, &nrhs, &copya[1], 
00640                                     &lda, &b[1], &ldb, &work[1], &ldwork, &
00641                                     work[m * nrhs + 1], &result[7]);
00642 
00643 /*                       Test 9:  Check norm of r'*A */
00644 /*                                workspace: NRHS*(M+N) */
00645 
00646                             result[8] = 0.f;
00647                             if (m > crank) {
00648                                 result[8] = sqrt17_("No transpose", &c__1, &m, 
00649                                          &n, &nrhs, &copya[1], &lda, &b[1], &
00650                                         ldb, &copyb[1], &ldb, &c__[1], &work[
00651                                         1], &lwork);
00652                             }
00653 
00654 /*                       Test 10:  Check if x is in the rowspace of A */
00655 /*                                workspace: (M+NRHS)*(N+2) */
00656 
00657                             result[9] = 0.f;
00658 
00659                             if (n > crank) {
00660                                 result[9] = sqrt14_("No transpose", &m, &n, &
00661                                         nrhs, &copya[1], &lda, &b[1], &ldb, &
00662                                         work[1], &lwork);
00663                             }
00664 
00665 /*                       Test SGELSS */
00666 
00667 /*                       SGELSS:  Compute the minimum-norm solution X */
00668 /*                       to min( norm( A * X - B ) ) */
00669 /*                       using the SVD. */
00670 
00671                             slacpy_("Full", &m, &n, &copya[1], &lda, &a[1], &
00672                                     lda);
00673                             slacpy_("Full", &m, &nrhs, &copyb[1], &ldb, &b[1], 
00674                                      &ldb);
00675                             s_copy(srnamc_1.srnamt, "SGELSS", (ftnlen)32, (
00676                                     ftnlen)6);
00677                             sgelss_(&m, &n, &nrhs, &a[1], &lda, &b[1], &ldb, &
00678                                     s[1], &rcond, &crank, &work[1], &lwork, &
00679                                     info);
00680                             if (info != 0) {
00681                                 alaerh_(path, "SGELSS", &info, &c__0, " ", &m, 
00682                                          &n, &nrhs, &c_n1, &nb, &itype, &
00683                                         nfail, &nerrs, nout);
00684                             }
00685 
00686 /*                       workspace used: 3*min(m,n) + */
00687 /*                                       max(2*min(m,n),nrhs,max(m,n)) */
00688 
00689 /*                       Test 11:  Compute relative error in svd */
00690 
00691                             if (rank > 0) {
00692                                 saxpy_(&mnmin, &c_b92, &copys[1], &c__1, &s[1]
00693 , &c__1);
00694                                 result[10] = sasum_(&mnmin, &s[1], &c__1) / 
00695                                         sasum_(&mnmin, &copys[1], &c__1) / (
00696                                         eps * (real) mnmin);
00697                             } else {
00698                                 result[10] = 0.f;
00699                             }
00700 
00701 /*                       Test 12:  Compute error in solution */
00702 
00703                             slacpy_("Full", &m, &nrhs, &copyb[1], &ldb, &work[
00704                                     1], &ldwork);
00705                             sqrt16_("No transpose", &m, &n, &nrhs, &copya[1], 
00706                                     &lda, &b[1], &ldb, &work[1], &ldwork, &
00707                                     work[m * nrhs + 1], &result[11]);
00708 
00709 /*                       Test 13:  Check norm of r'*A */
00710 
00711                             result[12] = 0.f;
00712                             if (m > crank) {
00713                                 result[12] = sqrt17_("No transpose", &c__1, &
00714                                         m, &n, &nrhs, &copya[1], &lda, &b[1], 
00715                                         &ldb, &copyb[1], &ldb, &c__[1], &work[
00716                                         1], &lwork);
00717                             }
00718 
00719 /*                       Test 14:  Check if x is in the rowspace of A */
00720 
00721                             result[13] = 0.f;
00722                             if (n > crank) {
00723                                 result[13] = sqrt14_("No transpose", &m, &n, &
00724                                         nrhs, &copya[1], &lda, &b[1], &ldb, &
00725                                         work[1], &lwork);
00726                             }
00727 
00728 /*                       Test SGELSD */
00729 
00730 /*                       SGELSD:  Compute the minimum-norm solution X */
00731 /*                       to min( norm( A * X - B ) ) using a */
00732 /*                       divide and conquer SVD. */
00733 
00734 /*                       Initialize vector IWORK. */
00735 
00736                             i__5 = n;
00737                             for (j = 1; j <= i__5; ++j) {
00738                                 iwork[j] = 0;
00739 /* L80: */
00740                             }
00741 
00742                             slacpy_("Full", &m, &n, &copya[1], &lda, &a[1], &
00743                                     lda);
00744                             slacpy_("Full", &m, &nrhs, &copyb[1], &ldb, &b[1], 
00745                                      &ldb);
00746 
00747                             s_copy(srnamc_1.srnamt, "SGELSD", (ftnlen)32, (
00748                                     ftnlen)6);
00749                             sgelsd_(&m, &n, &nrhs, &a[1], &lda, &b[1], &ldb, &
00750                                     s[1], &rcond, &crank, &work[1], &lwork, &
00751                                     iwork[1], &info);
00752                             if (info != 0) {
00753                                 alaerh_(path, "SGELSD", &info, &c__0, " ", &m, 
00754                                          &n, &nrhs, &c_n1, &nb, &itype, &
00755                                         nfail, &nerrs, nout);
00756                             }
00757 
00758 /*                       Test 15:  Compute relative error in svd */
00759 
00760                             if (rank > 0) {
00761                                 saxpy_(&mnmin, &c_b92, &copys[1], &c__1, &s[1]
00762 , &c__1);
00763                                 result[14] = sasum_(&mnmin, &s[1], &c__1) / 
00764                                         sasum_(&mnmin, &copys[1], &c__1) / (
00765                                         eps * (real) mnmin);
00766                             } else {
00767                                 result[14] = 0.f;
00768                             }
00769 
00770 /*                       Test 16:  Compute error in solution */
00771 
00772                             slacpy_("Full", &m, &nrhs, &copyb[1], &ldb, &work[
00773                                     1], &ldwork);
00774                             sqrt16_("No transpose", &m, &n, &nrhs, &copya[1], 
00775                                     &lda, &b[1], &ldb, &work[1], &ldwork, &
00776                                     work[m * nrhs + 1], &result[15]);
00777 
00778 /*                       Test 17:  Check norm of r'*A */
00779 
00780                             result[16] = 0.f;
00781                             if (m > crank) {
00782                                 result[16] = sqrt17_("No transpose", &c__1, &
00783                                         m, &n, &nrhs, &copya[1], &lda, &b[1], 
00784                                         &ldb, &copyb[1], &ldb, &c__[1], &work[
00785                                         1], &lwork);
00786                             }
00787 
00788 /*                       Test 18:  Check if x is in the rowspace of A */
00789 
00790                             result[17] = 0.f;
00791                             if (n > crank) {
00792                                 result[17] = sqrt14_("No transpose", &m, &n, &
00793                                         nrhs, &copya[1], &lda, &b[1], &ldb, &
00794                                         work[1], &lwork);
00795                             }
00796 
00797 /*                       Print information about the tests that did not */
00798 /*                       pass the threshold. */
00799 
00800                             for (k = 7; k <= 18; ++k) {
00801                                 if (result[k - 1] >= *thresh) {
00802                                     if (nfail == 0 && nerrs == 0) {
00803                                         alahd_(nout, path);
00804                                     }
00805                                     io___42.ciunit = *nout;
00806                                     s_wsfe(&io___42);
00807                                     do_fio(&c__1, (char *)&m, (ftnlen)sizeof(
00808                                             integer));
00809                                     do_fio(&c__1, (char *)&n, (ftnlen)sizeof(
00810                                             integer));
00811                                     do_fio(&c__1, (char *)&nrhs, (ftnlen)
00812                                             sizeof(integer));
00813                                     do_fio(&c__1, (char *)&nb, (ftnlen)sizeof(
00814                                             integer));
00815                                     do_fio(&c__1, (char *)&itype, (ftnlen)
00816                                             sizeof(integer));
00817                                     do_fio(&c__1, (char *)&k, (ftnlen)sizeof(
00818                                             integer));
00819                                     do_fio(&c__1, (char *)&result[k - 1], (
00820                                             ftnlen)sizeof(real));
00821                                     e_wsfe();
00822                                     ++nfail;
00823                                 }
00824 /* L90: */
00825                             }
00826                             nrun += 12;
00827 
00828 /* L100: */
00829                         }
00830 L110:
00831                         ;
00832                     }
00833 /* L120: */
00834                 }
00835 /* L130: */
00836             }
00837 /* L140: */
00838         }
00839 /* L150: */
00840     }
00841 
00842 /*     Print a summary of the results. */
00843 
00844     alasvm_(path, nout, &nfail, &nrun, &nerrs);
00845 
00846     return 0;
00847 
00848 /*     End of SDRVLS */
00849 
00850 } /* sdrvls_ */

swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:00