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

[Helmut Jarausch <jarausch@xxxxxxxxxxxxxxxxxxx>]

On 09/21/2012 09:53:19 AM, Hauke Heibel wrote:
> 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

I fear, that's too ambitious for me :

Trying

template <typename Derived, typename OtherDerived>
VectorXd Block_Solve(MatrixBase<Derived>& A, const  
MatrixBase<OtherDerived>& b )
{ int   N = A.rows(),  m = Teil_Dim, mrest= N-m;
  typedef typename Eigen::Block<Derived,Dynamic,Dynamic> PartialMatrix;
  typedef typename OtherDerived::SegmentReturnType PartialVector;

.....

seems to get the compiler into endless recursion since I get endless  
error messages.
It probably tries recursive definitions which have to fail.

Helmut.




> 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;
> > }
> >
> >