///////////////////////////////////////////////////////////////////////////////
//  Copyright 2011 John Maddock. Distributed under the Boost
//  Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

/* 1000d.f -- translated by f2c (version 20050501).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

http://www.netlib.org/f2c/libf2c.zip
*/
#include <iostream>
#include <iomanip>
#include <cmath>

#if defined(TEST_GMPXX)
#include <gmpxx.h>
typedef mpf_class real_type;
#elif defined(TEST_MPFRXX)
#include <gmpfrxx.h>
typedef mpfr_class real_type;
#elif defined(TEST_CPP_DEC_FLOAT)
#include <boost/multiprecision/cpp_dec_float.hpp>
typedef boost::multiprecision::cpp_dec_float_50 real_type;
#elif defined(TEST_MPFR_50)
#include <boost/multiprecision/mpfr.hpp>
typedef boost::multiprecision::mpfr_float_50 real_type;
#elif defined(TEST_MPF_50)
#include <boost/multiprecision/gmp.hpp>
typedef boost::multiprecision::mpf_float_50 real_type;
#elif defined(NATIVE_FLOAT128)
#include <boost/multiprecision/float128.hpp>
typedef __float128 real_type;

std::ostream& operator<<(std::ostream& os, const __float128& f)
{
   return os << boost::multiprecision::float128(f);
}

#include <boost/type_traits/has_left_shift.hpp>

namespace boost{

template<>
struct has_left_shift<std::basic_ostream<char>, __float128> : public mpl::true_ {};

template<>
double lexical_cast<double, __float128>(const __float128& f)
{ return f; }

}

#elif defined(TEST_FLOAT128)
#include <boost/multiprecision/float128.hpp>
typedef boost::multiprecision::float128 real_type;
#elif defined(TEST_CPP_BIN_FLOAT_QUAD)
#include <boost/multiprecision/cpp_bin_float.hpp>
typedef boost::multiprecision::cpp_bin_float_quad real_type;
#else
typedef double real_type;
#endif

#include <boost/lexical_cast.hpp>

#ifndef CAST_TO_RT
#  define CAST_TO_RT(x) x
#endif

extern "C" {
#include "f2c.h"
   integer s_wsfe(cilist *), e_wsfe(void), do_fio(integer *, char *, ftnlen),
      s_wsle(cilist *), do_lio(integer *, integer *, char *, ftnlen), 
      e_wsle(void);
   /* Subroutine */ int s_stop(char *, ftnlen);

#undef abs
#undef dabs
#define dabs abs
#undef dmin
#undef dmax
#define dmin min
#define dmax max

}
#include <time.h>

using std::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 dgefa_(real_type *, integer *, integer *, integer *, integer *), dgesl_(real_type *, integer *, integer *, integer *, real_type *, integer *);
int dmxpy_(integer *, real_type *, integer *, integer *, real_type *, real_type *);
int matgen_(real_type *, integer *, integer *, real_type *, real_type *);
real_type epslon_(real_type *);
real_type ran_(integer *);
int dscal_(integer *, real_type *, real_type *, integer *);
int daxpy_(integer *, real_type *, real_type *, integer *, real_type *, integer *);
integer idamax_(integer *, real_type *, integer *);
real_type ddot_(integer *, real_type *, integer *, real_type *, integer *);
int daxpy_(integer *, real_type *, real_type *, integer *, real_type *, integer *);
int dmxpy_(integer *, real_type *, integer *, integer *, real_type *, real_type *);

