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

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