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

[Helmut Jarausch <jarausch@xxxxxxxxxxxxxxxxxxx>]

On 09/21/2012 10:56:29 AM, Christoph Hertzberg wrote:
> On 21.09.2012 10:46, Helmut Jarausch wrote:
> > 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.
>
> That's essentially the same problem Norman Goldstein faces at the  
> moment (see the thread "Recursion and block matrices").
> The easiest solution would be to use the Ref-class from the  
> dev-branch (see Gael's mails in that thread), otherwise you need to  
> handcraft some kind of Map object.

Unfortunately, not easy enough for me :

#include <iomanip>
using std::hex;

#include <Eigen/Core>
#include <Eigen/LU>

using Eigen::MatrixXd; using Eigen::VectorXd;
using Eigen::Dynamic;  using Eigen::PartialPivLU;
using Eigen::ArrayXd;  using Eigen::MatrixBase;
using Eigen::Ref;

#define DBG 1

const int Teil_Dim = 4;

VectorXd Block_Solve(Ref<MatrixXd> A, const Ref<VectorXd>& b);
VectorXd Block_Solve(Ref<MatrixXd,0,Eigen::Stride<Dynamic,Dynamic> > A,  
const Ref<VectorXd>& b);

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

  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);
// error: no matching function for call to  
'Eigen::VectorBlock<Eigen::Ref<Eigen::Matrix<double,
//        -0x00000000000000001, 1> >,  
-0x00000000000000001>::VectorBlock(Eigen::VectorXd&, int, int&)'
/*
  x2  -= A21 * (  x1 / A11 );
  A22 -= A21 * ( A12 / A11 );
*/
  PartialPivLU<MatrixXd> A11_LU(A11);
  x2 -= A21 * A11_LU.solve(x1);
  A22-= A21 * A11_LU.solve(A12);

  x2 = Block_Solve( A22, x2 );
// error: call of overloaded 'Block_Solve(Block_Solve_Impl(TypeOfA&,  
const TypeOfb&)
//[with TypeOfA = Eigen::Ref<Eigen::Matrix<double,  
-0x00000000000000001, -0x00000000000000001> >,
// TypeOfb = Eigen::Ref<Eigen::Matrix<double, -0x00000000000000001,  
1> >,
// Eigen::VectorXd = Eigen::Matrix<double, -0x00000000000000001,  
1>]::PartialMatrix&,
// Block_Solve_Impl(TypeOfA&, const TypeOfb&)
// [with TypeOfA = Eigen::Ref<Eigen::Matrix<double,  
-0x00000000000000001, -0x00000000000000001> >,
// TypeOfb = Eigen::Ref<Eigen::Matrix<double, -0x00000000000000001,  
1> >,
// Eigen::VectorXd = Eigen::Matrix<double, -0x00000000000000001,  
1>]::PartialVector&)' is ambiguous

//  x1= A11_LU.solve( VectorXd(x1) - A12*x2 );  // aliasing problem !!!
  x1= A11_LU.solve( VectorXd(x1 - A12*x2) );  // aliasing problem !!!

  return x;
}

VectorXd Block_Solve(Ref<MatrixXd> A, const Ref<VectorXd>& b) {
  return Block_Solve_Impl(A,b);
// error: no matching function for call to  
'Eigen::VectorBlock<Eigen::Ref<Eigen::Matrix<double,
//       -0x00000000000000001, 1> >,  
-0x00000000000000001>::VectorBlock(Eigen::VectorXd&, int, int&)
}

VectorXd Block_Solve(Ref<MatrixXd,0,Eigen::Stride<Dynamic,Dynamic> > A,  
const Ref<VectorXd>& b) {
  return Block_Solve_Impl(A,b);
}


/* Teste Block_Solve an Hand einer (transponierten)
   Vandermonde Matrix mit Stuetzstellen an den Nullstellen
   der Tschebyscheff-Polynome
*/

int main()
{  int N;
   int i;
   cerr << "Dimension: "; cin >> N;

   MatrixXd  A(N,N);
   VectorXd  rhs(N);  // wird anfangs als temp. Vektor benutzt
   VectorXd  x(N);

   for (i= 0; i < N; i++ ) rhs[i]= i;
   rhs= cos( (M_PI/N)*rhs.array() );  // Nullst. T-Polynome
   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; //vorgegebene Loesung
   rhs= A*x;

   VectorXd  Lsg = Block_Solve( A, rhs );
// error: call of overloaded 'Block_Solve(Eigen::MatrixXd&,  
Eigen::VectorXd&)' is ambiguous
   cout << "Fehler: " << (Lsg-x).norm() << endl;
   return 0;
}