Merged latest updates from the Eigen trunk.

8 files changed, 209 insertions, 24 deletions

diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp
index 06a6a8654..3851f9df2 100644
--- a/test/eigensolver_selfadjoint.cpp
+++ b/test/eigensolver_selfadjoint.cpp
@@ -29,7 +29,21 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
   MatrixType a = MatrixType::Random(rows,cols);
   MatrixType a1 = MatrixType::Random(rows,cols);
   MatrixType symmA =  a.adjoint() * a + a1.adjoint() * a1;
+  MatrixType symmC = symmA;
+  
+  // randomly nullify some rows/columns
+  {
+    Index count = 1;//internal::random<Index>(-cols,cols);
+    for(Index k=0; k<count; ++k)
+    {
+      Index i = internal::random<Index>(0,cols-1);
+      symmA.row(i).setZero();
+      symmA.col(i).setZero();
+    }
+  }
+  
   symmA.template triangularView<StrictlyUpper>().setZero();
+  symmC.template triangularView<StrictlyUpper>().setZero();
 
   MatrixType b = MatrixType::Random(rows,cols);
   MatrixType b1 = MatrixType::Random(rows,cols);
@@ -40,7 +54,7 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
   SelfAdjointEigenSolver<MatrixType> eiDirect;
   eiDirect.computeDirect(symmA);
   // generalized eigen pb