extern "C" int MAIN__()
{
#ifdef TEST_MPF_50
   std::cout << "Testing number<mpf_float<50> >" << std::endl;
#elif defined(TEST_MPFR_50)
   std::cout << "Testing number<mpf_float<50> >" << std::endl;
#elif defined(TEST_GMPXX)
   std::cout << "Testing mpf_class at 50 decimal degits" << std::endl;
   mpf_set_default_prec(((50 + 1) * 1000L) / 301L);
#elif defined(TEST_MPFRXX)
   std::cout << "Testing mpfr_class at 50 decimal degits" << std::endl;
   mpfr_set_default_prec(((50 + 1) * 1000L) / 301L);
#elif defined(TEST_CPP_DEC_FLOAT)
   std::cout << "Testing number<cpp_dec_float<50> >" << std::endl;
#elif defined(NATIVE_FLOAT128)
   std::cout << "Testing __float128" << std::endl;
#elif defined(TEST_FLOAT128)
   std::cout << "Testing number<float128_backend, et_off>" << std::endl;
#else
   std::cout << "Testing double" << std::endl;
#endif


   /* Format strings */
   static char fmt_1[] = "(\002 Please send the results of this run to:\002"
      "//\002 Jack J. Dongarra\002/\002 Computer Science Department\002/"
      "\002 University of Tennessee\002/\002 Knoxville, Tennessee 37996"
      "-1300\002//\002 Fax: 615-974-8296\002//\002 Internet: dongarra@c"
      "s.utk.edu\002/)";
   static char fmt_40[] = "(\002     norm. resid      resid           mac"
      "hep\002,\002         x(1)          x(n)\002)";
   static char fmt_50[] = "(1p5e16.8)";
   static char fmt_60[] = "(//\002    times are reported for matrices of or"
      "der \002,i5)";
   static char fmt_70[] = "(6x,\002factor\002,5x,\002solve\002,6x,\002tota"
      "l\002,5x,\002mflops\002,7x,\002unit\002,6x,\002ratio\002)";
   static char fmt_80[] = "(\002 times for array with leading dimension o"
      "f\002,i4)";
   static char fmt_110[] = "(6(1pe11.3))";

   /* System generated locals */
   integer i__1;
   real_type d__1, d__2, d__3;

   /* Builtin functions */

   /* Local variables */
   static real_type a[1001000] /* was [1001][1000] */, b[1000];
   static integer i__, n;
   static real_type x[1000];
   static double t1;
   static integer lda;
   static double ops;
   static real_type eps;
   static integer info;
   static double time[6], cray, total;
   static integer ipvt[1000];
   static real_type resid, norma;
   static real_type normx;
   static real_type residn;

   /* Fortran I/O blocks */
   static cilist io___4 = { 0, 6, 0, fmt_1, 0 };
   static cilist io___20 = { 0, 6, 0, fmt_40, 0 };
   static cilist io___21 = { 0, 6, 0, fmt_50, 0 };
   static cilist io___22 = { 0, 6, 0, fmt_60, 0 };
   static cilist io___23 = { 0, 6, 0, fmt_70, 0 };
   static cilist io___24 = { 0, 6, 0, fmt_80, 0 };
   static cilist io___25 = { 0, 6, 0, fmt_110, 0 };
   static cilist io___26 = { 0, 6, 0, 0, 0 };


   lda = 1001;

   /*     this program was updated on 10/12/92 to correct a */
   /*     problem with the random number generator. The previous */
   /*     random number generator had a short period and produced */
   /*     singular matrices occasionally. */

   n = 1000;
   cray = .056f;
   s_wsfe(&io___4);
   e_wsfe();
   /* Computing 3rd power */
   d__1 = (real_type) n;
   /* Computing 2nd power */
   d__2 = (real_type) n;
   ops = boost::lexical_cast<double>(real_type(d__1 * (d__1 * d__1) * 2. / 3. + d__2 * d__2 * 2.));

   matgen_(a, &lda, &n, b, &norma);

   /* ****************************************************************** */
   /* ****************************************************************** */
   /*        you should replace the call to dgefa and dgesl */
   /*        by calls to your linear equation solver. */
   /* ****************************************************************** */
   /* ****************************************************************** */

   t1 = second_();
   dgefa_(a, &lda, &n, ipvt, &info);
   time[0] = second_() - t1;
   t1 = second_();
   dgesl_(a, &lda, &n, ipvt, b, &c__0);
   time[1] = second_() - t1;
   total = time[0] + time[1];
   /* ****************************************************************** */
   /* ****************************************************************** */

   /*     compute a residual to verify results. */

   i__1 = n;
   for (i__ = 1; i__ <= i__1; ++i__) {
      x[i__ - 1] = b[i__ - 1];
      /* L10: */
   }
   matgen_(a, &lda, &n, b, &norma);
   i__1 = n;
   for (i__ = 1; i__ <= i__1; ++i__) {
      b[i__ - 1] = -b[i__ - 1];
      /* L20: */
   }
   dmxpy_(&n, b, &n, &lda, x, a);
   resid = CAST_TO_RT(0);
   normx = CAST_TO_RT(0);
   i__1 = n;
   for (i__ = 1; i__ <= i__1; ++i__) {
      /* Computing MAX */
      d__2 = resid, d__3 = (d__1 = b[i__ - 1], abs(d__1));
      resid = (max)(d__2,d__3);
      /* Computing MAX */
      d__2 = normx, d__3 = (d__1 = x[i__ - 1], abs(d__1));
      normx = (max)(d__2,d__3);
      /* L30: */
   }
   eps = epslon_(&c_b7);
   residn = resid / (n * norma * normx * eps);
   s_wsfe(&io___20);
   e_wsfe();
   s_wsfe(&io___21);
   /*
   do_fio(&c__1, (char *)&residn, (ftnlen)sizeof(real_type));
   do_fio(&c__1, (char *)&resid, (ftnlen)sizeof(real_type));
   do_fio(&c__1, (char *)&eps, (ftnlen)sizeof(real_type));
   do_fio(&c__1, (char *)&x[0], (ftnlen)sizeof(real_type));
   do_fio(&c__1, (char *)&x[n - 1], (ftnlen)sizeof(real_type));
   */
   std::cout << std::setw(12) << std::setprecision(5) << residn << " " << resid << " " << eps << " " << x[0] << " " << x[n-1] << std::endl;
   e_wsfe();

   s_wsfe(&io___22);
   do_fio(&c__1, (char *)&n, (ftnlen)sizeof(integer));
   e_wsfe();
   s_wsfe(&io___23);
   e_wsfe();

   time[2] = total;
   time[3] = ops / (total * 1e6);
   time[4] = 2. / time[3];
   time[5] = total / cray;
   s_wsfe(&io___24);
   do_fio(&c__1, (char *)&lda, (ftnlen)sizeof(integer));
   e_wsfe();
   s_wsfe(&io___25);
   for (i__ = 1; i__ <= 6; ++i__) {
      // do_fio(&c__1, (char *)&time[i__ - 1], (ftnlen)sizeof(real_type));
      std::cout << std::setw(12) << std::setprecision(5) << time[i__ - 1];
   }
   e_wsfe();
   s_wsle(&io___26);
   do_lio(&c__9, &c__1, " end of tests -- this version dated 10/12/92", (
      ftnlen)44);
   e_wsle();

   s_stop("", (ftnlen)0);
   return 0;
} /* MAIN__ */

