Hi all,
I would like to construct a named LU decomposition and store this in a user defined class. For example,

file - solver.h
------------------------------------------------------
class solver
{
public:
solver(double dt, spline2D &spline);
~solver();

// Solver Variables
VectorXi BC;
MatrixXd f,MM,KK,VV;
!!!!! LU_decomp_of_MM !!!!!!
double dt,xMin,xMax,yMin,yMax,lx,ly;
std::ofstream gnuOut;

// Equation Variables
int neq,Vflag;
double mass,bx,by,EtransX,EtransY,x0,y0;

// Methods
double potential1D(double x);
double potential2D(double x, double y);
double getFreeWave1D(double x, double y, double time);

void readEqn();
void getInitialConditions1D(spline1D &spline);
void getInitialConditions2D(spline2D &spline);
void mkv(spline2D &spline);
void writeGnuForm(spline2D &spline, double time);
void getF(MatrixXd &f, MatrixXd &F);
void eulerStep();
void rk4Step();

private:
int p, q, ncp, ncp_l, ncp_x, ncp_y;

};
---------------------------------------------------------

I have no idea how to do this. Any help/direction would be vastly appreciated.

Thanks in advance,
-Corye

With Eigen 2.0:

Code: Select all

LU<MatrixXd> lu;


Then you can do:

Code: Select all

lu.compute(some_matrix);


Now your lu object stores the lu decomposition of the given matrix.

Thank you for quick reply.

I tried the following:
---------------------------------------
file solver.h
class solver
{
...
VectorXi BC;
MatrixXd f,MM,KK,VV;
LU<MatrixXd> lu;
...
};
----------------------------------------
file solver.cpp
void solver::mkv(spline2D &spline)
{
...
[define A]
lu.compute(A);
}

----------------------------------------
file tdse.cpp
voic solver::tdse()
{
...
[use lu]
...
}
----------------------------------------

And I get:
solver.h:20: error: ISO C++ forbids declaration of âLUâ with no type
solver.h:20: error: expected â;â before â<â token
solver.cpp: In member function âvoid solver::mkv(spline2D&)â:
solver.cpp:158: error: âluâ was not declared in this scope
make: *** [objs/solver.o] Error 1


My apologies if I am confused.

Did you check namespaces (Eigen::)? - ok, else MatrixXd won't work.
Did you include Eigen's LU header? - most algorithms are separate from Core.

regards
Sebastian

#ifndef SOLVER_H
FILE::solver.h

#include"main.h"
#include <Eigen/LU>// provides LU decomposition
USING_PART_OF_NAMESPACE_EIGEN
using namespace Eigen;

//*********************************************************************//
//*********************************************************************//

class solver
{
public:
...
LU<MatrixXf> lu;
...
}

FILE solver.cpp
//********************************************************************************//
//********************************************************************************//
// Method for calculating mass, stiffness and potential matrices
void solver::mkv(spline2D &spline)
{
...
lu.compute(LHS);
...
}

-------------------------------------------------------
-------------------------------------------------------

I now get a new error, so I must be getting close!

g++ -I ../EIGEN/ -c solver.cpp -o objs/solver.o
solver.cpp: In constructor âsolver::solver(double, double, double, double, double, spline2D&)â:
solver.cpp:25: error: no matching function for call to âEigen::LU<Eigen::Matrix<float, 10000, 10000, 2, 10000, 10000> >::LU()â
../EIGEN/Eigen/src/LU/LU.h:331: note: candidates are: Eigen::LU<MatrixType>::LU(const MatrixType&) [with MatrixType = Eigen::Matrix<float, 10000, 10000, 2, 10000, 10000>]
../EIGEN/Eigen/src/LU/LU.h:61: note: Eigen::LU<Eigen::Matrix<float, 10000, 10000, 2, 10000, 10000> >::LU(const Eigen::LU<Eigen::Matrix<float, 10000, 10000, 2, 10000, 10000> >&)
solver.cpp: In member function âvoid solver::mkv(spline2D&)â:
solver.cpp:163: error: âclass Eigen::LU<Eigen::Matrix<float, 10000, 10000, 2, 10000, 10000> >â has no member named âcomputeâ
make: *** [objs/solver.o] Error 1

I'm sorry, the default constructor LU(void) only exists in the development branch, not in 2.0. I had forgotten this "detail"...

Either you switch to the development branch (then, note that LU was renamed to FullPivLU) or you take a different approach so you don't need a default constructor.

I understand.

Are there cons to using the devel branch? I am a new Eigen user and need reliable results quickly w/o bugs.

By different approach, do you mean a code rewrite? or another approach w/in Eigen?

Thanks Jacob!

The devel branch can be relied on if you make sure to run the test suite before:
http://eigen.tuxfamily.org/index.php?title=Tests

Otherwise it's like a gamble, some revisions are bad.

By a different approach I mean any approach that does without using a default LU() constructor. For example your class could store a pointer

Code: Select all

LU<MatrixXd> *lu;


and when you have your matrix to decompose, you do lu = new LU<MatrixXd>(matrix);

Works fine now. Thank you for all your help guys. You have been great.
-Corey

I spoke prematurely. The code compiles using

LU<MatrixXd> *lu; in class header file
and
lu = new LU<MatrixXd>(matrix); in method

but,
lu.rank();
in an alternate method give the error::
tdse.cpp:146: error: request for member ârankâ in â((solver*)this)->solver::luâ, which is of non-class type âEigen::LU<Eigen::Matrix<double, 10000, 10000, 2, 10000, 10000> >*â

hey, this is a general C++ question, not Eigen specific

lu->rank() since lu is a pointer.

And make sure to 'delete' every pointer you have gotten through 'new'.

Ah, I see ... you must use
lu->rank();

Yep sorry, answered my own question ...