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

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]


I'm just learning using Eigen with partial matrices.

In the example below I try to implement a Block-Gauss-Elimination method. It works as shown, but my version needs to copy the parameters of function Block_Solve. I didn't manage to find a suitable type of the parameters which would avoid copying.
(The algorithm wouldn't need copies)

Many thanks for your help,

#include <iostream>
using std::cin;  using std::cout;  using std::cerr;  using std::endl;

#include <Eigen/Core>
#include <Eigen/LU>

using Eigen::MatrixXd; using Eigen::VectorXd;
using Eigen::Dynamic;  using Eigen::PartialPivLU;
using Eigen::ArrayXd;

const int Teil_Dim = 4;
typedef Eigen::Block<MatrixXd,Dynamic,Dynamic> PartialMatrix;
typedef VectorXd::SegmentReturnType PartialVector;

// how to avoid the copying for the paramters A and b ???

VectorXd Block_Solve( MatrixXd A, VectorXd b )
{ int   N = A.rows(),  m = Teil_Dim, mrest= N-m;

  if ( N <= m ) return  PartialPivLU<MatrixXd>(A).solve(b);
  VectorXd x = b;

  PartialMatrix  A11(A,0,0,m,m),

  PartialVector x1(x,0,m), x2(x,m,mrest);
  PartialPivLU<MatrixXd> A11_LU(A11);

  x2 -= A21 * A11_LU.solve(x1);
  A22-= A21 * A11_LU.solve(A12);

  x2 = Block_Solve( A22, x2 );

  VectorXd x1d = x1;  // copy needed against aliasing  ???
  x1= A11_LU.solve( x1d - A12*x2 );  // aliasing problem !!!

  return x;

int main()
{  int N;
   int i;
   cerr << "Dimension: "; cin >> N;

   MatrixXd  A(N,N);
   VectorXd  rhs(N);  // temp. Vektor initially
   VectorXd  x(N);

   for (i= 0; i < N; i++ ) rhs[i]= i;
   rhs= cos( (M_PI/N)*rhs.array() );
   ArrayXd   rhs_pow = rhs.array();

   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; // solution
   rhs= A*x;

   VectorXd  Lsg = Block_Solve( A, rhs );
   cout << "Fehler: " << (Lsg-x).norm() << endl;
   return 0;

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