/* Subroutine */ int matgen_(real_type *a, integer *lda, integer *n, 
   real_type *b, real_type *norma)
{
   /* System generated locals */
   integer a_dim1, a_offset, i__1, i__2;
   real_type d__1, d__2;

   /* Local variables */
   static integer i__, j;
   static integer init[4];


   /* Parameter adjustments */
   a_dim1 = *lda;
   a_offset = 1 + a_dim1;
   a -= a_offset;
   --b;

   /* Function Body */
   init[0] = 1;
   init[1] = 2;
   init[2] = 3;
   init[3] = 1325;
   *norma = CAST_TO_RT(0);
   i__1 = *n;
   for (j = 1; j <= i__1; ++j) {
      i__2 = *n;
      for (i__ = 1; i__ <= i__2; ++i__) {
         a[i__ + j * a_dim1] = ran_(init) - .5f;
         /* Computing MAX */
         d__2 = (d__1 = a[i__ + j * a_dim1], abs(d__1));
         *norma = (max)(d__2,*norma);
         /* L20: */
      }
      /* L30: */
   }
   i__1 = *n;
   for (i__ = 1; i__ <= i__1; ++i__) {
      b[i__] = CAST_TO_RT(0);
      /* L35: */
   }
   i__1 = *n;
   for (j = 1; j <= i__1; ++j) {
      i__2 = *n;
      for (i__ = 1; i__ <= i__2; ++i__) {
         b[i__] += a[i__ + j * a_dim1];
         /* L40: */
      }
      /* L50: */
   }
   return 0;
} /* matgen_ */

/* Subroutine */ int dgefa_(real_type *a, integer *lda, integer *n, integer *
   ipvt, integer *info)
{
   /* System generated locals */
   integer a_dim1, a_offset, i__1, i__2, i__3;

   /* Local variables */
   static integer j, k, l;
   static real_type t;
   static integer kp1, nm1;


   /*     dgefa factors a double precision matrix by gaussian elimination. */

   /*     dgefa is usually called by dgeco, but it can be called */
   /*     directly with a saving in time if  rcond  is not needed. */
   /*     (time for dgeco) = (1 + 9/n)*(time for dgefa) . */

   /*     on entry */

   /*        a       double precision(lda, n) */
   /*                the matrix to be factored. */

   /*        lda     integer */
   /*                the leading dimension of the array  a . */

   /*        n       integer */
   /*                the order of the matrix  a . */

   /*     on return */

   /*        a       an upper triangular matrix and the multipliers */
   /*                which were used to obtain it. */
   /*                the factorization can be written  a = l*u  where */
   /*                l  is a product of permutation and unit lower */
   /*                triangular matrices and  u  is upper triangular. */

   /*        ipvt    integer(n) */
   /*                an integer vector of pivot indices. */

   /*        info    integer */
   /*                = 0  normal value. */
   /*                = k  if  u(k,k) .eq. 0.0 .  this is not an error */
   /*                     condition for this subroutine, but it does */
   /*                     indicate that dgesl or dgedi will divide by zero */
   /*                     if called.  use  rcond  in dgeco for a reliable */
   /*                     indication of singularity. */

   /*     linpack. this version dated 08/14/78 . */
   /*     cleve moler, university of new mexico, argonne national lab. */

   /*     subroutines and functions */

   /*     blas daxpy,dscal,idamax */

   /*     internal variables */



   /*     gaussian elimination with partial pivoting */

   /* Parameter adjustments */
   a_dim1 = *lda;
   a_offset = 1 + a_dim1;
   a -= a_offset;
   --ipvt;

   /* Function Body */
   *info = 0;
   nm1 = *n - 1;
   if (nm1 < 1) {
      goto L70;
   }
   i__1 = nm1;
   for (k = 1; k <= i__1; ++k) {
      kp1 = k + 1;

      /*        find l = pivot index */

      i__2 = *n - k + 1;
      l = idamax_(&i__2, &a[k + k * a_dim1], &c__1) + k - 1;
      ipvt[k] = l;

      /*        zero pivot implies this column already triangularized */

      if (a[l + k * a_dim1] == 0.) {
         goto L40;
      }

      /*           interchange if necessary */

      if (l == k) {
         goto L10;
      }
      t = a[l + k * a_dim1];
      a[l + k * a_dim1] = a[k + k * a_dim1];
      a[k + k * a_dim1] = t;
L10:

      /*           compute multipliers */

      t = -1. / a[k + k * a_dim1];
      i__2 = *n - k;
      dscal_(&i__2, &t, &a[k + 1 + k * a_dim1], &c__1);

      /*           row elimination with column indexing */

      i__2 = *n;
      for (j = kp1; j <= i__2; ++j) {
         t = a[l + j * a_dim1];
         if (l == k) {
            goto L20;
         }
         a[l + j * a_dim1] = a[k + j * a_dim1];
         a[k + j * a_dim1] = t;
L20:
         i__3 = *n - k;
         daxpy_(&i__3, &t, &a[k + 1 + k * a_dim1], &c__1, &a[k + 1 + j * 
            a_dim1], &c__1);
         /* L30: */
      }
      goto L50;
L40:
      *info = k;
L50:
      /* L60: */
      ;
   }
L70:
   ipvt[*n] = *n;
   if (a[*n + *n * a_dim1] == 0.) {
      *info = *n;
   }
   return 0;
} /* dgefa_ */

