Add a packetmath unit test, re-enable the comma-initializer unit test,

and bug fix in PacketMath/SSE

6 files changed, 184 insertions, 7 deletions

diff --git a/Eigen/src/Core/arch/SSE/PacketMath.h b/Eigen/src/Core/arch/SSE/PacketMath.h
index 07aa84d69..f2744e340 100644
--- a/Eigen/src/Core/arch/SSE/PacketMath.h
+++ b/Eigen/src/Core/arch/SSE/PacketMath.h
@@ -219,7 +219,8 @@ struct ei_palign_impl<Offset,__m128>
 {
   inline static void run(__m128& first, const __m128& second)
   {
-    first = _mm_castsi128_ps(_mm_alignr_epi8(_mm_castps_si128(first), _mm_castps_si128(second), (4-Offset)*4));
+    if (Offset!=0)
+      first = _mm_castsi128_ps(_mm_alignr_epi8(_mm_castps_si128(second), _mm_castps_si128(first), (Offset)*4));
   }
 };
 
@@ -228,7 +229,18 @@ struct ei_palign_impl<Offset,__m128i>
 {
   inline static void run(__m128i& first, const __m128i& second)
   {
-    first = _mm_alignr_epi8(first, second, (4-Offset)*4);
+    if (Offset!=0)
+      first = _mm_alignr_epi8(second,first, (Offset)*4);
+  }
+};
+
+template<int Offset>
+struct ei_palign_impl<Offset,__m128d>
+{
+  inline static void run(__m128d& first, const __m128d& second)
+  {
+    if (Offset==1)
+      first = _mm_castsi128_pd(_mm_alignr_epi8(_mm_castpd_si128(second), _mm_castpd_si128(first), 8));
   }
 };
 #else
@@ -278,6 +290,19 @@ struct ei_palign_impl<Offset,__m128i>
     }
   }
 };
+
+template<int Offset>
+struct ei_palign_impl<Offset,__m128d>
+{
+  inline static void run(__m128d& first, const __m128d& second)
+  {
+    if (Offset==1)
+    {
+      first = _mm_castps_pd(_mm_movehl_ps(_mm_castpd_ps(first),_mm_castpd_ps(first)));
+      first = _mm_castps_pd(_mm_movelh_ps(_mm_castpd_ps(first),_mm_castpd_ps(second)));
+    }
+  }
+};
 #endif
 
 #endif // EIGEN_PACKET_MATH_SSE_H
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 30fe3e755..917715eb0 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -86,6 +86,7 @@ ENDIF(TEST_LIB)
 
 EI_ADD_TEST(sizeof)
 EI_ADD_TEST(nomalloc)
+EI_ADD_TEST(packetmath)
 EI_ADD_TEST(basicstuff)
 EI_ADD_TEST(linearstructure)
 EI_ADD_TEST(cwiseop)
@@ -95,6 +96,7 @@ EI_ADD_TEST(product_large ${EI_OFLAG})
 EI_ADD_TEST(adjoint)
 EI_ADD_TEST(submatrices)
 EI_ADD_TEST(miscmatrices)
+EI_ADD_TEST(commainitializer)
 EI_ADD_TEST(smallvectors)
 EI_ADD_TEST(map)
 EI_ADD_TEST(array)
diff --git a/test/cholesky.cpp b/test/cholesky.cpp
index 55a74c477..4bf28ef68 100644
--- a/test/cholesky.cpp
+++ b/test/cholesky.cpp
@@ -60,5 +60,6 @@ void test_cholesky()
     CALL_SUBTEST( cholesky(Matrix3f()) );
     CALL_SUBTEST( cholesky(Matrix4d()) );
     CALL_SUBTEST( cholesky(MatrixXcd(7,7)) );
+    CALL_SUBTEST( cholesky(MatrixXf(85,85)) );
   }
 }
