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

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

```

 Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/