/* Subroutine */ int dgesl_(real_type *a, integer *lda, integer *n, integer *
   ipvt, real_type *b, integer *job)
{
   /* System generated locals */
   integer a_dim1, a_offset, i__1, i__2;

   /* Local variables */
   static integer k, l;
   static real_type t;
   static integer kb, nm1;


   /*     dgesl solves the double precision system */
   /*     a * x = b  or  trans(a) * x = b */
   /*     using the factors computed by dgeco or dgefa. */

   /*     on entry */

   /*        a       double precision(lda, n) */
   /*                the output from dgeco or dgefa. */

   /*        lda     integer */
   /*                the leading dimension of the array  a . */

   /*        n       integer */
   /*                the order of the matrix  a . */

   /*        ipvt    integer(n) */
   /*                the pivot vector from dgeco or dgefa. */

   /*        b       double precision(n) */
   /*                the right hand side vector. */

   /*        job     integer */
   /*                = 0         to solve  a*x = b , */
   /*                = nonzero   to solve  trans(a)*x = b  where */
   /*                            trans(a)  is the transpose. */

   /*     on return */

   /*        b       the solution vector  x . */

   /*     error condition */

   /*        a division by zero will occur if the input factor contains a */
   /*        zero on the diagonal.  technically this indicates singularity */
   /*        but it is often caused by improper arguments or improper */
   /*        setting of lda .  it will not occur if the subroutines are */
   /*        called correctly and if dgeco has set rcond .gt. 0.0 */
   /*        or dgefa has set info .eq. 0 . */

   /*     to compute  inverse(a) * c  where  c  is a matrix */
   /*     with  p  columns */
   /*           call dgeco(a,lda,n,ipvt,rcond,z) */
   /*           if (rcond is too small) go to ... */
   /*           do 10 j = 1, p */
   /*              call dgesl(a,lda,n,ipvt,c(1,j),0) */
   /*        10 continue */

   /*     linpack. this version dated 08/14/78 . */
   /*     cleve moler, university of new mexico, argonne national lab. */

   /*     subroutines and functions */

   /*     blas daxpy,ddot */

   /*     internal variables */


   /* Parameter adjustments */
   a_dim1 = *lda;
   a_offset = 1 + a_dim1;
   a -= a_offset;
   --ipvt;
   --b;

   /* Function Body */
   nm1 = *n - 1;
   if (*job != 0) {
      goto L50;
   }

   /*        job = 0 , solve  a * x = b */
   /*        first solve  l*y = b */

   if (nm1 < 1) {
      goto L30;
   }
   i__1 = nm1;
   for (k = 1; k <= i__1; ++k) {
      l = ipvt[k];
      t = b[l];
      if (l == k) {
         goto L10;
      }
      b[l] = b[k];
      b[k] = t;
L10:
      i__2 = *n - k;
      daxpy_(&i__2, &t, &a[k + 1 + k * a_dim1], &c__1, &b[k + 1], &c__1);
      /* L20: */
   }
L30:

   /*        now solve  u*x = y */

   i__1 = *n;
   for (kb = 1; kb <= i__1; ++kb) {
      k = *n + 1 - kb;
      b[k] /= a[k + k * a_dim1];
      t = -b[k];
      i__2 = k - 1;
      daxpy_(&i__2, &t, &a[k * a_dim1 + 1], &c__1, &b[1], &c__1);
      /* L40: */
   }
   goto L100;
L50:

   /*        job = nonzero, solve  trans(a) * x = b */
   /*        first solve  trans(u)*y = b */

   i__1 = *n;
   for (k = 1; k <= i__1; ++k) {
      i__2 = k - 1;
      t = ddot_(&i__2, &a[k * a_dim1 + 1], &c__1, &b[1], &c__1);
      b[k] = (b[k] - t) / a[k + k * a_dim1];
      /* L60: */
   }

   /*        now solve trans(l)*x = y */

   if (nm1 < 1) {
      goto L90;
   }
   i__1 = nm1;
   for (kb = 1; kb <= i__1; ++kb) {
      k = *n - kb;
      i__2 = *n - k;
      b[k] += ddot_(&i__2, &a[k + 1 + k * a_dim1], &c__1, &b[k + 1], &c__1);
      l = ipvt[k];
      if (l == k) {
         goto L70;
      }
      t = b[l];
      b[l] = b[k];
      b[k] = t;
L70:
      /* L80: */
      ;
   }
L90:
L100:
   return 0;
} /* dgesl_ */

