Me again ... sorry for spamming the board.

I am using the SparseLDLT solver, but the following is true for triangular solvers (of course the LDLT solver is implemented using the triangular solvers).

I have a dense matrix and each column needs a different LDLT solver to act on it (there are alternatives which would not require this, but I think this will be the fastest). Therefore I thought I would do the following:

Code: Select all

 for ( ... ) {

    if ( ! Solvers[ k ].solveInPlace( field.col( k ) ) ) {
     /* do some stuff for failure */
     ...
    };
  };


Solvers is a std::vector<SparseLDLT>, which only has a solveInPlace, no solve -- I guess this is because the return value is a success flag. Anyway, I would think this would be easy, since treating a column of a column major matrix is easy.

However, field.col( ) returns a Block, which is not derived from MatrixBase.

Would it make sense for Block to be derived from MatrixBase? It seems like one should be able to treat a block of a matrix as a matrix.

There are some possibilities:
1) Copy field.col( k ) into a temporary, solveInPlace(), then copy temporary into result.
2) Add a solve() to SparseLDLT, but then there is no way to report failure, except throwing an error.
3) Make a solveInPlace which take Block, but this would need to be done up the chain to the triangular solvers
4) Make a Block derived from MatrixBase? This probably breaks a lot of things, right?

Any suggestions?

Thanks,
Leo

Block does inherit MatrixBase. Your problem you hit here is a limitation of C++ which does not allow passing non const reference to temporary. The solution is to name the temporary:

type_of_field::ColXpr col_k(field, k);
bla.solveInPlace(col_k);

In Eigen3 you can implicitly solves in place as follow:

field.col(k) = vla.solve(field.col(k));

where we autoamagically detect both "b" and "x" are the same...

Thanks, Gael. I had several errors in my thinking. I didn't notice that row() and col() have specialized return types different from block(), and I didn't realize that all of these are derived from MatrixBase. And I never knew about this C++ issue. Does this mean that the compiler was choosing the const version of col(), because it was a temporary?

You learn something every day ...

duetosymmetry wrote: Does this mean that the compiler was choosing the const version of col(), because it was a temporary?

yes, that's exactly it.

The rationale is that, at a first glance, it does not make sense to modify a temporary, but when a temporary is basically a reference to a non temporary (e.g., Block), then it makes sense to modify the referenced object. Future C++ standard solve this issue via rvalue references.