From 424c3ad26654c1cdaa526d73fda435ccfb0b5c79 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Thu, 17 Jul 2014 09:41:33 +0200
Subject: [PATCH] bug #842: fix specialized product for mpreal

---
 unsupported/Eigen/MPRealSupport     | 20 +++-----------------
 unsupported/test/mpreal_support.cpp |  4 ++--
 2 files changed, 5 insertions(+), 19 deletions(-)

diff --git a/unsupported/Eigen/MPRealSupport b/unsupported/Eigen/MPRealSupport
index d4b03647d..6a65a601b 100644
--- a/unsupported/Eigen/MPRealSupport
+++ b/unsupported/Eigen/MPRealSupport
@@ -139,9 +139,8 @@ int main()
     public:
       typedef mpfr::mpreal ResScalar;
       enum {
-        nr = 2, // must be 2 for proper packing...
+        nr = 1,
         mr = 1,
-        WorkSpaceFactor = nr,
         LhsProgress = 1,
         RhsProgress = 1
       };
@@ -154,9 +153,9 @@ int main()
 
       EIGEN_DONT_INLINE
       void operator()(mpreal* res, Index resStride, const mpreal* blockA, const mpreal* blockB, Index rows, Index depth, Index cols, mpreal alpha,
-                      Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0, mpreal* /*unpackedB*/ = 0)
+                      Index strideA=-1, Index strideB=-1, Index offsetA=0, Index offsetB=0)
       {
-        mpreal acc1, acc2, tmp;
+        mpreal acc1, tmp;
         
         if(strideA==-1) strideA = depth;
         if(strideB==-1) strideB = depth;
@@ -165,33 +164,20 @@ int main()
         {
           Index actual_nr = (std::min<Index>)(nr,cols-j);
           mpreal *C1 = res + j*resStride;
-          mpreal *C2 = res + (j+1)*resStride;
           for(Index i=0; i<rows; i++)
           {
             mpreal *B = const_cast<mpreal*>(blockB) + j*strideB + offsetB*actual_nr;
             mpreal *A = const_cast<mpreal*>(blockA) + i*strideA + offsetA;
             acc1 = 0;
-            acc2 = 0;
             for(Index k=0; k<depth; k++)
             {
               mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_ptr(), B[0].mpfr_ptr(), mpreal::get_default_rnd());
               mpfr_add(acc1.mpfr_ptr(), acc1.mpfr_ptr(), tmp.mpfr_ptr(),  mpreal::get_default_rnd());
-              
-              if(actual_nr==2) {
-                mpfr_mul(tmp.mpfr_ptr(), A[k].mpfr_ptr(), B[1].mpfr_ptr(), mpreal::get_default_rnd());
-                mpfr_add(acc2.mpfr_ptr(), acc2.mpfr_ptr(), tmp.mpfr_ptr(),  mpreal::get_default_rnd());
-              }
-              
               B+=actual_nr;
             }
             
             mpfr_mul(acc1.mpfr_ptr(), acc1.mpfr_ptr(), alpha.mpfr_ptr(), mpreal::get_default_rnd());
             mpfr_add(C1[i].mpfr_ptr(), C1[i].mpfr_ptr(), acc1.mpfr_ptr(),  mpreal::get_default_rnd());
-            
-            if(actual_nr==2) {
-              mpfr_mul(acc2.mpfr_ptr(), acc2.mpfr_ptr(), alpha.mpfr_ptr(), mpreal::get_default_rnd());
-              mpfr_add(C2[i].mpfr_ptr(), C2[i].mpfr_ptr(), acc2.mpfr_ptr(),  mpreal::get_default_rnd());
-            }
           }
         }
       }
diff --git a/unsupported/test/mpreal_support.cpp b/unsupported/test/mpreal_support.cpp
index bc00382be..13237886e 100644
--- a/unsupported/test/mpreal_support.cpp
+++ b/unsupported/test/mpreal_support.cpp
@@ -1,3 +1,4 @@
+#include <bench/BenchTimer.h>
 #include "main.h"
 #include <Eigen/MPRealSupport>
 #include <Eigen/LU>
@@ -32,12 +33,11 @@ void test_mpreal_support()
     VERIFY_IS_APPROX(A.array().abs2().sqrt(), A.array().abs());
     VERIFY_IS_APPROX(A.array().sin(),         sin(A.array()));
     VERIFY_IS_APPROX(A.array().cos(),         cos(A.array()));
-    
 
     // Cholesky
     X = S.selfadjointView<Lower>().llt().solve(B);
     VERIFY_IS_APPROX((S.selfadjointView<Lower>()*X).eval(),B);
-
+    
     // partial LU
     X = A.lu().solve(B);
     VERIFY_IS_APPROX((A*X).eval(),B);