-  GeneralizedSelfAdjointEigenSolver<MatrixType> eiSymmGen(symmA, symmB);
+  GeneralizedSelfAdjointEigenSolver<MatrixType> eiSymmGen(symmC, symmB);
 
   VERIFY_IS_EQUAL(eiSymm.info(), Success);
   VERIFY((symmA.template selfadjointView<Lower>() * eiSymm.eigenvectors()).isApprox(
@@ -57,27 +71,28 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
   VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmNoEivecs.eigenvalues());
   
   // generalized eigen problem Ax = lBx
-  eiSymmGen.compute(symmA, symmB,Ax_lBx);
+  eiSymmGen.compute(symmC, symmB,Ax_lBx);
   VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
-  VERIFY((symmA.template selfadjointView<Lower>() * eiSymmGen.eigenvectors()).isApprox(
+  VERIFY((symmC.template selfadjointView<Lower>() * eiSymmGen.eigenvectors()).isApprox(
           symmB.template selfadjointView<Lower>() * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
 
   // generalized eigen problem BAx = lx
-  eiSymmGen.compute(symmA, symmB,BAx_lx);
+  eiSymmGen.compute(symmC, symmB,BAx_lx);
   VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
-  VERIFY((symmB.template selfadjointView<Lower>() * (symmA.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox(
+  VERIFY((symmB.template selfadjointView<Lower>() * (symmC.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox(
          (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
 
   // generalized eigen problem ABx = lx
-  eiSymmGen.compute(symmA, symmB,ABx_lx);
+  eiSymmGen.compute(symmC, symmB,ABx_lx);
   VERIFY_IS_EQUAL(eiSymmGen.info(), Success);
-  VERIFY((symmA.template selfadjointView<Lower>() * (symmB.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox(
+  VERIFY((symmC.template selfadjointView<Lower>() * (symmB.template selfadjointView<Lower>() * eiSymmGen.eigenvectors())).isApprox(
          (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
 
 
+  eiSymm.compute(symmC);
   MatrixType sqrtSymmA = eiSymm.operatorSqrt();
-  VERIFY_IS_APPROX(MatrixType(symmA.template selfadjointView<Lower>()), sqrtSymmA*sqrtSymmA);
-  VERIFY_IS_APPROX(sqrtSymmA, symmA.template selfadjointView<Lower>()*eiSymm.operatorInverseSqrt());
+  VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), sqrtSymmA*sqrtSymmA);
+  VERIFY_IS_APPROX(sqrtSymmA, symmC.template selfadjointView<Lower>()*eiSymm.operatorInverseSqrt());
 
   MatrixType id = MatrixType::Identity(rows, cols);
   VERIFY_IS_APPROX(id.template selfadjointView<Lower>().operatorNorm(), RealScalar(1));
@@ -95,9 +110,9 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
   VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt());
 
   // test Tridiagonalization's methods
-  Tridiagonalization<MatrixType> tridiag(symmA);
+  Tridiagonalization<MatrixType> tridiag(symmC);
   // FIXME tridiag.matrixQ().adjoint() does not work
-  VERIFY_IS_APPROX(MatrixType(symmA.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT().eval() * MatrixType(tridiag.matrixQ()).adjoint());
+  VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT().eval() * MatrixType(tridiag.matrixQ()).adjoint());
   
   // Test computation of eigenvalues from tridiagonal matrix
   if(rows > 1)
@@ -111,8 +126,8 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
   if (rows > 1)
   {
     // Test matrix with NaN
-    symmA(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
-    SelfAdjointEigenSolver<MatrixType> eiSymmNaN(symmA);
+    symmC(0,0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();
+    SelfAdjointEigenSolver<MatrixType> eiSymmNaN(symmC);
     VERIFY_IS_EQUAL(eiSymmNaN.info(), NoConvergence);
   }
 }
@@ -122,8 +137,10 @@ void test_eigensolver_selfadjoint()
   int s = 0;
   for(int i = 0; i < g_repeat; i++) {
     // very important to test 3x3 and 2x2 matrices since we provide special paths for them
+    CALL_SUBTEST_1( selfadjointeigensolver(Matrix2f()) );
     CALL_SUBTEST_1( selfadjointeigensolver(Matrix2d()) );
     CALL_SUBTEST_1( selfadjointeigensolver(Matrix3f()) );
+    CALL_SUBTEST_1( selfadjointeigensolver(Matrix3d()) );
     CALL_SUBTEST_2( selfadjointeigensolver(Matrix4d()) );
     s = internal::random<int>(1,EIGEN_TEST_MAX_SIZE/4);
     CALL_SUBTEST_3( selfadjointeigensolver(MatrixXf(s,s)) );
diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp
index 36721b496..cd04db5be 100644
--- a/test/jacobisvd.cpp
+++ b/test/jacobisvd.cpp
@@ -315,16 +315,30 @@ void jacobisvd_inf_nan()
   VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf));
   svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV);
 
-  Scalar some_nan = zero<Scalar>() / zero<Scalar>();
-  VERIFY(some_nan != some_nan);
-  svd.compute(MatrixType::Constant(10,10,some_nan), ComputeFullU | ComputeFullV);
+  Scalar nan = std::numeric_limits<Scalar>::quiet_NaN();
+  VERIFY(nan != nan);
+  svd.compute(MatrixType::Constant(10,10,nan), ComputeFullU | ComputeFullV);
 
   MatrixType m = MatrixType::Zero(10,10);
   m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_inf;
   svd.compute(m, ComputeFullU | ComputeFullV);
 
   m = MatrixType::Zero(10,10);
-  m(internal::random<int>(0,9), internal::random<int>(0,9)) = some_nan;
+  m(internal::random<int>(0,9), internal::random<int>(0,9)) = nan;
+  svd.compute(m, ComputeFullU | ComputeFullV);
+  
+  // regression test for bug 791
+  m.resize(3,3);
+  m << 0,    2*NumTraits<Scalar>::epsilon(),  0.5,
+       0,   -0.5,                             0,
+       nan,  0,                               0;
+  svd.compute(m, ComputeFullU | ComputeFullV);
+  
+  m.resize(4,4);
+  m <<  1, 0, 0, 0,
+        0, 3, 1, 2e-308,
+        1, 0, 1, nan,
+        0, nan, nan, 0;
   svd.compute(m, ComputeFullU | ComputeFullV);
 }
 
@@ -340,11 +354,33 @@ void jacobisvd_underoverflow()
   Matrix2d M;
   M << -7.90884e-313, -4.94e-324,
                  0, 5.60844e-313;
+  JacobiSVD<Matrix2d> svd;
+  svd.compute(M,ComputeFullU|ComputeFullV);
+  jacobisvd_check_full(M,svd);
+  
+  VectorXd value_set(9);
+  value_set << 0, 1, -1, 5.60844e-313, -5.60844e-313, 4.94e-324, -4.94e-324, -4.94e-223, 4.94e-223;
+  Array4i id(0,0,0,0);
+  int k = 0;
+  do
+  {
+    M << value_set(id(0)), value_set(id(1)), value_set(id(2)), value_set(id(3));
+    svd.compute(M,ComputeFullU|ComputeFullV);
+    jacobisvd_check_full(M,svd);
+    
+    id(k)++;
+    if(id(k)>=value_set.size())
+    {
+      while(k<3 && id(k)>=value_set.size()) id(++k)++;
+      id.head(k).setZero();
+      k=0;
+    }
+    
+  } while((id<int(value_set.size())).all());
+  
 #if defined __INTEL_COMPILER
 #pragma warning pop
 #endif
-  JacobiSVD<Matrix2d> svd;
-  svd.compute(M); // just check we don't loop indefinitely
   
   // Check for overflow:
   Matrix3d M3;
@@ -353,7 +389,8 @@ void jacobisvd_underoverflow()
        -8.7190887618028355e+307, -7.3453213709232193e+307, -2.4367363684472105e+307;
 
   JacobiSVD<Matrix3d> svd3;
-  svd3.compute(M3); // just check we don't loop indefinitely
+  svd3.compute(M3,ComputeFullU|ComputeFullV); // just check we don't loop indefinitely
+  jacobisvd_check_full(M3,svd3);
 }
 
 void jacobisvd_preallocate()
@@ -437,6 +474,7 @@ void test_jacobisvd()
 
     // Test on inf/nan matrix
     CALL_SUBTEST_7( jacobisvd_inf_nan<MatrixXf>() );
+    CALL_SUBTEST_10( jacobisvd_inf_nan<MatrixXd>() );
   }
 
   CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) ));
diff --git a/test/linearstructure.cpp b/test/linearstructure.cpp
index 618984d5c..b627915ce 100644
--- a/test/linearstructure.cpp
+++ b/test/linearstructure.cpp
@@ -2,11 +2,16 @@
 // for linear algebra.
 //
 // Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
+static bool g_called;
+
+#define EIGEN_SPECIAL_SCALAR_MULTIPLE_PLUGIN { g_called = true; }
+
 #include "main.h"
 
 template<typename MatrixType> void linearStructure(const MatrixType& m)
@@ -68,6 +73,24 @@ template<typename MatrixType> void linearStructure(const MatrixType& m)
   VERIFY_IS_APPROX(m1.block(0,0,rows,cols) * s1, m1 * s1);
 }
 
+// Make sure that complex * real and real * complex are properly optimized
+template<typename MatrixType> void real_complex(DenseIndex rows = MatrixType::RowsAtCompileTime, DenseIndex cols = MatrixType::ColsAtCompileTime)
+{
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename MatrixType::RealScalar RealScalar;
+  
+  RealScalar s = internal::random<RealScalar>();
+  MatrixType m1 = MatrixType::Random(rows, cols);
+  
+  g_called = false;
+  VERIFY_IS_APPROX(s*m1, Scalar(s)*m1);
+  VERIFY(g_called && "real * matrix<complex> not properly optimized");
+  
+  g_called = false;
+  VERIFY_IS_APPROX(m1*s, m1*Scalar(s));
+  VERIFY(g_called && "matrix<complex> * real not properly optimized");
+}
+
 void test_linearstructure()
 {
   for(int i = 0; i < g_repeat; i++) {
@@ -80,5 +103,8 @@ void test_linearstructure()
     CALL_SUBTEST_7( linearStructure(MatrixXi (internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
     CALL_SUBTEST_8( linearStructure(MatrixXcd(internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2), internal::random<int>(1,EIGEN_TEST_MAX_SIZE/2))) );
     CALL_SUBTEST_9( linearStructure(ArrayXXf (internal::random<int>(1,EIGEN_TEST_MAX_SIZE), internal::random<int>(1,EIGEN_TEST_MAX_SIZE))) );
+    
+    CALL_SUBTEST_10( real_complex<Matrix4cd>() );
+    CALL_SUBTEST_10( real_complex<MatrixXcf>(10,10) );
   }
 }
diff --git a/test/main.h b/test/main.h
index 763cec8f9..b504970f3 100644
--- a/test/main.h
+++ b/test/main.h
@@ -17,13 +17,36 @@
 #include <sstream>
 #include <vector>
 #include <typeinfo>
+
+// The following includes of STL headers have to be done _before_ the
+// definition of macros min() and max().  The reason is that many STL
+// implementations will not work properly as the min and max symbols collide
+// with the STL functions std:min() and std::max().  The STL headers may check
+// for the macro definition of min/max and issue a warning or undefine the
+// macros.
+//
+// Still, Windows defines min() and max() in windef.h as part of the regular
+// Windows system interfaces and many other Windows APIs depend on these
+// macros being available.  To prevent the macro expansion of min/max and to
+// make Eigen compatible with the Windows environment all function calls of
+// std::min() and std::max() have to be written with parenthesis around the
+// function name.
+//
+// All STL headers used by Eigen should be included here.  Because main.h is
+// included before any Eigen header and because the STL headers are guarded
+// against multiple inclusions, no STL header will see our own min/max macro
+// definitions.
 #include <limits>
 #include <algorithm>
-#include <sstream>
 #include <complex>
 #include <deque>
 #include <queue>
+#include <list>
 
+// To test that all calls from Eigen code to std::min() and std::max() are
+// protected by parenthesis against macro expansion, the min()/max() macros
+// are defined here and any not-parenthesized min/max call will cause a
+// compiler error.
 #define min(A,B) please_protect_your_min_with_parentheses
 #define max(A,B) please_protect_your_max_with_parentheses
 
@@ -76,6 +99,10 @@ namespace Eigen
 
 #define EIGEN_DEFAULT_IO_FORMAT IOFormat(4, 0, "  ", "\n", "", "", "", "")
 
+#if (defined(_CPPUNWIND) || defined(__EXCEPTIONS)) && !defined(__CUDA_ARCH__)
+  #define EIGEN_EXCEPTIONS
+#endif
+
 #ifndef EIGEN_NO_ASSERTION_CHECKING
 
   namespace Eigen
@@ -172,7 +199,7 @@ namespace Eigen
 
 #ifndef VERIFY_RAISES_ASSERT
   #define VERIFY_RAISES_ASSERT(a) \
-    std::cout << "Can't VERIFY_RAISES_ASSERT( " #a " ) with exceptions disabled";
+    std::cout << "Can't VERIFY_RAISES_ASSERT( " #a " ) with exceptions disabled\n";
 #endif
     
   #if !defined(__CUDACC__)
diff --git a/test/product.h b/test/product.h
index 856b234ac..0b3abe402 100644
--- a/test/product.h
+++ b/test/product.h
@@ -139,4 +139,12 @@ template<typename MatrixType> void product(const MatrixType& m)
   // inner product
   Scalar x = square2.row(c) * square2.col(c2);
   VERIFY_IS_APPROX(x, square2.row(c).transpose().cwiseProduct(square2.col(c2)).sum());
+  
+  // outer product
+  VERIFY_IS_APPROX(m1.col(c) * m1.row(r), m1.block(0,c,rows,1) * m1.block(r,0,1,cols));
+  VERIFY_IS_APPROX(m1.row(r).transpose() * m1.col(c).transpose(), m1.block(r,0,1,cols).transpose() * m1.block(0,c,rows,1).transpose());
+  VERIFY_IS_APPROX(m1.block(0,c,rows,1) * m1.row(r), m1.block(0,c,rows,1) * m1.block(r,0,1,cols));
+  VERIFY_IS_APPROX(m1.col(c) * m1.block(r,0,1,cols), m1.block(0,c,rows,1) * m1.block(r,0,1,cols));
+  VERIFY_IS_APPROX(m1.leftCols(1) * m1.row(r), m1.block(0,0,rows,1) * m1.block(r,0,1,cols));
+  VERIFY_IS_APPROX(m1.col(c) * m1.topRows(1), m1.block(0,c,rows,1) * m1.block(0,0,1,cols));  
 }
diff --git a/test/sparseqr.cpp b/test/sparseqr.cpp
index 1fe4a98ee..8e6887dd3 100644
--- a/test/sparseqr.cpp
+++ b/test/sparseqr.cpp
@@ -54,6 +54,8 @@ template<typename Scalar> void test_sparseqr_scalar()
   
   b = dA * DenseVector::Random(A.cols());
   solver.compute(A);
+  if(internal::random<float>(0,1)>0.5)
+    solver.factorize(A);  // this checks that calling analyzePattern is not needed if the pattern do not change.
   if (solver.info() != Success)
   {
     std::cerr << "sparse QR factorization failed\n";
diff --git a/test/stable_norm.cpp b/test/stable_norm.cpp
index 549f91fbf..f76919af4 100644
--- a/test/stable_norm.cpp
+++ b/test/stable_norm.cpp
@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2009-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
@@ -14,6 +14,21 @@ template<typename T> bool isNotNaN(const T& x)
   return x==x;
 }
 
+template<typename T> bool isNaN(const T& x)
+{
+  return x!=x;
+}
+
+template<typename T> bool isInf(const T& x)
+{
+  return x > NumTraits<T>::highest();
+}
+
+template<typename T> bool isMinusInf(const T& x)
+{
+  return x < NumTraits<T>::lowest();
+}
+
 // workaround aggressive optimization in ICC
 template<typename T> EIGEN_DONT_INLINE  T sub(T a, T b) { return a - b; }
 
@@ -106,6 +121,58 @@ template<typename MatrixType> void stable_norm(const MatrixType& m)
   VERIFY_IS_APPROX(vrand.rowwise().stableNorm(),      vrand.rowwise().norm());
   VERIFY_IS_APPROX(vrand.rowwise().blueNorm(),        vrand.rowwise().norm());
   VERIFY_IS_APPROX(vrand.rowwise().hypotNorm(),       vrand.rowwise().norm());
+  
+  // test NaN, +inf, -inf 
+  MatrixType v;
+  Index i = internal::random<Index>(0,rows-1);
+  Index j = internal::random<Index>(0,cols-1);
+
+  // NaN
+  {
+    v = vrand;
+    v(i,j) = RealScalar(0)/RealScalar(0);
+    VERIFY(!isFinite(v.squaredNorm()));   VERIFY(isNaN(v.squaredNorm()));
+    VERIFY(!isFinite(v.norm()));          VERIFY(isNaN(v.norm()));
+    VERIFY(!isFinite(v.stableNorm()));    VERIFY(isNaN(v.stableNorm()));
+    VERIFY(!isFinite(v.blueNorm()));      VERIFY(isNaN(v.blueNorm()));
+    VERIFY(!isFinite(v.hypotNorm()));     VERIFY(isNaN(v.hypotNorm()));
+  }
+  
+  // +inf
+  {
+    v = vrand;
+    v(i,j) = RealScalar(1)/RealScalar(0);
+    VERIFY(!isFinite(v.squaredNorm()));   VERIFY(isInf(v.squaredNorm()));
+    VERIFY(!isFinite(v.norm()));          VERIFY(isInf(v.norm()));
+    VERIFY(!isFinite(v.stableNorm()));    VERIFY(isInf(v.stableNorm()));
+    VERIFY(!isFinite(v.blueNorm()));      VERIFY(isInf(v.blueNorm()));
+    VERIFY(!isFinite(v.hypotNorm()));     VERIFY(isInf(v.hypotNorm()));
+  }
+  
+  // -inf
+  {
+    v = vrand;
+    v(i,j) = RealScalar(-1)/RealScalar(0);
+    VERIFY(!isFinite(v.squaredNorm()));   VERIFY(isInf(v.squaredNorm()));
+    VERIFY(!isFinite(v.norm()));          VERIFY(isInf(v.norm()));
+    VERIFY(!isFinite(v.stableNorm()));    VERIFY(isInf(v.stableNorm()));
+    VERIFY(!isFinite(v.blueNorm()));      VERIFY(isInf(v.blueNorm()));
+    VERIFY(!isFinite(v.hypotNorm()));     VERIFY(isInf(v.hypotNorm()));
+  }
+  
+  // mix
+  {
+    Index i2 = internal::random<Index>(0,rows-1);
+    Index j2 = internal::random<Index>(0,cols-1);
+    v = vrand;
+    v(i,j) = RealScalar(-1)/RealScalar(0);
+    v(i2,j2) = RealScalar(0)/RealScalar(0);
+    VERIFY(!isFinite(v.squaredNorm()));   VERIFY(isNaN(v.squaredNorm()));
+    VERIFY(!isFinite(v.norm()));          VERIFY(isNaN(v.norm()));
+    VERIFY(!isFinite(v.stableNorm()));    VERIFY(isNaN(v.stableNorm()));
+    VERIFY(!isFinite(v.blueNorm()));      VERIFY(isNaN(v.blueNorm()));
+    VERIFY(!isFinite(v.hypotNorm()));     VERIFY(isNaN(v.hypotNorm()));
+  }
 }
 
 void test_stable_norm()
diff --git a/test/upperbidiagonalization.cpp b/test/upperbidiagonalization.cpp
index d15bf588b..847b34b55 100644
--- a/test/upperbidiagonalization.cpp
+++ b/test/upperbidiagonalization.cpp
@@ -35,7 +35,7 @@ void test_upperbidiagonalization()
    CALL_SUBTEST_1( upperbidiag(MatrixXf(3,3)) );
    CALL_SUBTEST_2( upperbidiag(MatrixXd(17,12)) );
    CALL_SUBTEST_3( upperbidiag(MatrixXcf(20,20)) );
-   CALL_SUBTEST_4( upperbidiag(MatrixXcd(16,15)) );
+   CALL_SUBTEST_4( upperbidiag(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(16,15)) );
    CALL_SUBTEST_5( upperbidiag(Matrix<float,6,4>()) );
    CALL_SUBTEST_6( upperbidiag(Matrix<float,5,5>()) );
    CALL_SUBTEST_7( upperbidiag(Matrix<double,4,3>()) );