/* Subroutine */ int daxpy_(integer *n, real_type *da, real_type *dx, 
   integer *incx, real_type *dy, integer *incy)
{
   /* System generated locals */
   integer i__1;

   /* Local variables */
   static integer i__, m, ix, iy, mp1;


   /*     constant times a vector plus a vector. */
   /*     uses unrolled loops for increments equal to one. */
   /*     jack dongarra, linpack, 3/11/78. */


   /* Parameter adjustments */
   --dy;
   --dx;

   /* Function Body */
   if (*n <= 0) {
      return 0;
   }
   if (*da == 0.) {
      return 0;
   }
   if (*incx == 1 && *incy == 1) {
      goto L20;
   }

   /*        code for unequal increments or equal increments */
   /*          not equal to 1 */

   ix = 1;
   iy = 1;
   if (*incx < 0) {
      ix = (-(*n) + 1) * *incx + 1;
   }
   if (*incy < 0) {
      iy = (-(*n) + 1) * *incy + 1;
   }
   i__1 = *n;
   for (i__ = 1; i__ <= i__1; ++i__) {
      dy[iy] += *da * dx[ix];
      ix += *incx;
      iy += *incy;
      /* L10: */
   }
   return 0;

   /*        code for both increments equal to 1 */


   /*        clean-up loop */

L20:
   m = *n % 4;
   if (m == 0) {
      goto L40;
   }
   i__1 = m;
   for (i__ = 1; i__ <= i__1; ++i__) {
      dy[i__] += *da * dx[i__];
      /* L30: */
   }
   if (*n < 4) {
      return 0;
   }
L40:
   mp1 = m + 1;
   i__1 = *n;
   for (i__ = mp1; i__ <= i__1; i__ += 4) {
      dy[i__] += *da * dx[i__];
      dy[i__ + 1] += *da * dx[i__ + 1];
      dy[i__ + 2] += *da * dx[i__ + 2];
      dy[i__ + 3] += *da * dx[i__ + 3];
      /* L50: */
   }
   return 0;
} /* daxpy_ */

real_type ddot_(integer *n, real_type *dx, integer *incx, real_type *dy, 
   integer *incy)
{
   /* System generated locals */
   integer i__1;
   real_type ret_val;

   /* Local variables */
   static integer i__, m, ix, iy, mp1;
   static real_type dtemp;


   /*     forms the dot product of two vectors. */
   /*     uses unrolled loops for increments equal to one. */
   /*     jack dongarra, linpack, 3/11/78. */


   /* Parameter adjustments */
   --dy;
   --dx;

   /* Function Body */
   ret_val = CAST_TO_RT(0);
   dtemp = CAST_TO_RT(0);
   if (*n <= 0) {
      return ret_val;
   }
   if (*incx == 1 && *incy == 1) {
      goto L20;
   }

   /*        code for unequal increments or equal increments */
   /*          not equal to 1 */

   ix = 1;
   iy = 1;
   if (*incx < 0) {
      ix = (-(*n) + 1) * *incx + 1;
   }
   if (*incy < 0) {
      iy = (-(*n) + 1) * *incy + 1;
   }
   i__1 = *n;
   for (i__ = 1; i__ <= i__1; ++i__) {
      dtemp += dx[ix] * dy[iy];
      ix += *incx;
      iy += *incy;
      /* L10: */
   }
   ret_val = dtemp;
   return ret_val;

   /*        code for both increments equal to 1 */


   /*        clean-up loop */

L20:
   m = *n % 5;
   if (m == 0) {
      goto L40;
   }
   i__1 = m;
   for (i__ = 1; i__ <= i__1; ++i__) {
      dtemp += dx[i__] * dy[i__];
      /* L30: */
   }
   if (*n < 5) {
      goto L60;
   }
L40:
   mp1 = m + 1;
   i__1 = *n;
   for (i__ = mp1; i__ <= i__1; i__ += 5) {
      dtemp = dtemp + dx[i__] * dy[i__] + dx[i__ + 1] * dy[i__ + 1] + dx[
         i__ + 2] * dy[i__ + 2] + dx[i__ + 3] * dy[i__ + 3] + dx[i__ + 
            4] * dy[i__ + 4];
         /* L50: */
   }
L60:
   ret_val = dtemp;
   return ret_val;
} /* ddot_ */

