Hi,

I'm working on porting a library of mine from LAPACK++ to Eigen. In the past, I've relied heavily on LAPACK++'s submatrix views (via the LaIndex class), which are represented by blocks in Eigen. I also created sub-blocks of blocks, which was trivial in LAPACK++ because all submatrix views also are objects of type LaGenMatDouble, the default matrix type. As VectorXd and Block<VectorXd> seem to be not interoperable (as stated in http://eigen.tuxfamily.org/dox/TopicFun ... Types.html), I figured I could just use Block<VectorXd> as arguments for some of my subroutines. However, there's still a problem with recursion: It seems like creating a block from a block generates an object of type Block< Block<VectorXd> > rather than a new Block<VectorXd>. Is there any way to avoid this in order to e.g. restrict blocks to an even smaller sub-range of the underlying vector? In this case, using argument templates wouldn't help because it would lead to an "infinite recursion" of templates, i.e. Block<VectorXd>, Block< Block<VectorXd> >, Block< Block<Block<VectorXd> > > etc.

A similar question: Some of my functions take a vector (or a Block<VectorXd>, for that matter) as argument, but I'd like to use them on single columns of matrices, too. So is there a straightforward way to convert a single column (or row) of a matrix to a Block<VectorXd>?

Best,

Daniel

Currently the only general way to accept any kind of expression is to 'templatize' your functions, but there is a plane to address your use case which appears to be quite common:
http://listengine.tuxfamily.org/lists.t ... 00062.html

Thank you for the tip. It appears to me that the following code

Code: Select all

typedef Stride< Dynamic, Dynamic > DynamicStride;

template <typename OtherDerived>
Map< Matrix< typename internal::traits<OtherDerived>::Scalar, Dynamic, Dynamic>, 0, DynamicStride > toMap(const OtherDerived &source)
{
   return Map< Matrix<typename internal::traits<OtherDerived>::Scalar, Dynamic, Dynamic>, 0, DynamicStride >((typename internal::traits<OtherDerived>::Scalar *)source.data(), source.rows(), source.cols(), DynamicStride(source.outerStride(), source.innerStride()));
}

typedef Map< MatrixXd, 0, DynamicStride > MapXd;

should allow me to convert any dense object to a Map, and it at least seems to work for Matrices, Vectors, Maps and Blocks.

There still remains the problem of using such a map as an lvalue: Some of my subroutines need to write to their arguments, so I pass those as a reference to MapXd (e.g. "void writeData(MapXd & output)" rather than "void writeData(const MapXd & output)"). Then, however, calling "writeData(toMap(myMatrix))" fails because "toMap(myMatrix)" can't be used as an lvalue (error code "Non-const lvalue reference to type 'MapXd'..." in Clang). So I have to write

Code: Select all

MapXd tmp = toMap(myMatrix);
writeData(tmp);

instead. Is there any way to avoid this?

Another question:

In some of my modules, the code

Code: Select all

MatrixXd A;
   VectorXd b;
   A.row(0) = b;

gives the following errors (with Clang as well as with LLVM GCC 4.2):

Code: Select all

Use of overloaded operator '=' is ambiguous (with operand types 'RowXpr' (aka 'Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, internal::traits<Matrix<double, -1, -1, 0, -1, -1> >::ColsAtCompileTime, IsRowMajor>') and 'VectorXd' (aka 'Matrix<double, Dynamic, 1>'))

with candidate functions

Code: Select all

template <typename OtherDerived>
    Derived& operator=(const DenseBase<OtherDerived>& other);

    template <typename OtherDerived>
    Derived& operator=(const EigenBase<OtherDerived>& other);

from MatrixBase. In other modules, this error doesn't appear. Have you got any idea what could cause the issue?

For your first question: that's a general C++ issue, and there is no ideal solution. With C++11 new features you can probably use rvalue references and/or move semantic for that purpose.

A simple solution is to pass the Map<> object by value rather than by reference. Map objects are very lightweight.

Another solution is to use a 'const Map<Type>&' and const cast the map object. Note that in Eigen a Map<Type> is different to a Map<const Type>. So if you update your toMap wrapper to properly take care of constness you should be fine. Your const arguments should then be declared using 'const Map<const Type>&'.

I encountered exactly the same abiguous operator= with Clang 3.1 and col instead of row like MrMage did (see below). Is there any workaround to fix this?

MrMage wrote:Another question:

In some of my modules, the code

Code: Select all

MatrixXd A;
   VectorXd b;
   A.row(0) = b;

gives the following errors (with Clang as well as with LLVM GCC 4.2):

Code: Select all

Use of overloaded operator '=' is ambiguous (with operand types 'RowXpr' (aka 'Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, internal::traits<Matrix<double, -1, -1, 0, -1, -1> >::ColsAtCompileTime, IsRowMajor>') and 'VectorXd' (aka 'Matrix<double, Dynamic, 1>'))

with candidate functions

Code: Select all

template <typename OtherDerived>
    Derived& operator=(const DenseBase<OtherDerived>& other);

    template <typename OtherDerived>
    Derived& operator=(const EigenBase<OtherDerived>& other);

from MatrixBase. In other modules, this error doesn't appear. Have you got any idea what could cause the issue?

sbirk wrote:I encountered exactly the same abiguous operator= with Clang 3.1 and col instead of row like MrMage did (see below). Is there any workaround to fix this?

My solution was to comment out in MatrixBase.h:

Code: Select all

    /*template <typename OtherDerived>
    Derived& operator=(const EigenBase<OtherDerived>& other);*/


and modify Matrix.h:

Code: Select all

     template<typename OtherDerived>
     EIGEN_STRONG_INLINE Matrix& operator=(const EigenBase<OtherDerived> &other)
     {
-      return Base::operator=(other);
+      other.derived().evalTo(derived());
+      return derived();
     }


