[Boost-commit] svn:boost r73649 - in sandbox/big_number: boost/math boost/math/big_number boost/math/concepts libs/math/test

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: [phone]\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();
 }