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

[Helmut Jarausch <jarausch@xxxxxxxxxxxxxxxxxxx>]

Hi,

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


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

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