as well as Assign.h:

Code: Select all

-template<typename Derived>
+/*template<typename Derived>
 template <typename OtherDerived>
 EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::operator=(const EigenBase<OtherDerived>& other)
 {
   other.derived().evalTo(derived());
   return derived();
-}
+}*/


It compiles well and runs flawlessly for me, but I don't know which side effects this could have.

For completeness, here's the final code arrived at in order to obtain maps from any kind of Eigen objects:

Code: Select all

typedef Stride< Dynamic, Dynamic > DynamicStride;

template <typename OtherDerived>
inline Map< Matrix< typename internal::traits<OtherDerived>::Scalar, Dynamic, Dynamic>, 0, DynamicStride > toMap(const OtherDerived &source)
{
   return Map< Matrix<typename internal::traits<OtherDerived>::Scalar, Dynamic, Dynamic>, 0, DynamicStride >
   ((typename internal::traits<OtherDerived>::Scalar *)source.data(),
    source.rows(),
    source.cols(),
    DynamicStride((OtherDerived::RowsAtCompileTime != 1 ? source.outerStride() : source.innerStride()),
               (OtherDerived::RowsAtCompileTime != 1 ? source.innerStride() : source.outerStride())));
}

template <typename OtherDerived>
inline Map< Matrix< typename internal::traits<OtherDerived>::Scalar, Dynamic, 1>, 0, DynamicStride > toVectorMap(const OtherDerived &source)
{
   assert(source.rows() == 1 || source.cols() == 1);
   bool swapStrides = (source.rows() == 1
                  && !OtherDerived::IsRowMajor);
   return Map< Matrix<typename internal::traits<OtherDerived>::Scalar, Dynamic, 1>, 0, DynamicStride >
   ((typename internal::traits<OtherDerived>::Scalar *)source.data(),
    source.size(),
    DynamicStride((swapStrides != 1 ? source.outerStride() : source.innerStride()),
               (swapStrides != 1 ? source.innerStride() : source.outerStride())));
}

typedef Map< MatrixXd, 0, DynamicStride > MapXd;
typedef Map< VectorXd, 0, DynamicStride > VectorMapXd;

sbirk: I cannot reproduce. Here is my test program:

Code: Select all

#include <Eigen/Dense>
using namespace Eigen;

int main()
{
  MatrixXd A(10,10);
  VectorXd b(10);
     A.row(0) = b;
}


ggael wrote:sbirk: I cannot reproduce. Here is my test program:

I could only reproduce the issue within templated classes with Clang. Here's a test case that doesn't compile for me:

Code: Select all

#include <Eigen/Dense>
using namespace Eigen;

template< typename TestType >
class TemplateTest
{
public:
   TemplateTest() { }
   
   void otherTest()
   {
      MatrixXd A(4, 4);
      VectorXd b(4);
      A.row(0) = b;
   }
};

int main()
{
    TemplateTest<double>().otherTest();
}

With error:

Code: Select all

.../eigenTests.cpp:22:12: Use of overloaded operator '=' is ambiguous (with operand types 'RowXpr' (aka 'Block<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 1, internal::traits<Matrix<double, -1, -1, 0, -1, -1> >::ColsAtCompileTime, IsRowMajor>') and 'VectorXd' (aka 'Matrix<double, Dynamic, 1>'))

ggael wrote:sbirk: I cannot reproduce. Here is my test program:

Sorry, I should have supplied you with an example in the first place. The code below is a stripped down example of the code causing the ambiguity.

Compiling with

Code: Select all

g++ -I../../libs/eigen -c eigencol.cpp -o eigencol.o

works well, but

Code: Select all

clang++ -I../../libs/eigen -c eigencol.cpp -o eigencol.o

fails.

Here is the code:

Code: Select all

#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <Eigen/Core>

using namespace std;

typedef Eigen::MatrixXcd  MATDENSE;
typedef Eigen::SparseMatrix<std::complex<double>, Eigen::ColMajor>  MATSPARSE;

template <typename MAT>
void testcall_impl( MAT &M )
{
    MATDENSE A(3,3);
    Eigen::VectorXcd v(3);
    A.col(2) = v;
}

void testcall(MATSPARSE &M)
{
    testcall_impl(M);
}


MrMage wrote:My solution was to comment out in MatrixBase.h: (...)

Thanks. This seems to be working well. I just used

Code: Select all

#ifndef __clang__

, s.t. the gcc version isn't changed at all.

ok, now I can reproduce. However the workaround will be more complicated.

We should also report to clang because this is clearly a compiler bug. Any volunteer?

Maybe we could narrow down the test case so that it doesn't have to include the whole Eigen library in order to fail? I would then file a bug with Apple (in fact, the more people file the bug, the quicker it gets fixed).

I agree, that for this bug report the Eigen library should be condensed to only the necessary parts to illustrate the bug. But I'm afraid that I don't know the internals well enough to do this.

ggael wrote:ok, now I can reproduce. However the workaround will be more complicated.

Any update on a proper workaround? My workaround leads to an error "/Users/daniel/Documents/Xcode/Uni/eigen/Eigen/src/Core/DenseBase.h:508:65: No member named 'THE_EVAL_EVALTO_FUNCTION_SHOULD_NEVER_BE_CALLED_FOR_DENSE_OBJECTS' in 'Eigen::internal::static_assertion<false>'" with the following code:

Code: Select all

VectorXd u, f;
ConjugateGradient< SparseMatrix<double> > solver(C);
u = solver.solve(f);

I would just like to add that I am also experiencing this problem. Was it reported to Apple/Clang? I don't understand Eigen's internals well enough to file a meaningful bug report myself.