123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302 |
- ///////////////////////////////////////////////////////////////////////////////
- // Copyright 2011 John Maddock. Distributed under the Boost
- // Software License, Version 1.0. (See accompanying file
- // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
- /* 1000d.f -- translated by f2c (version 20050501).
- You must link the resulting object file with libf2c:
- on Microsoft Windows system, link with libf2c.lib;
- on Linux or Unix systems, link with .../path/to/libf2c.a -lm
- or, if you install libf2c.a in a standard place, with -lf2c -lm
- -- in that order, at the end of the command line, as in
- cc *.o -lf2c -lm
- Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
- http://www.netlib.org/f2c/libf2c.zip
- */
- #include <iostream>
- #include <iomanip>
- #include <cmath>
- #if defined(TEST_GMPXX)
- #include <gmpxx.h>
- typedef mpf_class real_type;
- #elif defined(TEST_MPFRXX)
- #include <gmpfrxx.h>
- typedef mpfr_class real_type;
- #elif defined(TEST_CPP_DEC_FLOAT)
- #include <boost/multiprecision/cpp_dec_float.hpp>
- typedef boost::multiprecision::cpp_dec_float_50 real_type;
- #elif defined(TEST_MPFR_50)
- #include <boost/multiprecision/mpfr.hpp>
- typedef boost::multiprecision::mpfr_float_50 real_type;
- #elif defined(TEST_MPF_50)
- #include <boost/multiprecision/gmp.hpp>
- typedef boost::multiprecision::mpf_float_50 real_type;
- #elif defined(NATIVE_FLOAT128)
- #include <boost/multiprecision/float128.hpp>
- typedef __float128 real_type;
- std::ostream& operator<<(std::ostream& os, const __float128& f)
- {
- return os << boost::multiprecision::float128(f);
- }
- #include <boost/type_traits/has_left_shift.hpp>
- namespace boost {
- template <>
- struct has_left_shift<std::basic_ostream<char>, __float128> : public mpl::true_
- {};
- template <>
- double lexical_cast<double, __float128>(const __float128& f)
- {
- return f;
- }
- } // namespace boost
- #elif defined(TEST_FLOAT128)
- #include <boost/multiprecision/float128.hpp>
- typedef boost::multiprecision::float128 real_type;
- #elif defined(TEST_CPP_BIN_FLOAT_QUAD)
- #include <boost/multiprecision/cpp_bin_float.hpp>
- typedef boost::multiprecision::cpp_bin_float_quad real_type;
- #elif defined(TEST_CPP_BIN_FLOAT_OCT)
- #include <boost/multiprecision/cpp_bin_float.hpp>
- typedef boost::multiprecision::cpp_bin_float_oct real_type;
- #else
- typedef double real_type;
- #endif
- #include <boost/lexical_cast.hpp>
- #ifndef CAST_TO_RT
- #define CAST_TO_RT(x) x
- #endif
- extern "C" {
- #include "f2c.h"
- integer s_wsfe(cilist*), e_wsfe(void), do_fio(integer*, char*, ftnlen),
- s_wsle(cilist*), do_lio(integer*, integer*, char*, ftnlen),
- e_wsle(void);
- /* Subroutine */ int s_stop(char*, ftnlen);
- #undef abs
- #undef dabs
- #define dabs abs
- #undef dmin
- #undef dmax
- #define dmin min
- #define dmax max
- }
- #include <time.h>
- using std::max;
- using std::min;
- /* Table of constant values */
- static integer c__0 = 0;
- static real_type c_b7 = CAST_TO_RT(1);
- static integer c__1 = 1;
- static integer c__9 = 9;
- inline double second_(void)
- {
- return ((double)(clock())) / CLOCKS_PER_SEC;
- }
- int dgefa_(real_type*, integer*, integer*, integer*, integer*), dgesl_(real_type*, integer*, integer*, integer*, real_type*, integer*);
- int dmxpy_(integer*, real_type*, integer*, integer*, real_type*, real_type*);
- int matgen_(real_type*, integer*, integer*, real_type*, real_type*);
- real_type epslon_(real_type*);
- real_type ran_(integer*);
- int dscal_(integer*, real_type*, real_type*, integer*);
- int daxpy_(integer*, real_type*, real_type*, integer*, real_type*, integer*);
- integer idamax_(integer*, real_type*, integer*);
- real_type ddot_(integer*, real_type*, integer*, real_type*, integer*);
- int daxpy_(integer*, real_type*, real_type*, integer*, real_type*, integer*);
- int dmxpy_(integer*, real_type*, integer*, integer*, real_type*, real_type*);
- extern "C" int MAIN__()
- {
- #ifdef TEST_MPF_50
- std::cout << "Testing number<mpf_float<50> >" << std::endl;
- #elif defined(TEST_MPFR_50)
- std::cout << "Testing number<mpf_float<50> >" << std::endl;
- #elif defined(TEST_GMPXX)
- std::cout << "Testing mpf_class at 50 decimal degits" << std::endl;
- mpf_set_default_prec(((50 + 1) * 1000L) / 301L);
- #elif defined(TEST_MPFRXX)
- std::cout << "Testing mpfr_class at 50 decimal degits" << std::endl;
- mpfr_set_default_prec(((50 + 1) * 1000L) / 301L);
- #elif defined(TEST_CPP_DEC_FLOAT)
- std::cout << "Testing number<cpp_dec_float<50> >" << std::endl;
- #elif defined(NATIVE_FLOAT128)
- std::cout << "Testing __float128" << std::endl;
- #elif defined(TEST_FLOAT128)
- std::cout << "Testing number<float128_backend, et_off>" << std::endl;
- #else
- std::cout << "Testing double" << std::endl;
- #endif
- /* Format strings */
- static char fmt_1[] = "(\002 Please send the results of this run to:\002"
- "//\002 Jack J. Dongarra\002/\002 Computer Science Department\002/"
- "\002 University of Tennessee\002/\002 Knoxville, Tennessee 37996"
- "-1300\002//\002 Fax: 615-974-8296\002//\002 Internet: dongarra@c"
- "s.utk.edu\002/)";
- static char fmt_40[] = "(\002 norm. resid resid mac"
- "hep\002,\002 x(1) x(n)\002)";
- static char fmt_50[] = "(1p5e16.8)";
- static char fmt_60[] = "(//\002 times are reported for matrices of or"
- "der \002,i5)";
- static char fmt_70[] = "(6x,\002factor\002,5x,\002solve\002,6x,\002tota"
- "l\002,5x,\002mflops\002,7x,\002unit\002,6x,\002ratio\002)";
- static char fmt_80[] = "(\002 times for array with leading dimension o"
- "f\002,i4)";
- static char fmt_110[] = "(6(1pe11.3))";
- /* System generated locals */
- integer i__1;
- real_type d__1, d__2, d__3;
- /* Builtin functions */
- /* Local variables */
- static real_type a[1001000] /* was [1001][1000] */, b[1000];
- static integer i__, n;
- static real_type x[1000];
- static double t1;
- static integer lda;
- static double ops;
- static real_type eps;
- static integer info;
- static double time[6], cray, total;
- static integer ipvt[1000];
- static real_type resid, norma;
- static real_type normx;
- static real_type residn;
- /* Fortran I/O blocks */
- static cilist io___4 = {0, 6, 0, fmt_1, 0};
- static cilist io___20 = {0, 6, 0, fmt_40, 0};
- static cilist io___21 = {0, 6, 0, fmt_50, 0};
- static cilist io___22 = {0, 6, 0, fmt_60, 0};
- static cilist io___23 = {0, 6, 0, fmt_70, 0};
- static cilist io___24 = {0, 6, 0, fmt_80, 0};
- static cilist io___25 = {0, 6, 0, fmt_110, 0};
- static cilist io___26 = {0, 6, 0, 0, 0};
- lda = 1001;
- /* this program was updated on 10/12/92 to correct a */
- /* problem with the random number generator. The previous */
- /* random number generator had a short period and produced */
- /* singular matrices occasionally. */
- n = 1000;
- cray = .056f;
- s_wsfe(&io___4);
- e_wsfe();
- /* Computing 3rd power */
- d__1 = (real_type)n;
- /* Computing 2nd power */
- d__2 = (real_type)n;
- ops = boost::lexical_cast<double>(real_type(d__1 * (d__1 * d__1) * 2. / 3. + d__2 * d__2 * 2.));
- matgen_(a, &lda, &n, b, &norma);
- /* ****************************************************************** */
- /* ****************************************************************** */
- /* you should replace the call to dgefa and dgesl */
- /* by calls to your linear equation solver. */
- /* ****************************************************************** */
- /* ****************************************************************** */
- t1 = second_();
- dgefa_(a, &lda, &n, ipvt, &info);
- time[0] = second_() - t1;
- t1 = second_();
- dgesl_(a, &lda, &n, ipvt, b, &c__0);
- time[1] = second_() - t1;
- total = time[0] + time[1];
- /* ****************************************************************** */
- /* ****************************************************************** */
- /* compute a residual to verify results. */
- i__1 = n;
- for (i__ = 1; i__ <= i__1; ++i__)
- {
- x[i__ - 1] = b[i__ - 1];
- /* L10: */
- }
- matgen_(a, &lda, &n, b, &norma);
- i__1 = n;
- for (i__ = 1; i__ <= i__1; ++i__)
- {
- b[i__ - 1] = -b[i__ - 1];
- /* L20: */
- }
- dmxpy_(&n, b, &n, &lda, x, a);
- resid = CAST_TO_RT(0);
- normx = CAST_TO_RT(0);
- i__1 = n;
- for (i__ = 1; i__ <= i__1; ++i__)
- {
- /* Computing MAX */
- d__2 = resid, d__3 = (d__1 = b[i__ - 1], abs(d__1));
- resid = (max)(d__2, d__3);
- /* Computing MAX */
- d__2 = normx, d__3 = (d__1 = x[i__ - 1], abs(d__1));
- normx = (max)(d__2, d__3);
- /* L30: */
- }
- eps = epslon_(&c_b7);
- residn = resid / (n * norma * normx * eps);
- s_wsfe(&io___20);
- e_wsfe();
- s_wsfe(&io___21);
- /*
- do_fio(&c__1, (char *)&residn, (ftnlen)sizeof(real_type));
- do_fio(&c__1, (char *)&resid, (ftnlen)sizeof(real_type));
- do_fio(&c__1, (char *)&eps, (ftnlen)sizeof(real_type));
- do_fio(&c__1, (char *)&x[0], (ftnlen)sizeof(real_type));
- do_fio(&c__1, (char *)&x[n - 1], (ftnlen)sizeof(real_type));
- */
- std::cout << std::setw(12) << std::setprecision(5) << residn << " " << resid << " " << eps << " " << x[0] << " " << x[n - 1] << std::endl;
- e_wsfe();
- s_wsfe(&io___22);
- do_fio(&c__1, (char*)&n, (ftnlen)sizeof(integer));
- e_wsfe();
- s_wsfe(&io___23);
- e_wsfe();
- time[2] = total;
- time[3] = ops / (total * 1e6);
- time[4] = 2. / time[3];
- time[5] = total / cray;
- s_wsfe(&io___24);
- do_fio(&c__1, (char*)&lda, (ftnlen)sizeof(integer));
- e_wsfe();
- s_wsfe(&io___25);
- for (i__ = 1; i__ <= 6; ++i__)
- {
- // do_fio(&c__1, (char *)&time[i__ - 1], (ftnlen)sizeof(real_type));
- std::cout << std::setw(12) << std::setprecision(5) << time[i__ - 1];
- }
- e_wsfe();
- s_wsle(&io___26);
- do_lio(&c__9, &c__1, " end of tests -- this version dated 10/12/92", (ftnlen)44);
- e_wsle();
- s_stop("", (ftnlen)0);
- return 0;
- } /* MAIN__ */
- /* Subroutine */ int matgen_(real_type* a, integer* lda, integer* n,
- real_type* b, real_type* norma)
- {
- /* System generated locals */
- integer a_dim1, a_offset, i__1, i__2;
- real_type d__1, d__2;
- /* Local variables */
- static integer i__, j;
- static integer init[4];
- /* Parameter adjustments */
- a_dim1 = *lda;
- a_offset = 1 + a_dim1;
- a -= a_offset;
- --b;
- /* Function Body */
- init[0] = 1;
- init[1] = 2;
- init[2] = 3;
- init[3] = 1325;
- *norma = CAST_TO_RT(0);
- i__1 = *n;
- for (j = 1; j <= i__1; ++j)
- {
- i__2 = *n;
- for (i__ = 1; i__ <= i__2; ++i__)
- {
- a[i__ + j * a_dim1] = ran_(init) - .5f;
- /* Computing MAX */
- d__2 = (d__1 = a[i__ + j * a_dim1], abs(d__1));
- *norma = (max)(d__2, *norma);
- /* L20: */
- }
- /* L30: */
- }
- i__1 = *n;
- for (i__ = 1; i__ <= i__1; ++i__)
- {
- b[i__] = CAST_TO_RT(0);
- /* L35: */
- }
- i__1 = *n;
- for (j = 1; j <= i__1; ++j)
- {
- i__2 = *n;
- for (i__ = 1; i__ <= i__2; ++i__)
- {
- b[i__] += a[i__ + j * a_dim1];
- /* L40: */
- }
- /* L50: */
- }
- return 0;
- } /* matgen_ */
- /* Subroutine */ int dgefa_(real_type* a, integer* lda, integer* n, integer* ipvt, integer* info)
- {
- /* System generated locals */
- integer a_dim1, a_offset, i__1, i__2, i__3;
- /* Local variables */
- static integer j, k, l;
- static real_type t;
- static integer kp1, nm1;
- /* dgefa factors a double precision matrix by gaussian elimination. */
- /* dgefa is usually called by dgeco, but it can be called */
- /* directly with a saving in time if rcond is not needed. */
- /* (time for dgeco) = (1 + 9/n)*(time for dgefa) . */
- /* on entry */
- /* a double precision(lda, n) */
- /* the matrix to be factored. */
- /* lda integer */
- /* the leading dimension of the array a . */
- /* n integer */
- /* the order of the matrix a . */
- /* on return */
- /* a an upper triangular matrix and the multipliers */
- /* which were used to obtain it. */
- /* the factorization can be written a = l*u where */
- /* l is a product of permutation and unit lower */
- /* triangular matrices and u is upper triangular. */
- /* ipvt integer(n) */
- /* an integer vector of pivot indices. */
- /* info integer */
- /* = 0 normal value. */
- /* = k if u(k,k) .eq. 0.0 . this is not an error */
- /* condition for this subroutine, but it does */
- /* indicate that dgesl or dgedi will divide by zero */
- /* if called. use rcond in dgeco for a reliable */
- /* indication of singularity. */
- /* linpack. this version dated 08/14/78 . */
- /* cleve moler, university of new mexico, argonne national lab. */
- /* subroutines and functions */
- /* blas daxpy,dscal,idamax */
- /* internal variables */
- /* gaussian elimination with partial pivoting */
- /* Parameter adjustments */
- a_dim1 = *lda;
- a_offset = 1 + a_dim1;
- a -= a_offset;
- --ipvt;
- /* Function Body */
- *info = 0;
- nm1 = *n - 1;
- if (nm1 < 1)
- {
- goto L70;
- }
- i__1 = nm1;
- for (k = 1; k <= i__1; ++k)
- {
- kp1 = k + 1;
- /* find l = pivot index */
- i__2 = *n - k + 1;
- l = idamax_(&i__2, &a[k + k * a_dim1], &c__1) + k - 1;
- ipvt[k] = l;
- /* zero pivot implies this column already triangularized */
- if (a[l + k * a_dim1] == 0.)
- {
- goto L40;
- }
- /* interchange if necessary */
- if (l == k)
- {
- goto L10;
- }
- t = a[l + k * a_dim1];
- a[l + k * a_dim1] = a[k + k * a_dim1];
- a[k + k * a_dim1] = t;
- L10:
- /* compute multipliers */
- t = -1. / a[k + k * a_dim1];
- i__2 = *n - k;
- dscal_(&i__2, &t, &a[k + 1 + k * a_dim1], &c__1);
- /* row elimination with column indexing */
- i__2 = *n;
- for (j = kp1; j <= i__2; ++j)
- {
- t = a[l + j * a_dim1];
- if (l == k)
- {
- goto L20;
- }
- a[l + j * a_dim1] = a[k + j * a_dim1];
- a[k + j * a_dim1] = t;
- L20:
- i__3 = *n - k;
- daxpy_(&i__3, &t, &a[k + 1 + k * a_dim1], &c__1, &a[k + 1 + j * a_dim1], &c__1);
- /* L30: */
- }
- goto L50;
- L40:
- *info = k;
- L50:
- /* L60: */
- ;
- }
- L70:
- ipvt[*n] = *n;
- if (a[*n + *n * a_dim1] == 0.)
- {
- *info = *n;
- }
- return 0;
- } /* dgefa_ */
- /* Subroutine */ int dgesl_(real_type* a, integer* lda, integer* n, integer* ipvt, real_type* b, integer* job)
- {
- /* System generated locals */
- integer a_dim1, a_offset, i__1, i__2;
- /* Local variables */
- static integer k, l;
- static real_type t;
- static integer kb, nm1;
- /* dgesl solves the double precision system */
- /* a * x = b or trans(a) * x = b */
- /* using the factors computed by dgeco or dgefa. */
- /* on entry */
- /* a double precision(lda, n) */
- /* the output from dgeco or dgefa. */
- /* lda integer */
- /* the leading dimension of the array a . */
- /* n integer */
- /* the order of the matrix a . */
- /* ipvt integer(n) */
- /* the pivot vector from dgeco or dgefa. */
- /* b double precision(n) */
- /* the right hand side vector. */
- /* job integer */
- /* = 0 to solve a*x = b , */
- /* = nonzero to solve trans(a)*x = b where */
- /* trans(a) is the transpose. */
- /* on return */
- /* b the solution vector x . */
- /* error condition */
- /* a division by zero will occur if the input factor contains a */
- /* zero on the diagonal. technically this indicates singularity */
- /* but it is often caused by improper arguments or improper */
- /* setting of lda . it will not occur if the subroutines are */
- /* called correctly and if dgeco has set rcond .gt. 0.0 */
- /* or dgefa has set info .eq. 0 . */
- /* to compute inverse(a) * c where c is a matrix */
- /* with p columns */
- /* call dgeco(a,lda,n,ipvt,rcond,z) */
- /* if (rcond is too small) go to ... */
- /* do 10 j = 1, p */
- /* call dgesl(a,lda,n,ipvt,c(1,j),0) */
- /* 10 continue */
- /* linpack. this version dated 08/14/78 . */
- /* cleve moler, university of new mexico, argonne national lab. */
- /* subroutines and functions */
- /* blas daxpy,ddot */
- /* internal variables */
- /* Parameter adjustments */
- a_dim1 = *lda;
- a_offset = 1 + a_dim1;
- a -= a_offset;
- --ipvt;
- --b;
- /* Function Body */
- nm1 = *n - 1;
- if (*job != 0)
- {
- goto L50;
- }
- /* job = 0 , solve a * x = b */
- /* first solve l*y = b */
- if (nm1 < 1)
- {
- goto L30;
- }
- i__1 = nm1;
- for (k = 1; k <= i__1; ++k)
- {
- l = ipvt[k];
- t = b[l];
- if (l == k)
- {
- goto L10;
- }
- b[l] = b[k];
- b[k] = t;
- L10:
- i__2 = *n - k;
- daxpy_(&i__2, &t, &a[k + 1 + k * a_dim1], &c__1, &b[k + 1], &c__1);
- /* L20: */
- }
- L30:
- /* now solve u*x = y */
- i__1 = *n;
- for (kb = 1; kb <= i__1; ++kb)
- {
- k = *n + 1 - kb;
- b[k] /= a[k + k * a_dim1];
- t = -b[k];
- i__2 = k - 1;
- daxpy_(&i__2, &t, &a[k * a_dim1 + 1], &c__1, &b[1], &c__1);
- /* L40: */
- }
- goto L100;
- L50:
- /* job = nonzero, solve trans(a) * x = b */
- /* first solve trans(u)*y = b */
- i__1 = *n;
- for (k = 1; k <= i__1; ++k)
- {
- i__2 = k - 1;
- t = ddot_(&i__2, &a[k * a_dim1 + 1], &c__1, &b[1], &c__1);
- b[k] = (b[k] - t) / a[k + k * a_dim1];
- /* L60: */
- }
- /* now solve trans(l)*x = y */
- if (nm1 < 1)
- {
- goto L90;
- }
- i__1 = nm1;
- for (kb = 1; kb <= i__1; ++kb)
- {
- k = *n - kb;
- i__2 = *n - k;
- b[k] += ddot_(&i__2, &a[k + 1 + k * a_dim1], &c__1, &b[k + 1], &c__1);
- l = ipvt[k];
- if (l == k)
- {
- goto L70;
- }
- t = b[l];
- b[l] = b[k];
- b[k] = t;
- L70:
- /* L80: */
- ;
- }
- L90:
- L100:
- return 0;
- } /* dgesl_ */
- /* Subroutine */ int daxpy_(integer* n, real_type* da, real_type* dx,
- integer* incx, real_type* dy, integer* incy)
- {
- /* System generated locals */
- integer i__1;
- /* Local variables */
- static integer i__, m, ix, iy, mp1;
- /* constant times a vector plus a vector. */
- /* uses unrolled loops for increments equal to one. */
- /* jack dongarra, linpack, 3/11/78. */
- /* Parameter adjustments */
- --dy;
- --dx;
- /* Function Body */
- if (*n <= 0)
- {
- return 0;
- }
- if (*da == 0.)
- {
- return 0;
- }
- if (*incx == 1 && *incy == 1)
- {
- goto L20;
- }
- /* code for unequal increments or equal increments */
- /* not equal to 1 */
- ix = 1;
- iy = 1;
- if (*incx < 0)
- {
- ix = (-(*n) + 1) * *incx + 1;
- }
- if (*incy < 0)
- {
- iy = (-(*n) + 1) * *incy + 1;
- }
- i__1 = *n;
- for (i__ = 1; i__ <= i__1; ++i__)
- {
- dy[iy] += *da * dx[ix];
- ix += *incx;
- iy += *incy;
- /* L10: */
- }
- return 0;
- /* code for both increments equal to 1 */
- /* clean-up loop */
- L20:
- m = *n % 4;
- if (m == 0)
- {
- goto L40;
- }
- i__1 = m;
- for (i__ = 1; i__ <= i__1; ++i__)
- {
- dy[i__] += *da * dx[i__];
- /* L30: */
- }
- if (*n < 4)
- {
- return 0;
- }
- L40:
- mp1 = m + 1;
- i__1 = *n;
- for (i__ = mp1; i__ <= i__1; i__ += 4)
- {
- dy[i__] += *da * dx[i__];
- dy[i__ + 1] += *da * dx[i__ + 1];
- dy[i__ + 2] += *da * dx[i__ + 2];
- dy[i__ + 3] += *da * dx[i__ + 3];
- /* L50: */
- }
- return 0;
- } /* daxpy_ */
- real_type ddot_(integer* n, real_type* dx, integer* incx, real_type* dy,
- integer* incy)
- {
- /* System generated locals */
- integer i__1;
- real_type ret_val;
- /* Local variables */
- static integer i__, m, ix, iy, mp1;
- static real_type dtemp;
- /* forms the dot product of two vectors. */
- /* uses unrolled loops for increments equal to one. */
- /* jack dongarra, linpack, 3/11/78. */
- /* Parameter adjustments */
- --dy;
- --dx;
- /* Function Body */
- ret_val = CAST_TO_RT(0);
- dtemp = CAST_TO_RT(0);
- if (*n <= 0)
- {
- return ret_val;
- }
- if (*incx == 1 && *incy == 1)
- {
- goto L20;
- }
- /* code for unequal increments or equal increments */
- /* not equal to 1 */
- ix = 1;
- iy = 1;
- if (*incx < 0)
- {
- ix = (-(*n) + 1) * *incx + 1;
- }
- if (*incy < 0)
- {
- iy = (-(*n) + 1) * *incy + 1;
- }
- i__1 = *n;
- for (i__ = 1; i__ <= i__1; ++i__)
- {
- dtemp += dx[ix] * dy[iy];
- ix += *incx;
- iy += *incy;
- /* L10: */
- }
- ret_val = dtemp;
- return ret_val;
- /* code for both increments equal to 1 */
- /* clean-up loop */
- L20:
- m = *n % 5;
- if (m == 0)
- {
- goto L40;
- }
- i__1 = m;
- for (i__ = 1; i__ <= i__1; ++i__)
- {
- dtemp += dx[i__] * dy[i__];
- /* L30: */
- }
- if (*n < 5)
- {
- goto L60;
- }
- L40:
- mp1 = m + 1;
- i__1 = *n;
- for (i__ = mp1; i__ <= i__1; i__ += 5)
- {
- dtemp = dtemp + dx[i__] * dy[i__] + dx[i__ + 1] * dy[i__ + 1] + dx[i__ + 2] * dy[i__ + 2] + dx[i__ + 3] * dy[i__ + 3] + dx[i__ + 4] * dy[i__ + 4];
- /* L50: */
- }
- L60:
- ret_val = dtemp;
- return ret_val;
- } /* ddot_ */
- /* Subroutine */ int dscal_(integer* n, real_type* da, real_type* dx,
- integer* incx)
- {
- /* System generated locals */
- integer i__1, i__2;
- /* Local variables */
- static integer i__, m, mp1, nincx;
- /* scales a vector by a constant. */
- /* uses unrolled loops for increment equal to one. */
- /* jack dongarra, linpack, 3/11/78. */
- /* Parameter adjustments */
- --dx;
- /* Function Body */
- if (*n <= 0)
- {
- return 0;
- }
- if (*incx == 1)
- {
- goto L20;
- }
- /* code for increment not equal to 1 */
- nincx = *n * *incx;
- i__1 = nincx;
- i__2 = *incx;
- for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2)
- {
- dx[i__] = *da * dx[i__];
- /* L10: */
- }
- return 0;
- /* code for increment equal to 1 */
- /* clean-up loop */
- L20:
- m = *n % 5;
- if (m == 0)
- {
- goto L40;
- }
- i__2 = m;
- for (i__ = 1; i__ <= i__2; ++i__)
- {
- dx[i__] = *da * dx[i__];
- /* L30: */
- }
- if (*n < 5)
- {
- return 0;
- }
- L40:
- mp1 = m + 1;
- i__2 = *n;
- for (i__ = mp1; i__ <= i__2; i__ += 5)
- {
- dx[i__] = *da * dx[i__];
- dx[i__ + 1] = *da * dx[i__ + 1];
- dx[i__ + 2] = *da * dx[i__ + 2];
- dx[i__ + 3] = *da * dx[i__ + 3];
- dx[i__ + 4] = *da * dx[i__ + 4];
- /* L50: */
- }
- return 0;
- } /* dscal_ */
- integer idamax_(integer* n, real_type* dx, integer* incx)
- {
- /* System generated locals */
- integer ret_val, i__1;
- real_type d__1;
- /* Local variables */
- static integer i__, ix;
- static real_type dmax__;
- /* finds the index of element having max. dabsolute value. */
- /* jack dongarra, linpack, 3/11/78. */
- /* Parameter adjustments */
- --dx;
- /* Function Body */
- ret_val = 0;
- if (*n < 1)
- {
- return ret_val;
- }
- ret_val = 1;
- if (*n == 1)
- {
- return ret_val;
- }
- if (*incx == 1)
- {
- goto L20;
- }
- /* code for increment not equal to 1 */
- ix = 1;
- dmax__ = abs(dx[1]);
- ix += *incx;
- i__1 = *n;
- for (i__ = 2; i__ <= i__1; ++i__)
- {
- if ((d__1 = dx[ix], abs(d__1)) <= dmax__)
- {
- goto L5;
- }
- ret_val = i__;
- dmax__ = (d__1 = dx[ix], abs(d__1));
- L5:
- ix += *incx;
- /* L10: */
- }
- return ret_val;
- /* code for increment equal to 1 */
- L20:
- dmax__ = abs(dx[1]);
- i__1 = *n;
- for (i__ = 2; i__ <= i__1; ++i__)
- {
- if ((d__1 = dx[i__], abs(d__1)) <= dmax__)
- {
- goto L30;
- }
- ret_val = i__;
- dmax__ = (d__1 = dx[i__], abs(d__1));
- L30:;
- }
- return ret_val;
- } /* idamax_ */
- real_type epslon_(real_type* x)
- {
- #if defined(TEST_MPF_100) || defined(TEST_MPFR_100) || defined(TEST_GMPXX) || defined(TEST_MPFRXX)
- return std::ldexp(1.0, 1 - ((100 + 1) * 1000L) / 301L);
- #elif defined(TEST_CPP_DEC_FLOAT_BN)
- return std::pow(10.0, 1 - std::numeric_limits<efx::cpp_dec_float_50>::digits10);
- #elif defined(NATIVE_FLOAT128)
- return FLT128_EPSILON;
- #else
- return CAST_TO_RT(std::numeric_limits<real_type>::epsilon());
- #endif
- } /* epslon_ */
- /* Subroutine */ int mm_(real_type* a, integer* lda, integer* n1, integer* n3, real_type* b, integer* ldb, integer* n2, real_type* c__,
- integer* ldc)
- {
- /* System generated locals */
- integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2;
- /* Local variables */
- static integer i__, j;
- /* purpose: */
- /* multiply matrix b times matrix c and store the result in matrix a. */
- /* parameters: */
- /* a double precision(lda,n3), matrix of n1 rows and n3 columns */
- /* lda integer, leading dimension of array a */
- /* n1 integer, number of rows in matrices a and b */
- /* n3 integer, number of columns in matrices a and c */
- /* b double precision(ldb,n2), matrix of n1 rows and n2 columns */
- /* ldb integer, leading dimension of array b */
- /* n2 integer, number of columns in matrix b, and number of rows in */
- /* matrix c */
- /* c double precision(ldc,n3), matrix of n2 rows and n3 columns */
- /* ldc integer, leading dimension of array c */
- /* ---------------------------------------------------------------------- */
- /* Parameter adjustments */
- a_dim1 = *lda;
- a_offset = 1 + a_dim1;
- a -= a_offset;
- b_dim1 = *ldb;
- b_offset = 1 + b_dim1;
- b -= b_offset;
- c_dim1 = *ldc;
- c_offset = 1 + c_dim1;
- c__ -= c_offset;
- /* Function Body */
- i__1 = *n3;
- for (j = 1; j <= i__1; ++j)
- {
- i__2 = *n1;
- for (i__ = 1; i__ <= i__2; ++i__)
- {
- a[i__ + j * a_dim1] = CAST_TO_RT(0);
- /* L10: */
- }
- dmxpy_(n2, &a[j * a_dim1 + 1], n1, ldb, &c__[j * c_dim1 + 1], &b[b_offset]);
- /* L20: */
- }
- return 0;
- } /* mm_ */
- /* Subroutine */ int dmxpy_(integer* n1, real_type* y, integer* n2, integer* ldm, real_type* x, real_type* m)
- {
- /* System generated locals */
- integer m_dim1, m_offset, i__1, i__2;
- /* Local variables */
- static integer i__, j, jmin;
- /* purpose: */
- /* multiply matrix m times vector x and add the result to vector y. */
- /* parameters: */
- /* n1 integer, number of elements in vector y, and number of rows in */
- /* matrix m */
- /* y double precision(n1), vector of length n1 to which is added */
- /* the product m*x */
- /* n2 integer, number of elements in vector x, and number of columns */
- /* in matrix m */
- /* ldm integer, leading dimension of array m */
- /* x double precision(n2), vector of length n2 */
- /* m double precision(ldm,n2), matrix of n1 rows and n2 columns */
- /* ---------------------------------------------------------------------- */
- /* cleanup odd vector */
- /* Parameter adjustments */
- --y;
- m_dim1 = *ldm;
- m_offset = 1 + m_dim1;
- m -= m_offset;
- --x;
- /* Function Body */
- j = *n2 % 2;
- if (j >= 1)
- {
- i__1 = *n1;
- for (i__ = 1; i__ <= i__1; ++i__)
- {
- y[i__] += x[j] * m[i__ + j * m_dim1];
- /* L10: */
- }
- }
- /* cleanup odd group of two vectors */
- j = *n2 % 4;
- if (j >= 2)
- {
- i__1 = *n1;
- for (i__ = 1; i__ <= i__1; ++i__)
- {
- y[i__] = y[i__] + x[j - 1] * m[i__ + (j - 1) * m_dim1] + x[j] * m[i__ + j * m_dim1];
- /* L20: */
- }
- }
- /* cleanup odd group of four vectors */
- j = *n2 % 8;
- if (j >= 4)
- {
- i__1 = *n1;
- for (i__ = 1; i__ <= i__1; ++i__)
- {
- y[i__] = y[i__] + x[j - 3] * m[i__ + (j - 3) * m_dim1] + x[j - 2] * m[i__ + (j - 2) * m_dim1] + x[j - 1] * m[i__ + (j - 1) * m_dim1] + x[j] * m[i__ + j * m_dim1];
- /* L30: */
- }
- }
- /* cleanup odd group of eight vectors */
- j = *n2 % 16;
- if (j >= 8)
- {
- i__1 = *n1;
- for (i__ = 1; i__ <= i__1; ++i__)
- {
- y[i__] = y[i__] + x[j - 7] * m[i__ + (j - 7) * m_dim1] + x[j - 6] * m[i__ + (j - 6) * m_dim1] + x[j - 5] * m[i__ + (j - 5) * m_dim1] + x[j - 4] * m[i__ + (j - 4) * m_dim1] + x[j - 3] * m[i__ + (j - 3) * m_dim1] + x[j - 2] * m[i__ + (j - 2) * m_dim1] + x[j - 1] * m[i__ + (j - 1) * m_dim1] + x[j] * m[i__ + j * m_dim1];
- /* L40: */
- }
- }
- /* main loop - groups of sixteen vectors */
- jmin = j + 16;
- i__1 = *n2;
- for (j = jmin; j <= i__1; j += 16)
- {
- i__2 = *n1;
- for (i__ = 1; i__ <= i__2; ++i__)
- {
- y[i__] = y[i__] + x[j - 15] * m[i__ + (j - 15) * m_dim1] + x[j - 14] * m[i__ + (j - 14) * m_dim1] + x[j - 13] * m[i__ + (j - 13) * m_dim1] + x[j - 12] * m[i__ + (j - 12) * m_dim1] + x[j - 11] * m[i__ + (j - 11) * m_dim1] + x[j - 10] * m[i__ + (j - 10) * m_dim1] + x[j - 9] * m[i__ + (j - 9) * m_dim1] + x[j - 8] * m[i__ + (j - 8) * m_dim1] + x[j - 7] * m[i__ + (j - 7) * m_dim1] + x[j - 6] * m[i__ + (j - 6) * m_dim1] + x[j - 5] * m[i__ + (j - 5) * m_dim1] + x[j - 4] * m[i__ + (j - 4) * m_dim1] + x[j - 3] * m[i__ + (j - 3) * m_dim1] + x[j - 2] * m[i__ + (j - 2) * m_dim1] + x[j - 1] * m[i__ + (j - 1) * m_dim1] + x[j] * m[i__ + j * m_dim1];
- /* L50: */
- }
- /* L60: */
- }
- return 0;
- } /* dmxpy_ */
- real_type ran_(integer* iseed)
- {
- /* System generated locals */
- real_type ret_val;
- /* Local variables */
- static integer it1, it2, it3, it4;
- /* modified from the LAPACK auxiliary routine 10/12/92 JD */
- /* -- LAPACK auxiliary routine (version 1.0) -- */
- /* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
- /* Courant Institute, Argonne National Lab, and Rice University */
- /* February 29, 1992 */
- /* .. Array Arguments .. */
- /* .. */
- /* Purpose */
- /* ======= */
- /* DLARAN returns a random double number from a uniform (0,1) */
- /* distribution. */
- /* Arguments */
- /* ========= */
- /* ISEED (input/output) INTEGER array, dimension (4) */
- /* On entry, the seed of the random number generator; the array */
- /* elements must be between 0 and 4095, and ISEED(4) must be */
- /* odd. */
- /* On exit, the seed is updated. */
- /* Further Details */
- /* =============== */
- /* This routine uses a multiplicative congruential method with modulus */
- /* 2**48 and multiplier 33952834046453 (see G.S.Fishman, */
- /* 'Multiplicative congruential random number generators with modulus */
- /* 2**b: an exhaustive analysis for b = 32 and a partial analysis for */
- /* b = 48', Math. Comp. 189, pp 331-344, 1990). */
- /* 48-bit integers are stored in 4 integer array elements with 12 bits */
- /* per element. Hence the routine is portable across machines with */
- /* integers of 32 bits or more. */
- /* .. Parameters .. */
- /* .. */
- /* .. Local Scalars .. */
- /* .. */
- /* .. Intrinsic Functions .. */
- /* .. */
- /* .. Executable Statements .. */
- /* multiply the seed by the multiplier modulo 2**48 */
- /* Parameter adjustments */
- --iseed;
- /* Function Body */
- it4 = iseed[4] * 2549;
- it3 = it4 / 4096;
- it4 -= it3 << 12;
- it3 = it3 + iseed[3] * 2549 + iseed[4] * 2508;
- it2 = it3 / 4096;
- it3 -= it2 << 12;
- it2 = it2 + iseed[2] * 2549 + iseed[3] * 2508 + iseed[4] * 322;
- it1 = it2 / 4096;
- it2 -= it1 << 12;
- it1 = it1 + iseed[1] * 2549 + iseed[2] * 2508 + iseed[3] * 322 + iseed[4] * 494;
- it1 %= 4096;
- /* return updated seed */
- iseed[1] = it1;
- iseed[2] = it2;
- iseed[3] = it3;
- iseed[4] = it4;
- /* convert 48-bit integer to a double number in the interval (0,1) */
- ret_val = ((real_type)it1 + ((real_type)it2 + ((real_type)it3 + (real_type)it4 * 2.44140625e-4) * 2.44140625e-4) * 2.44140625e-4) * 2.44140625e-4;
- return ret_val;
- /* End of RAN */
- } /* ran_ */
- /*
- Double results:
- ~~~~~~~~~~~~~~
- norm. resid resid machep x(1) x(n)
- 6.4915 7.207e-013 2.2204e-016 1 1
- times are reported for matrices of order 1000
- factor solve total mflops unit ratio
- times for array with leading dimension of1001
- 1.443 0.003 1.446 462.43 0.004325 25.821
- mpf_class results:
- ~~~~~~~~~~~~~~~~~~
- norm. resid resid machep x(1) x(n)
- 3.6575e-05 5.2257e-103 2.8575e-101 1 1
- times are reported for matrices of order 1000
- factor solve total mflops unit ratio
- times for array with leading dimension of1001
- 266.45 0.798 267.24 2.5021 0.79933 4772.2
- number<gmp_float<100> >:
- ~~~~~~~~~~~~~~~~~~~~~~~~~~~
- norm. resid resid machep x(1) x(n)
- 0.36575e-4 0.52257e-102 0.28575e-100 0.1e1 0.1e1
- times are reported for matrices of order 1000
- factor solve total mflops unit ratio
- times for array with leading dimension of1001
- 279.96 0.84 280.8 2.3813 0.83988 5014.3
- boost::multiprecision::ef::cpp_dec_float_50:
- ~~~~~~~~~~~~~~~~~~~~~~~~~
- norm. resid resid machep x(1) x(n)
- 2.551330735e-16 1.275665107e-112 1e-99 1 1
- times are reported for matrices of order 1000
- factor solve total mflops unit ratio
- times for array with leading dimension of1001
- 363.89 1.074 364.97 1.8321 1.0916 6517.3
- */
|