Subject: [Boost-commit] svn:boost r73649 - in sandbox/big_number: boost/math boost/math/big_number boost/math/concepts libs/math/test
From: john_at_[hidden]
Date: 2011-08-11 07:28:12
Author: johnmaddock
Date: 2011-08-11 07:28:11 EDT (Thu, 11 Aug 2011)
New Revision: 73649
URL: http://svn.boost.org/trac/boost/changeset/73649
Log:
Add rvalue reference support.
Add LINPACK benchmark test.
Update arithmetic tests to work with types other than big_number.
Fix precision of ostream& operator<<.
Added:
sandbox/big_number/libs/math/test/linpack-benchmark.cpp (contents, props changed)
Text files modified:
sandbox/big_number/boost/math/big_number.hpp | 26 +++++++++++++--------
sandbox/big_number/boost/math/big_number/gmp.hpp | 46 +++++++++++++++++++++++++++++++++------
sandbox/big_number/boost/math/concepts/big_number_architypes.hpp | 2
sandbox/big_number/libs/math/test/test_arithmetic.cpp | 34 ++++++++++++++++++++++++++--
4 files changed, 87 insertions(+), 21 deletions(-)
Modified: sandbox/big_number/boost/math/big_number.hpp
==============================================================================
--- sandbox/big_number/boost/math/big_number.hpp (original)
+++ sandbox/big_number/boost/math/big_number.hpp 2011-08-11 07:28:11 EDT (Thu, 11 Aug 2011)
@@ -190,13 +190,6 @@
BOOST_ASSERT(proto::value(*this) == this);
m_backend = canonical_value(v);
}
- /*
- big_number(unsigned digits10, typename enable_if_c<false == Backend::is_fixed>::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 @@
do_assign(e, typename proto::tag_of<Exp>::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 <class Exp>
big_number& operator+=(const detail::big_number_exp<Exp>& e)
{
@@ -367,9 +373,9 @@
//
// 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 @@
template <class Backend>
std::ostream& operator << (std::ostream& os, const big_number<Backend>& r)
{
- return os << r.str();
+ return os << r.str(static_cast<unsigned>(os.precision()));
}
template <class Backend>
Modified: sandbox/big_number/boost/math/big_number/gmp.hpp
==============================================================================
--- sandbox/big_number/boost/math/big_number/gmp.hpp (original)
+++ sandbox/big_number/boost/math/big_number/gmp.hpp 2011-08-11 07:28:11 EDT (Thu, 11 Aug 2011)
@@ -33,11 +33,25 @@
{
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<unsigned>::digits) - 1);
@@ -233,13 +247,13 @@
{
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 @@
}
~gmp_real_imp()
{
- mpf_clear(m_data);
+ if(m_data[0]._mp_d)
+ mpf_clear(m_data);
}
void negate()
{
@@ -299,13 +314,21 @@
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<detail::gmp_real_imp<digits10>*>(this) = static_cast<detail::gmp_real_imp<digits10> const&>(o);
return *this;
}
-
+#ifndef BOOST_NO_RVALUE_REFERENCES
+ gmp_real& operator=(gmp_real&& o)
+ {
+ *static_cast<detail::gmp_real_imp<digits10>*>(this) = static_cast<detail::gmp_real_imp<digits10>&&>(o);
+ return *this;
+ }
+#endif
template <class V>
gmp_real& operator=(const V& v)
{
@@ -326,6 +349,9 @@
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 @@
*static_cast<detail::gmp_real_imp<0>*>(this) = static_cast<detail::gmp_real_imp<0> const&>(o);
return *this;
}
-
+#ifndef BOOST_NO_RVALUE_REFERENCES
+ gmp_real& operator=(gmp_real&& o)
+ {
+ *static_cast<detail::gmp_real_imp<0>*>(this) = static_cast<detail::gmp_real_imp<0> &&>(o);
+ return *this;
+ }
+#endif
template <class V>
gmp_real& operator=(const V& v)
{
@@ -706,7 +738,7 @@
{
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);
Modified: sandbox/big_number/boost/math/concepts/big_number_architypes.hpp
==============================================================================
--- sandbox/big_number/boost/math/concepts/big_number_architypes.hpp (original)
+++ sandbox/big_number/boost/math/concepts/big_number_architypes.hpp 2011-08-11 07:28:11 EDT (Thu, 11 Aug 2011)
@@ -252,7 +252,7 @@
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<std::string>(m_value));
std::cout << "Converting to string (" << s << ")" << std::endl;
Added: sandbox/big_number/libs/math/test/linpack-benchmark.cpp
==============================================================================
--- (empty file)
+++ sandbox/big_number/libs/math/test/linpack-benchmark.cpp 2011-08-11 07:28:11 EDT (Thu, 11 Aug 2011)
@@ -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 <boost/lexical_cast.hpp>
+#include <iostream>
+#include <iomanip>
+#ifdef TEST_BIG_NUMBER
+#include <boost/math/big_number/gmp.hpp>
+typedef boost::math::mpf_real_100 real_type;
+#elif defined(TEST_GMPXX)
+#include <gmpxx.h>
+typedef mpf_class real_type;
+#elif defined(TEST_EF_GMP)
+#define E_FLOAT_TYPE_GMP
+#include <e_float/e_float.h>
+typedef e_float real_type;
+#define CAST_TO_RT(x) real_type(x)
+#elif defined(TEST_MATH_EF)
+#define E_FLOAT_TYPE_GMP
+#include <boost/math/bindings/e_float.hpp>
+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 <time.h>
+
+
+#if defined(TEST_EF_GMP)
+#include <functions/functions.h>
+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<mpf_real<100> >" << 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_at_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<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;
+ 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<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;
+ 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<gmp_real<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::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
+*/
Modified: sandbox/big_number/libs/math/test/test_arithmetic.cpp
==============================================================================
--- sandbox/big_number/libs/math/test/test_arithmetic.cpp (original)
+++ sandbox/big_number/libs/math/test/test_arithmetic.cpp 2011-08-11 07:28:11 EDT (Thu, 11 Aug 2011)
@@ -4,10 +4,8 @@
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_
#include <boost/detail/lightweight_test.hpp>
-#include <boost/math/concepts/big_number_architypes.hpp>
-#include <boost/math/big_number/gmp.hpp>
-#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 <boost/math/big_number/gmp.hpp>
+#endif
+#ifdef TEST_BACKEND
+#include <boost/math/concepts/big_number_architypes.hpp>
+#endif
+#ifdef TEST_E_FLOAT
+#include <boost/math/big_number.hpp>
+#include <boost/math/bindings/e_float.hpp>
+#endif
+
template <class Real>
void test_integer_ops(const boost::mpl::false_&){}
@@ -83,9 +92,12 @@
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 @@
r = frexp(v, &exp);
BOOST_TEST(r == 0.5);
BOOST_TEST(exp == -8);
+#endif
}
template <class Real, class Num>
@@ -347,8 +360,10 @@
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 @@
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 @@
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 @@
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 @@
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:
//
@@ -485,5 +510,8 @@
#ifdef TEST_MPZ
test<boost::math::mpz_int>();
#endif
+#ifdef TEST_E_FLOAT
+ test<boost::math::ef::e_float>();
+#endif
return boost::report_errors();
}