|Re: [eigen] RFC - how to avoid (impossible) recursive declarations|
[ Thread Index |
| More lists.tuxfamily.org/eigen Archives
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
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
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;
PartialVector x1(x,0,m), x2(x,m,mrest);
x2 -= A21 * ( x1 / A11 );
A22 -= A21 * ( A12 / 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) );
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