From 30fbb4d528ad0b46880086dfbd4fbec1624f4746 Mon Sep 17 00:00:00 2001 From: John Maddock Date: Thu, 11 Aug 2011 11:28:11 +0000 Subject: [PATCH] Add rvalue reference support. Add LINPACK benchmark test. Update arithmetic tests to work with types other than big_number. Fix precision of ostream& operator<<. [SVN r73649] --- include/boost/math/big_number.hpp | 26 +- include/boost/math/big_number/gmp.hpp | 46 +- .../math/concepts/big_number_architypes.hpp | 2 +- math/test/linpack-benchmark.cpp | 1258 +++++++++++++++++ math/test/test_arithmetic.cpp | 34 +- 5 files changed, 1345 insertions(+), 21 deletions(-) create mode 100644 math/test/linpack-benchmark.cpp diff --git a/include/boost/math/big_number.hpp b/include/boost/math/big_number.hpp index af00ec74..f9a7c98b 100644 --- a/include/boost/math/big_number.hpp +++ b/include/boost/math/big_number.hpp @@ -190,13 +190,6 @@ public: BOOST_ASSERT(proto::value(*this) == this); m_backend = canonical_value(v); } - /* - big_number(unsigned digits10, typename enable_if_c::type* dummy = 0) - : base_type(digits10) - { - proto::value(*this) = this; - BOOST_ASSERT(proto::value(*this) == this); - }*/ big_number(const big_number& e, unsigned digits10) : m_backend(e.m_backend, digits10) { proto::value(*this) = this; @@ -240,6 +233,19 @@ public: do_assign(e, typename proto::tag_of::type()); } +#ifndef BOOST_NO_RVALUE_REFERENCES + big_number(big_number&& r) : m_backend(r.m_backend) + { + proto::value(*this) = this; + BOOST_ASSERT(proto::value(*this) == this); + } + big_number& operator=(big_number&& r) + { + m_backend.swap(r.m_backend); + return *this; + } +#endif + template big_number& operator+=(const detail::big_number_exp& e) { @@ -367,9 +373,9 @@ public: // // String conversion functions: // - std::string str()const + std::string str(unsigned digits = 0)const { - return m_backend.str(); + return m_backend.str(digits); } // // Default precision: @@ -1160,7 +1166,7 @@ inline typename boost::enable_if, bool>: template std::ostream& operator << (std::ostream& os, const big_number& r) { - return os << r.str(); + return os << r.str(static_cast(os.precision())); } template diff --git a/include/boost/math/big_number/gmp.hpp b/include/boost/math/big_number/gmp.hpp index aaa5a068..cba6ed56 100644 --- a/include/boost/math/big_number/gmp.hpp +++ b/include/boost/math/big_number/gmp.hpp @@ -33,11 +33,25 @@ struct gmp_real_imp { mpf_init_set(m_data, o.m_data); } +#ifndef BOOST_NO_RVALUE_REFERENCES + gmp_real_imp(gmp_real_imp&& o) + { + m_data[0] = o.m_data[0]; + o.m_data[0]._mp_d = 0; + } +#endif gmp_real_imp& operator = (const gmp_real_imp& o) { mpf_set(m_data, o.m_data); return *this; } +#ifndef BOOST_NO_RVALUE_REFERENCES + gmp_real_imp& operator = (gmp_real_imp&& o) + { + mpf_swap(m_data, o.m_data); + return *this; + } +#endif gmp_real_imp& operator = (boost::uintmax_t i) { boost::uintmax_t mask = ((1uLL << std::numeric_limits::digits) - 1); @@ -233,13 +247,13 @@ struct gmp_real_imp { mpf_swap(m_data, o.m_data); } - std::string str()const + std::string str(unsigned digits)const { mp_exp_t e; void *(*alloc_func_ptr) (size_t); void *(*realloc_func_ptr) (void *, size_t, size_t); void (*free_func_ptr) (void *, size_t); - const char* ps = mpf_get_str (0, &e, 10, 0, m_data); + const char* ps = mpf_get_str (0, &e, 10, digits, m_data); std::string s("0."); if(ps[0] == '-') { @@ -258,7 +272,8 @@ struct gmp_real_imp } ~gmp_real_imp() { - mpf_clear(m_data); + if(m_data[0]._mp_d) + mpf_clear(m_data); } void negate() { @@ -299,13 +314,21 @@ struct gmp_real : public detail::gmp_real_imp mpf_init2(this->m_data, ((digits10 + 1) * 1000L) / 301L); } gmp_real(const gmp_real& o) : gmp_real_imp(o) {} - +#ifndef BOOST_NO_RVALUE_REFERENCES + gmp_real(gmp_real&& o) : gmp_real_imp(o) {} +#endif gmp_real& operator=(const gmp_real& o) { *static_cast*>(this) = static_cast const&>(o); return *this; } - +#ifndef BOOST_NO_RVALUE_REFERENCES + gmp_real& operator=(gmp_real&& o) + { + *static_cast*>(this) = static_cast&&>(o); + return *this; + } +#endif template gmp_real& operator=(const V& v) { @@ -326,6 +349,9 @@ struct gmp_real<0> : public detail::gmp_real_imp<0> mpf_init2(this->m_data, ((digits10 + 1) * 1000L) / 301L); } gmp_real(const gmp_real& o) : gmp_real_imp(o) {} +#ifndef BOOST_NO_RVALUE_REFERENCES + gmp_real(gmp_real&& o) : gmp_real_imp(o) {} +#endif gmp_real(const gmp_real& o, unsigned digits10) { mpf_init2(this->m_data, ((digits10 + 1) * 1000L) / 301L); @@ -337,7 +363,13 @@ struct gmp_real<0> : public detail::gmp_real_imp<0> *static_cast*>(this) = static_cast const&>(o); return *this; } - +#ifndef BOOST_NO_RVALUE_REFERENCES + gmp_real& operator=(gmp_real&& o) + { + *static_cast*>(this) = static_cast &&>(o); + return *this; + } +#endif template gmp_real& operator=(const V& v) { @@ -706,7 +738,7 @@ struct gmp_int { mpz_swap(m_data, o.m_data); } - std::string str()const + std::string str(unsigned)const { void *(*alloc_func_ptr) (size_t); void *(*realloc_func_ptr) (void *, size_t, size_t); diff --git a/include/boost/math/concepts/big_number_architypes.hpp b/include/boost/math/concepts/big_number_architypes.hpp index f7c196e1..31b4d992 100644 --- a/include/boost/math/concepts/big_number_architypes.hpp +++ b/include/boost/math/concepts/big_number_architypes.hpp @@ -252,7 +252,7 @@ struct big_number_backend_real_architype std::cout << "Swapping (" << m_value << " with " << o.m_value << ")" << std::endl; std::swap(m_value, o.m_value); } - std::string str()const + std::string str(unsigned)const { std::string s(boost::lexical_cast(m_value)); std::cout << "Converting to string (" << s << ")" << std::endl; diff --git a/math/test/linpack-benchmark.cpp b/math/test/linpack-benchmark.cpp new file mode 100644 index 00000000..c2c1d054 --- /dev/null +++ b/math/test/linpack-benchmark.cpp @@ -0,0 +1,1258 @@ +/////////////////////////////////////////////////////////////////////////////// +// 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 +#include +#include +#ifdef TEST_BIG_NUMBER +#include +typedef boost::math::mpf_real_100 real_type; +#elif defined(TEST_GMPXX) +#include +typedef mpf_class real_type; +#elif defined(TEST_EF_GMP) +#define E_FLOAT_TYPE_GMP +#include +typedef e_float real_type; +#define CAST_TO_RT(x) real_type(x) +#elif defined(TEST_MATH_EF) +#define E_FLOAT_TYPE_GMP +#include +typedef boost::math::ef::e_float real_type; +//#define CAST_TO_RT(x) real_type(x) +#else +typedef double real_type; +#endif + +#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 min +#undef max +#undef dmin +#undef dmax +#define dmin min +#define dmax max + +} +#include + + +#if defined(TEST_EF_GMP) +#include +using namespace ef; +#endif + +using std::min; +using std::max; + +/* 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 main() +{ +#ifdef TEST_BIG_NUMBER + std::cout << "Testing big_number >" << std::endl; +#elif defined(TEST_GMPXX) + std::cout << "Testing mpfr_class at 100 decimal degits" << std::endl; + mpf_set_default_prec(((100 + 1) * 1000L) / 301L); +#elif defined(TEST_EF_GMP) + std::cout << "Testing gmp::e_float" << std::endl; +#elif defined(TEST_MATH_EF) + std::cout << "Testing boost::math::ef::e_float" << 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]; + extern /* Subroutine */ int dgefa_(real_type *, integer *, integer *, + integer *, integer *), dgesl_(real_type *, integer *, integer *, + integer *, real_type *, integer *); + static real_type resid, norma; + extern /* Subroutine */ int dmxpy_(integer *, real_type *, integer *, + integer *, real_type *, real_type *); + static real_type normx; + extern /* Subroutine */ int matgen_(real_type *, integer *, integer *, + real_type *, real_type *); + static real_type residn; + extern real_type epslon_(real_type *); + + /* 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(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; + extern real_type ran_(integer *); + 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; + extern /* Subroutine */ int dscal_(integer *, real_type *, real_type *, + integer *), daxpy_(integer *, real_type *, real_type *, integer + *, real_type *, integer *); + extern integer idamax_(integer *, real_type *, integer *); + + + /* 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; + extern real_type ddot_(integer *, real_type *, integer *, real_type *, + integer *); + extern /* Subroutine */ int daxpy_(integer *, real_type *, real_type *, + integer *, real_type *, integer *); + + + /* 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) +{ +#ifdef TEST_BIG_NUMBER + return ldexp(1.0, 1 - ((100 + 1) * 1000L) / 301L); +#elif defined(TEST_GMPXX) + return ldexp(1.0, 1 - ((100 + 1) * 1000L) / 301L); +#else + return CAST_TO_RT(std::numeric_limits::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; + extern /* Subroutine */ int dmxpy_(integer *, real_type *, integer *, + integer *, real_type *, real_type *); + + + /* 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 + + +big_number >: +~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + 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::math::ef::e_float: +~~~~~~~~~~~~~~~~~~~~~~~~~ + + 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 +*/ diff --git a/math/test/test_arithmetic.cpp b/math/test/test_arithmetic.cpp index fa2fc9c6..99436325 100644 --- a/math/test/test_arithmetic.cpp +++ b/math/test/test_arithmetic.cpp @@ -4,10 +4,8 @@ // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_ #include -#include -#include -#if !defined(TEST_MPF50) && !defined(TEST_MPF) && !defined(TEST_BACKEND) && !defined(TEST_MPZ) +#if !defined(TEST_MPF50) && !defined(TEST_MPF) && !defined(TEST_BACKEND) && !defined(TEST_MPZ) && !defined(TEST_E_FLOAT) # define TEST_MPF50 # define TEST_MPF # define TEST_BACKEND @@ -22,6 +20,17 @@ #endif +#if defined(TEST_MPF50) || defined(TEST_MPF) || defined(TEST_MPZ) +#include +#endif +#ifdef TEST_BACKEND +#include +#endif +#ifdef TEST_E_FLOAT +#include +#include +#endif + template void test_integer_ops(const boost::mpl::false_&){} @@ -83,9 +92,12 @@ void test_real_ops(const boost::mpl::true_&) BOOST_TEST(ceil(Real(5) / 2) == 3); BOOST_TEST(floor(Real(-5) / 2) == -3); BOOST_TEST(ceil(Real(-5) / 2) == -2); +#ifndef TEST_E_FLOAT BOOST_TEST(trunc(Real(5) / 2) == 2); BOOST_TEST(trunc(Real(-5) / 2) == -2); +#endif +#ifndef TEST_E_FLOAT // // ldexp and frexp, these pretty much have to implemented by each backend: // @@ -100,6 +112,7 @@ void test_real_ops(const boost::mpl::true_&) r = frexp(v, &exp); BOOST_TEST(r == 0.5); BOOST_TEST(exp == -8); +#endif } template @@ -347,8 +360,10 @@ void test() ac = a * ac; BOOST_TEST(ac == 8*8); ac = a; +#ifndef TEST_E_FLOAT ac = ac + "8"; BOOST_TEST(ac == 16); +#endif ac = a; ac += +a; BOOST_TEST(ac == 16); @@ -362,8 +377,10 @@ void test() ac += b*c; BOOST_TEST(ac == 8 + 64 * 500); ac = a; +#ifndef TEST_E_FLOAT ac = ac - "8"; BOOST_TEST(ac == 0); +#endif ac = a; ac -= +a; BOOST_TEST(ac == 0); @@ -382,8 +399,12 @@ void test() ac = a; ac -= ac * b; BOOST_TEST(ac == 8 - 8 * 64); +#ifndef TEST_E_FLOAT ac = a * "8"; BOOST_TEST(ac == 64); +#else + ac = a * 8; +#endif ac *= +a; BOOST_TEST(ac == 64 * 8); ac = a; @@ -398,8 +419,10 @@ void test() ac = a; ac *= b + c; BOOST_TEST(ac == 8 * (64 + 500)); +#ifndef TEST_E_FLOAT ac = b / "8"; BOOST_TEST(ac == 8); +#endif ac = b; ac /= +a; BOOST_TEST(ac == 8); @@ -412,8 +435,10 @@ void test() ac = b; ac /= a + Real(0); BOOST_TEST(ac == 8); +#ifndef TEST_E_FLOAT ac = a + std::string("8"); BOOST_TEST(ac == 16); +#endif // // Comparisons: // @@ -484,6 +509,9 @@ int main() #endif #ifdef TEST_MPZ test(); +#endif +#ifdef TEST_E_FLOAT + test(); #endif return boost::report_errors(); }