Re: [eigen] RFC - how to avoid (impossible) recursive declarations

[Hauke Heibel <hauke.heibel@xxxxxxxxx>]

Hi,

first, you should take a look at this part of the documentation:

http://eigen.tuxfamily.org/dox/TopicFunctionTakingEigenTypes.html

Your code compiles because since you pass by value, you allow implicit
conversions. I.e. the function will create temporary objects, as you
already observed. When you look at your code, you want to pass
MatrixXd, VectorXd as well as Eigen::Block<MatrixXd,Dynamic,Dynamic>
and VectorXd::SegmentReturnType. The most specialized common base
class of these objects is MatrixBase. In order to prevent copies and
code duplication you have write Block_Solve as a template function. I
guess the interface would look like this:

template <typename Derived, typename OtherDerived>
VectorXd Block_Solve( const MatrixBase<Derived>& A, const
MatrixBase<OtherDerived>& b );

Regards,
Hauke

On Fri, Sep 21, 2012 at 9:33 AM, Helmut Jarausch
<jarausch@xxxxxxxxxxxxxxxxxxx> wrote:
> Hi,
>
> I'm just learning using Eigen with partial matrices.
>
> In the example below I try to implement a Block-Gauss-Elimination method.
> It works as shown, but my version needs to copy the parameters of function
> Block_Solve.
> I didn't manage to find a suitable type of the parameters which would avoid
> copying.
> (The algorithm wouldn't need copies)
>
> Many thanks for your help,
> Helmut.
>
>
> #include <iostream>
> using std::cin;  using std::cout;  using std::cerr;  using std::endl;
>
> #include <Eigen/Core>
> #include <Eigen/LU>
>
> using Eigen::MatrixXd; using Eigen::VectorXd;
> using Eigen::Dynamic;  using Eigen::PartialPivLU;
> using Eigen::ArrayXd;
>
> const int Teil_Dim = 4;
> typedef Eigen::Block<MatrixXd,Dynamic,Dynamic> PartialMatrix;
> typedef VectorXd::SegmentReturnType PartialVector;
>
> // how to avoid the copying for the paramters A and b ???
>
> VectorXd Block_Solve( MatrixXd A, VectorXd b )
> { int   N = A.rows(),  m = Teil_Dim, mrest= N-m;
>
>   if ( N <= m ) return  PartialPivLU<MatrixXd>(A).solve(b);
>   VectorXd x = b;
>
>   PartialMatrix  A11(A,0,0,m,m),
>                  A12(A,0,m,m,mrest),
>                  A21(A,m,0,mrest,m),
>                  A22(A,m,m,mrest,mrest);
>
>   PartialVector x1(x,0,m), x2(x,m,mrest);
>   PartialPivLU<MatrixXd> A11_LU(A11);
>
>   x2 -= A21 * A11_LU.solve(x1);
>   A22-= A21 * A11_LU.solve(A12);
>
>   x2 = Block_Solve( A22, x2 );
>
>   VectorXd x1d = x1;  // copy needed against aliasing  ???
>   x1= A11_LU.solve( x1d - A12*x2 );  // aliasing problem !!!
>
>   return x;
> }
>
> int main()
> {  int N;
>    int i;
>    cerr << "Dimension: "; cin >> N;
>
>    MatrixXd  A(N,N);
>    VectorXd  rhs(N);  // temp. Vektor initially
>    VectorXd  x(N);
>
>    for (i= 0; i < N; i++ ) rhs[i]= i;
>    rhs= cos( (M_PI/N)*rhs.array() );
>    ArrayXd   rhs_pow = rhs.array();
>
>    A.col(0).fill(1.0);
>    A.col(1)= rhs;
>    for ( i= 2; i < N; i++ ) {
>      rhs_pow*= rhs.array();
>      A.col(i)= rhs_pow;
>
>    }
>    for ( i= 0; i < N; i++ ) x[i]= i+1; // solution
>    rhs= A*x;
>
>    VectorXd  Lsg = Block_Solve( A, rhs );
>    cout << "Fehler: " << (Lsg-x).norm() << endl;
>    return 0;
> }
>
>