/* Subroutine */ int dscal_(integer *n, real_type *da, real_type *dx, 
   integer *incx)
{
   /* System generated locals */
   integer i__1, i__2;

   /* Local variables */
   static integer i__, m, mp1, nincx;


   /*     scales a vector by a constant. */
   /*     uses unrolled loops for increment equal to one. */
   /*     jack dongarra, linpack, 3/11/78. */


   /* Parameter adjustments */
   --dx;

   /* Function Body */
   if (*n <= 0) {
      return 0;
   }
   if (*incx == 1) {
      goto L20;
   }

   /*        code for increment not equal to 1 */

   nincx = *n * *incx;
   i__1 = nincx;
   i__2 = *incx;
   for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
      dx[i__] = *da * dx[i__];
      /* L10: */
   }
   return 0;

   /*        code for increment equal to 1 */


   /*        clean-up loop */

L20:
   m = *n % 5;
   if (m == 0) {
      goto L40;
   }
   i__2 = m;
   for (i__ = 1; i__ <= i__2; ++i__) {
      dx[i__] = *da * dx[i__];
      /* L30: */
   }
   if (*n < 5) {
      return 0;
   }
L40:
   mp1 = m + 1;
   i__2 = *n;
   for (i__ = mp1; i__ <= i__2; i__ += 5) {
      dx[i__] = *da * dx[i__];
      dx[i__ + 1] = *da * dx[i__ + 1];
      dx[i__ + 2] = *da * dx[i__ + 2];
      dx[i__ + 3] = *da * dx[i__ + 3];
      dx[i__ + 4] = *da * dx[i__ + 4];
      /* L50: */
   }
   return 0;
} /* dscal_ */

integer idamax_(integer *n, real_type *dx, integer *incx)
{
   /* System generated locals */
   integer ret_val, i__1;
   real_type d__1;

   /* Local variables */
   static integer i__, ix;
   static real_type dmax__;


   /*     finds the index of element having max. dabsolute value. */
   /*     jack dongarra, linpack, 3/11/78. */


   /* Parameter adjustments */
   --dx;

   /* Function Body */
   ret_val = 0;
   if (*n < 1) {
      return ret_val;
   }
   ret_val = 1;
   if (*n == 1) {
      return ret_val;
   }
   if (*incx == 1) {
      goto L20;
   }

   /*        code for increment not equal to 1 */

   ix = 1;
   dmax__ = abs(dx[1]);
   ix += *incx;
   i__1 = *n;
   for (i__ = 2; i__ <= i__1; ++i__) {
      if ((d__1 = dx[ix], abs(d__1)) <= dmax__) {
         goto L5;
      }
      ret_val = i__;
      dmax__ = (d__1 = dx[ix], abs(d__1));
L5:
      ix += *incx;
      /* L10: */
   }
   return ret_val;

   /*        code for increment equal to 1 */

L20:
   dmax__ = abs(dx[1]);
   i__1 = *n;
   for (i__ = 2; i__ <= i__1; ++i__) {
      if ((d__1 = dx[i__], abs(d__1)) <= dmax__) {
         goto L30;
      }
      ret_val = i__;
      dmax__ = (d__1 = dx[i__], abs(d__1));
L30:
      ;
   }
   return ret_val;
} /* idamax_ */

real_type epslon_(real_type *x)
{
#if defined(TEST_MPF_100) || defined(TEST_MPFR_100) || defined(TEST_GMPXX) || defined(TEST_MPFRXX)
   return std::ldexp(1.0, 1 - ((100 + 1) * 1000L) / 301L);
#elif defined(TEST_CPP_DEC_FLOAT_BN)
   return std::pow(10.0, 1-std::numeric_limits<efx::cpp_dec_float_50>::digits10);
#elif defined(NATIVE_FLOAT128)
   return FLT128_EPSILON;
#else
   return CAST_TO_RT(std::numeric_limits<real_type>::epsilon());
#endif
} /* epslon_ */

/* Subroutine */ int mm_(real_type *a, integer *lda, integer *n1, integer *
   n3, real_type *b, integer *ldb, integer *n2, real_type *c__, 
   integer *ldc)
{
   /* System generated locals */
   integer a_dim1, a_offset, b_dim1, b_offset, c_dim1, c_offset, i__1, i__2;

   /* Local variables */
   static integer i__, j;


   /*   purpose: */
   /*     multiply matrix b times matrix c and store the result in matrix a. */

   /*   parameters: */

   /*     a double precision(lda,n3), matrix of n1 rows and n3 columns */

   /*     lda integer, leading dimension of array a */

   /*     n1 integer, number of rows in matrices a and b */

   /*     n3 integer, number of columns in matrices a and c */

   /*     b double precision(ldb,n2), matrix of n1 rows and n2 columns */

   /*     ldb integer, leading dimension of array b */

   /*     n2 integer, number of columns in matrix b, and number of rows in */
   /*         matrix c */

   /*     c double precision(ldc,n3), matrix of n2 rows and n3 columns */

   /*     ldc integer, leading dimension of array c */

   /* ---------------------------------------------------------------------- */

   /* Parameter adjustments */
   a_dim1 = *lda;
   a_offset = 1 + a_dim1;
   a -= a_offset;
   b_dim1 = *ldb;
   b_offset = 1 + b_dim1;
   b -= b_offset;
   c_dim1 = *ldc;
   c_offset = 1 + c_dim1;
   c__ -= c_offset;

   /* Function Body */
   i__1 = *n3;
   for (j = 1; j <= i__1; ++j) {
      i__2 = *n1;
      for (i__ = 1; i__ <= i__2; ++i__) {
         a[i__ + j * a_dim1] = CAST_TO_RT(0);
         /* L10: */
      }
      dmxpy_(n2, &a[j * a_dim1 + 1], n1, ldb, &c__[j * c_dim1 + 1], &b[
         b_offset]);
         /* L20: */
   }

   return 0;
} /* mm_ */

