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

[Helmut Jarausch <jarausch@xxxxxxxxxxxxxxxxxxx>]

Hi Gael,

I'm glad you showed me such an elegant and simple looking solution.
This can even be shown to my Mathematics students with yet a limited  
experience in C++.
It's the beauty of C++ that many solutions look pretty simple while the  
implementation
is just the opposite ... unless one doesn't make any errors.

Many, many thanks,
Helmut.


On 09/21/2012 05:32:25 PM, Gael Guennebaud 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