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

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

BTW, the low level Map<..> strategy suggested by Hauke or Christoph is
what I did in PartialPivLu.h.

gael

On Fri, Sep 21, 2012 at 5:32 PM, Gael Guennebaud
<gael.guennebaud@xxxxxxxxx> wrote:
> 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;
>> }