diff --git a/test/commainitializer.cpp b/test/commainitializer.cpp
index fa3e3e348..257ecc21f 100644
--- a/test/commainitializer.cpp
+++ b/test/commainitializer.cpp
@@ -28,16 +28,16 @@ void test_commainitializer()
 {
   Matrix3d m3;
   Matrix4d m4;
-  VERIFY_RAISES_ASSERT(m4 = m3);
 
   VERIFY_RAISES_ASSERT( (m3 << 1, 2, 3, 4, 5, 6, 7, 8) );
   VERIFY_RAISES_ASSERT( (m3 << 1, 2, 3, 4, 5, 6, 7, 8, 9, 10) );
 
   double data[] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
+  Matrix3d ref = Map<Matrix<double,3,3,3,3,RowMajorBit> >(data);
 
   m3 = Matrix3d::Random();
   m3 << 1, 2, 3, 4, 5, 6, 7, 8, 9;
-  VERIFY_IS_APPROX(m3, (Matrix<double,3,3,RowMajorBit>::map(data)) );
+  VERIFY_IS_APPROX(m3, ref );
 
   Vector3d vec[3];
   vec[0] << 1, 4, 7;
@@ -45,7 +45,7 @@ void test_commainitializer()
   vec[2] << 3, 6, 9;
   m3 = Matrix3d::Random();
   m3 << vec[0], vec[1], vec[2];
-  VERIFY_IS_APPROX(m3, (Matrix<double,3,3,RowMajorBit>::map(data)) );
+  VERIFY_IS_APPROX(m3, ref);
 
   vec[0] << 1, 2, 3;
   vec[1] << 4, 5, 6;
@@ -54,5 +54,5 @@ void test_commainitializer()
   m3 << vec[0].transpose(),
         4, 5, 6,
         vec[2].transpose();
-  VERIFY_IS_APPROX(m3, (Matrix<double,3,3,RowMajorBit>::map(data)) );
+  VERIFY_IS_APPROX(m3, ref);
 }
