Re: [eigen] RFC - how to avoid (impossible) recursive declarations
• To: eigen@xxxxxxxxxxxxxxxxxxx
• Subject: Re: [eigen] RFC - how to avoid (impossible) recursive declarations
• From: Gael Guennebaud <gael.guennebaud@xxxxxxxxx>
• Date: Fri, 21 Sep 2012 17:32:25 +0200
• Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=20120113; h=mime-version:in-reply-to:references:from:date:message-id:subject:to :content-type; bh=lX5dzW75Ac6nd6DXrb55bRh8odTuv/cv28xvTTxrIRk=; b=M1hArCSabsMvaTwhksU4l92VXwbM0x52Hi0QgOEbN/mqKHI8ld559gMxwaAeOIy4o5 zK2q1PfisMYWzq3DnYTaogSwkcwihbJxl1bGrlYORG0/lZKDTtUeM9+wwOr52Jb+f8aB 4haRJXxM5kxK3tM/FRICmhHcEqaZXM02YGJaNe65LhH7ZZZx/D/jBvjYBjcbjV62NMMC 4ktxI+CZ8BVWNqaFBUkt7K6Wl5lMiyjWmMDXpzvPvXjkfByc1kyarvRKm0mZOl/15qfO R2evyOXNXx8d5HsWi2ho0PNxwky08zXHVNH4Ubdvp3Sfx6TjNfOldg4FUYiTFfdy+vve lGBg==

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

```

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