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

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

Hi,

that's an interesting example. Here is one possible way to deal with Ref<>:

VectorXd Block_Solve(Ref<MatrixXd> A, const Ref<const VectorXd>& b )
{
  int   N = A.rows(),  m = Teil_Dim, mrest= N-m;
  typedef typename VectorXd::SegmentReturnType PartialVector;


  if ( N <= m ) return  PartialPivLU<MatrixXd>(A).solve(b);
  VectorXd x = b;

  Ref<MatrixXd>  A11(A.block(0,0,m,m)),
                 A12(A.block(0,m,m,mrest)),
                 A21(A.block(m,0,mrest,m)),
                 A22(A.block(m,m,mrest,mrest));

  PartialVector x1(x,0,m), x2(x,m,mrest);
/*
  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 );
  x1 = A11_LU.solve( VectorXd(x1 - A12*x2) );

  return x;
}

By default Ref<MatrixXd> can map any MatrixXd or Block<MatrixXd> or
Map<MatrixXd> as long as the "inner" stride (or increment) is 1.

You can also declare the A11 and the likes as Eigen::Block<Ref<MatrixXd> >.

gael

On Fri, Sep 21, 2012 at 1:42 PM, Helmut Jarausch
<jarausch@xxxxxxxxxxxxxxxxxxx> wrote:
> #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;
> }