diff --git a/test/packetmath.cpp b/test/packetmath.cpp
new file mode 100644
index 000000000..c282762d6
--- /dev/null
+++ b/test/packetmath.cpp
@@ -0,0 +1,149 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra. Eigen itself is part of the KDE project.
+//
+// Copyright (C) 2006-2008 Benoit Jacob <jacob@math.jussieu.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#include "main.h"
+
+// using namespace Eigen;
+
+template<typename Scalar> bool areApprox(const Scalar* a, const Scalar* b, int size)
+{
+  for (int i=0; i<size; ++i)
+    if (!ei_isApprox(a[i],b[i])) return false;
+  return true;
+}
+
+#define CHECK_CWISE(REFOP, POP) { \
+  for (int i=0; i<PacketSize; ++i) \
+    ref[i] = REFOP(data1[i], data1[i+PacketSize]); \
+  ei_pstore(data2, POP(ei_pload(data1), ei_pload(data1+PacketSize))); \
+  VERIFY(areApprox(ref, data2, PacketSize) && #POP); \
+}
+
+#define REF_ADD(a,b) ((a)+(b))
+#define REF_SUB(a,b) ((a)-(b))
+#define REF_MUL(a,b) ((a)*(b))
+#define REF_DIV(a,b) ((a)/(b))
+
+namespace std {
+
+template<> const complex<float>& min(const complex<float>& a, const complex<float>& b)
+{ return a.real() < b.real() ? a : b; }
+
+template<> const complex<float>& max(const complex<float>& a, const complex<float>& b)
+{ return a.real() < b.real() ? b : a; }
+
+}
+
+template<typename Scalar> void packetmath()
+{
+  typedef typename ei_packet_traits<Scalar>::type Packet;
+  const int PacketSize = ei_packet_traits<Scalar>::size;
+
+  const int size = PacketSize*4;
+  Scalar data1[ei_packet_traits<Scalar>::size*4];
+  Scalar data2[ei_packet_traits<Scalar>::size*4];
+  Packet packets[PacketSize];
+  Scalar ref[ei_packet_traits<Scalar>::size*4];
+  for (int i=0; i<size; ++i)
+  {
+    data1[i] = ei_random<Scalar>();
+    data2[i] = ei_random<Scalar>();
+  }
+
+  ei_pstore(data2, ei_pload(data1));
+  VERIFY(areApprox(data1, data2, PacketSize) && "aligned load/store");
+
+  for (int offset=0; offset<PacketSize; ++offset)
+  {
+    ei_pstore(data2, ei_ploadu(data1+offset));
+    VERIFY(areApprox(data1+offset, data2, PacketSize) && "ei_ploadu");
+  }
+
+  for (int offset=0; offset<PacketSize; ++offset)
+  {
+    ei_pstoreu(data2+offset, ei_pload(data1));
+    VERIFY(areApprox(data1, data2+offset, PacketSize) && "ei_pstoreu");
+  }
+
+  if (!ei_is_same_type<Scalar,double>::ret)
+  {
+    for (int offset=0; offset<PacketSize; ++offset)
+    {
+      packets[0] = ei_pload(data1);
+      packets[1] = ei_pload(data1+PacketSize);
+          if (offset==0) ei_palign<0>(packets[0], packets[1]);
+      else if (offset==1) ei_palign<1>(packets[0], packets[1]);
+      else if (offset==2) ei_palign<2>(packets[0], packets[1]);
+      else if (offset==3) ei_palign<3>(packets[0], packets[1]);
+      ei_pstore(data2, packets[0]);
+
+      for (int i=0; i<PacketSize; ++i)
+        ref[i] = data1[i+offset];
+
+      VERIFY(areApprox(ref, data2, PacketSize) && "ei_palign");
+    }
+  }
+
+  CHECK_CWISE(REF_ADD,  ei_padd);
+  CHECK_CWISE(REF_SUB,  ei_psub);
+  CHECK_CWISE(REF_MUL,  ei_pmul);
+  if (!ei_is_same_type<Scalar,int>::ret)
+    CHECK_CWISE(REF_DIV,  ei_pdiv);
+  CHECK_CWISE(std::min, ei_pmin);
+  CHECK_CWISE(std::max, ei_pmax);
+
+  for (int i=0; i<PacketSize; ++i)
+    ref[i] = data1[0];
+  ei_pstore(data2, ei_pset1(data1[0]));
+  VERIFY(areApprox(ref, data2, PacketSize) && "ei_pset1");
+
+  VERIFY(ei_isApprox(data1[0], ei_pfirst(ei_pload(data1))) && "ei_pfirst");
+
+  ref[0] = 0;
+  for (int i=0; i<PacketSize; ++i)
+    ref[0] += data1[i];
+  VERIFY(ei_isApprox(ref[0], ei_predux(ei_pload(data1))) && "ei_predux");
+
+  for (int j=0; j<PacketSize; ++j)
+  {
+    ref[j] = 0;
+    for (int i=0; i<PacketSize; ++i)
+      ref[j] += data1[i+j*PacketSize];
+    packets[j] = ei_pload(data1+j*PacketSize);
+  }
+  ei_pstore(data2, ei_preduxp(packets));
+  VERIFY(areApprox(ref, data2, PacketSize) && "ei_preduxp");
+
+  
+}
+
+void test_packetmath()
+{
+  for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST( packetmath<float>() );
+    CALL_SUBTEST( packetmath<double>() );
+    CALL_SUBTEST( packetmath<int>() );
+    packetmath<std::complex<float> >();
+  }
+}
diff --git a/test/triangular.cpp b/test/triangular.cpp
index de3b85537..fd744114c 100644
--- a/test/triangular.cpp
+++ b/test/triangular.cpp
@@ -101,6 +101,6 @@ void test_triangular()
     CALL_SUBTEST( triangular(Matrix3d()) );
     CALL_SUBTEST( triangular(MatrixXcf(4, 4)) );
     CALL_SUBTEST( triangular(Matrix<std::complex<float>,8, 8>()) );
-    CALL_SUBTEST( triangular(MatrixXf(12,12)) );
+    CALL_SUBTEST( triangular(MatrixXf(85,85)) );
   }
 }