/* Subroutine */ int dmxpy_(integer *n1, real_type *y, integer *n2, integer *
   ldm, real_type *x, real_type *m)
{
   /* System generated locals */
   integer m_dim1, m_offset, i__1, i__2;

   /* Local variables */
   static integer i__, j, jmin;


   /*   purpose: */
   /*     multiply matrix m times vector x and add the result to vector y. */

   /*   parameters: */

   /*     n1 integer, number of elements in vector y, and number of rows in */
   /*         matrix m */

   /*     y double precision(n1), vector of length n1 to which is added */
   /*         the product m*x */

   /*     n2 integer, number of elements in vector x, and number of columns */
   /*         in matrix m */

   /*     ldm integer, leading dimension of array m */

   /*     x double precision(n2), vector of length n2 */

   /*     m double precision(ldm,n2), matrix of n1 rows and n2 columns */

   /* ---------------------------------------------------------------------- */

   /*   cleanup odd vector */

   /* Parameter adjustments */
   --y;
   m_dim1 = *ldm;
   m_offset = 1 + m_dim1;
   m -= m_offset;
   --x;

   /* Function Body */
   j = *n2 % 2;
   if (j >= 1) {
      i__1 = *n1;
      for (i__ = 1; i__ <= i__1; ++i__) {
         y[i__] += x[j] * m[i__ + j * m_dim1];
         /* L10: */
      }
   }

   /*   cleanup odd group of two vectors */

   j = *n2 % 4;
   if (j >= 2) {
      i__1 = *n1;
      for (i__ = 1; i__ <= i__1; ++i__) {
         y[i__] = y[i__] + x[j - 1] * m[i__ + (j - 1) * m_dim1] + x[j] * m[
            i__ + j * m_dim1];
            /* L20: */
      }
   }

   /*   cleanup odd group of four vectors */

   j = *n2 % 8;
   if (j >= 4) {
      i__1 = *n1;
      for (i__ = 1; i__ <= i__1; ++i__) {
         y[i__] = y[i__] + x[j - 3] * m[i__ + (j - 3) * m_dim1] + x[j - 2] 
         * m[i__ + (j - 2) * m_dim1] + x[j - 1] * m[i__ + (j - 1) *
            m_dim1] + x[j] * m[i__ + j * m_dim1];
         /* L30: */
      }
   }

   /*   cleanup odd group of eight vectors */

   j = *n2 % 16;
   if (j >= 8) {
      i__1 = *n1;
      for (i__ = 1; i__ <= i__1; ++i__) {
         y[i__] = y[i__] + x[j - 7] * m[i__ + (j - 7) * m_dim1] + x[j - 6] 
         * m[i__ + (j - 6) * m_dim1] + x[j - 5] * m[i__ + (j - 5) *
            m_dim1] + x[j - 4] * m[i__ + (j - 4) * m_dim1] + x[j - 3]
         * m[i__ + (j - 3) * m_dim1] + x[j - 2] * m[i__ + (j - 2) 
            * m_dim1] + x[j - 1] * m[i__ + (j - 1) * m_dim1] + x[j] * 
            m[i__ + j * m_dim1];
         /* L40: */
      }
   }

   /*   main loop - groups of sixteen vectors */

   jmin = j + 16;
   i__1 = *n2;
   for (j = jmin; j <= i__1; j += 16) {
      i__2 = *n1;
      for (i__ = 1; i__ <= i__2; ++i__) {
         y[i__] = y[i__] + x[j - 15] * m[i__ + (j - 15) * m_dim1] + x[j - 
            14] * m[i__ + (j - 14) * m_dim1] + x[j - 13] * m[i__ + (j 
            - 13) * m_dim1] + x[j - 12] * m[i__ + (j - 12) * m_dim1] 
         + x[j - 11] * m[i__ + (j - 11) * m_dim1] + x[j - 10] * m[
            i__ + (j - 10) * m_dim1] + x[j - 9] * m[i__ + (j - 9) * 
               m_dim1] + x[j - 8] * m[i__ + (j - 8) * m_dim1] + x[j - 7] 
            * m[i__ + (j - 7) * m_dim1] + x[j - 6] * m[i__ + (j - 6) *
               m_dim1] + x[j - 5] * m[i__ + (j - 5) * m_dim1] + x[j - 4]
            * m[i__ + (j - 4) * m_dim1] + x[j - 3] * m[i__ + (j - 3) 
               * m_dim1] + x[j - 2] * m[i__ + (j - 2) * m_dim1] + x[j - 
               1] * m[i__ + (j - 1) * m_dim1] + x[j] * m[i__ + j * 
               m_dim1];
            /* L50: */
      }
      /* L60: */
   }
   return 0;
} /* dmxpy_ */

