Good day!

I am trying to solve 1D adv-diff-reaction ODE(exactly, ODE not PDE) with Eigen sparse matrices but I can not get correct solution.

Code: Select all

        VectorXd *mpRhsVec;
        SparseMatrix<double> *mpLhsMat;
for (int i=1; i<mNumNodes-1; i++) {       
        double xm = mpGrid->mNodes[i-1].coordinate;
        double x = mpGrid->mNodes[i].coordinate;
        double xp = mpGrid->mNodes[i+1].coordinate;
        double diffusion_alpha = 2.0/(xp-xm)/(x-xm);
        double diffusion_beta = -2.0/(xp-x)/(x-xm);
        double diffusion_gamma = 2.0/(xp-xm)/(xp-x);
        (*mpLhsMat).insert(i,i-1) = (mpOde->mCoeffOfUxx)*diffusion_alpha - (mpOde->mCoeffOfUx)/(xp-xm);
        (*mpLhsMat).insert(i,i) = (mpOde->mCoeffOfUxx)*diffusion_beta + mpOde->mCoeffOfU;
        (*mpLhsMat).insert(i,i+1) = (mpOde->mCoeffOfUxx)*diffusion_gamma +(mpOde->mCoeffOfUx)/(xp-xm);
}
        .....
       SimplicialLDLT<SparseMatrix<double> > solver;
       solver.compute(*mpLhsMat);
       mpSolVec  = solver.solve(*mpRhsVec);


code could be found here: https://gist.github.com/pycckuu/10640287

I form sparse matrices and XD vectors. I solve the problem but get incorrect solution (I checked with Mathematica, Matlab(used the same matrices and made A\b) as well as I rewrote the same code in PETSC and got correct results). But Eigen give me the same incorrect result for 3 different solvers (CG, SimplicialLDLT, SimplicialLLT). I realize that it is something with how I use it but I can not understand where the problem...
Thank you in advance

By the way, Great toolbox! I would like to use it instead of PETSC due to EIGEN is very intuitive! Great work!

How much incorrect? Are you sure that your problem is symmetric positive definite? Have you checked that solver.info()==Eigen::Success after calling .compute()? Can you try a not too large problem with a dense solver (use MatrixXd(*mpLhsMat) to construct a dense matrix from a sparse one )?

Thank you very much for such quick reply.

How much incorrect?

Completely incorrect.

Are you sure that your problem is symmetric positive definite?

No, I am not sure about this. How can I check it? Is there any way to convert the matrix to symmetric positive definite?

Have you checked that solver.info()==Eigen::Success after calling .compute()?

This test didnt passed I think because of the above reason.

Can you try a not too large problem with a dense solver (use MatrixXd(*mpLhsMat) to construct a dense matrix from a sparse one )?

Thank you for great suggestion. FullPivHouseholderQR worked for me.

indeed, looking at your code, the matrix is not even symmetric. So for large scale problem, you can try SparseLU or BiCGSTAB.

ggael wrote:indeed, looking at your code, the matrix is not even symmetric. So for large scale problem, you can try SparseLU or BiCGSTAB.

I've got BiCGSTAB working but can not force SparseLU to work for me:

Code: Select all

SparseLU<SparseMatrix<double> > slu;
slu.compute(*mpLhsMat);
mpSolVec=slu.solve(*mpRhsVec);

but it tells me:

Code: Select all

Assertion failed: (info && "COLAMD failed "), function operator()


I was trying to find the problem on the forum but didn't have a success. Thank you in advance.

This is because you matrix contains empty rows or columns.

Code: Select all

   SparseLU<SparseMatrix<double,ColMajor>, AMDOrdering<int> > slu;

this solved the problem. But I didn't understand how