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


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