real_type ran_(integer *iseed)
{
   /* System generated locals */
   real_type ret_val;

   /* Local variables */
   static integer it1, it2, it3, it4;


   /*     modified from the LAPACK auxiliary routine 10/12/92 JD */
   /*  -- LAPACK auxiliary routine (version 1.0) -- */
   /*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
   /*     Courant Institute, Argonne National Lab, and Rice University */
   /*     February 29, 1992 */

   /*     .. Array Arguments .. */
   /*     .. */

   /*  Purpose */
   /*  ======= */

   /*  DLARAN returns a random double number from a uniform (0,1) */
   /*  distribution. */

   /*  Arguments */
   /*  ========= */

   /*  ISEED   (input/output) INTEGER array, dimension (4) */
   /*          On entry, the seed of the random number generator; the array */
   /*          elements must be between 0 and 4095, and ISEED(4) must be */
   /*          odd. */
   /*          On exit, the seed is updated. */

   /*  Further Details */
   /*  =============== */

   /*  This routine uses a multiplicative congruential method with modulus */
   /*  2**48 and multiplier 33952834046453 (see G.S.Fishman, */
   /*  'Multiplicative congruential random number generators with modulus */
   /*  2**b: an exhaustive analysis for b = 32 and a partial analysis for */
   /*  b = 48', Math. Comp. 189, pp 331-344, 1990). */

   /*  48-bit integers are stored in 4 integer array elements with 12 bits */
   /*  per element. Hence the routine is portable across machines with */
   /*  integers of 32 bits or more. */

   /*     .. Parameters .. */
   /*     .. */
   /*     .. Local Scalars .. */
   /*     .. */
   /*     .. Intrinsic Functions .. */
   /*     .. */
   /*     .. Executable Statements .. */

   /*     multiply the seed by the multiplier modulo 2**48 */

   /* Parameter adjustments */
   --iseed;

   /* Function Body */
   it4 = iseed[4] * 2549;
   it3 = it4 / 4096;
   it4 -= it3 << 12;
   it3 = it3 + iseed[3] * 2549 + iseed[4] * 2508;
   it2 = it3 / 4096;
   it3 -= it2 << 12;
   it2 = it2 + iseed[2] * 2549 + iseed[3] * 2508 + iseed[4] * 322;
   it1 = it2 / 4096;
   it2 -= it1 << 12;
   it1 = it1 + iseed[1] * 2549 + iseed[2] * 2508 + iseed[3] * 322 + iseed[4] 
   * 494;
   it1 %= 4096;

   /*     return updated seed */

   iseed[1] = it1;
   iseed[2] = it2;
   iseed[3] = it3;
   iseed[4] = it4;

   /*     convert 48-bit integer to a double number in the interval (0,1) */

   ret_val = ((real_type) it1 + ((real_type) it2 + ((real_type) it3 + (
      real_type) it4 * 2.44140625e-4) * 2.44140625e-4) * 2.44140625e-4)
      * 2.44140625e-4;
   return ret_val;

   /*     End of RAN */

} /* ran_ */

/*

Double results:
~~~~~~~~~~~~~~

norm. resid      resid           machep         x(1)          x(n)
6.4915           7.207e-013      2.2204e-016    1             1



times are reported for matrices of order  1000
factor     solve      total     mflops       unit      ratio
times for array with leading dimension of1001
1.443     0.003      1.446     462.43       0.004325  25.821


mpf_class results:
~~~~~~~~~~~~~~~~~~

norm. resid      resid           machep         x(1)          x(n)
3.6575e-05       5.2257e-103     2.8575e-101    1             1



times are reported for matrices of order  1000
factor     solve      total     mflops       unit      ratio
times for array with leading dimension of1001
266.45     0.798      267.24    2.5021       0.79933   4772.2


number<gmp_float<100> >:
~~~~~~~~~~~~~~~~~~~~~~~~~~~

     norm. resid      resid           machep         x(1)          x(n)
  0.36575e-4          0.52257e-102   0.28575e-100    0.1e1         0.1e1



    times are reported for matrices of order  1000
      factor     solve      total     mflops       unit      ratio
 times for array with leading dimension of1001
      279.96        0.84       280.8      2.3813     0.83988      5014.3

boost::multiprecision::ef::cpp_dec_float_50:
~~~~~~~~~~~~~~~~~~~~~~~~~

     norm. resid      resid           machep         x(1)          x(n)
     2.551330735e-16  1.275665107e-112 1e-99         1             1



    times are reported for matrices of order  1000
      factor     solve      total     mflops       unit      ratio
 times for array with leading dimension of1001
      363.89      1.074     364.97    1.8321       1.0916    6517.3
*/