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

[ Thread Index | Date Index | More Archives ]

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,

On 09/21/2012 05:32:25 PM, Gael Guennebaud wrote:

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)),

  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> >.


Mail converted by MHonArc 2.6.19+