HarmonicSmoothing code uses CG for non-symmetric matrix

[Prof Garth Wells]

Eigen backend fails for HarmonicSmoothing test because HarmonicSmoothing code uses CG for non-symmetric matrix.

[Jan Blechta]

CG is justifiable here because right-hand side is in the image of symmetric part of the operator and initial approximation is in the kernel of non-symmetric part of the operator. This may fail at least from two reasons:

1. the reasoning is perfect under precise arithmetics but may fail in floating-point
2. Eigen uses zero initial guess

I'll look into it.

[Prof Garth Wells]

It fails because Eigen only uses the lower or upper part of the matrix.

Using CG for non-symmetric matrices is unsafe because backends may exploit symmetry for performance reasons.

[Jan Blechta]

Pushed some fixes to branch:jan/feature-424-eigen-backend.

I'd rather say that such backends exploit wrong assumptions. CG is in this case mathematically correct and PETSc respects it. The same trick is actually employed in navier-stokes demo; it breaks also with Eigen CG. Eigen should have Whole apart from Upper and Lower.

[Jan Blechta]

I was looking today if I can propose a fix to Eigen. Little problem is that I cannot compile DOLFIN with Eigen dev version

In file included from /usr/users/blechta/fenics/eigen/Eigen/SparseCore:37:0,
                 from /usr/users/blechta/fenics/eigen/Eigen/Sparse:19,
                 from /usr/users/blechta/fenics/fenics-dev-dbg/src/dolfin/dolfin/la/EigenMatrix.h:31,
                 from /usr/users/blechta/fenics/fenics-dev-dbg/src/dolfin/dolfin/la/EigenMatrix.cpp:18:
/usr/users/blechta/fenics/eigen/Eigen/src/SparseCore/SparseCwiseUnaryOp.h: In member function 'Derived& Eigen::SparseMatrixBase<Derived>::operator*=(const Scalar&) [with Derived = Eigen::Block<Eigen::SparseMatrix<double, 1>, 1, -0x00000000000000001, true>, Eigen::SparseMatrixBase<Derived>::Scalar = double]':
/usr/users/blechta/fenics/fenics-dev-dbg/src/dolfin/dolfin/la/EigenMatrix.cpp:229:24:   instantiated from here
/usr/users/blechta/fenics/eigen/Eigen/src/SparseCore/SparseCwiseUnaryOp.h:169:42: error: no type named 'InnerIterator' in 'class Eigen::Block<Eigen::SparseMatrix<double, 1>, 1, -0x00000000000000001, true>'
/usr/users/blechta/fenics/eigen/Eigen/src/SparseCore/SparseCwiseUnaryOp.h:169:42: error: no type named 'InnerIterator' in 'class Eigen::Block<Eigen::SparseMatrix<double, 1>, 1, -0x00000000000000001, true>'
/usr/users/blechta/fenics/eigen/Eigen/src/SparseCore/SparseCwiseUnaryOp.h:169:42: error: no type named 'InnerIterator' in 'class Eigen::Block<Eigen::SparseMatrix<double, 1>, 1, -0x00000000000000001, true>'
/usr/users/blechta/fenics/eigen/Eigen/src/SparseCore/SparseCwiseUnaryOp.h:169:42: error: no type named 'InnerIterator' in 'class Eigen::Block<Eigen::SparseMatrix<double, 1>, 1, -0x00000000000000001, true>'

@garth-wells, Do you know whether this "non-symmetric CG" trick works with TPetra? Do HarmonicSmoothing unit test and navier-stokes demo pass?

[Prof Garth Wells]

@blechta I don't know if it works with Epetra. I don't like to the 'trick' - it's bound to cause problems sooner or later.

I don't know what's up with the Eigen dev version. @chris_richardson might know.

[Jan Blechta]

Fix proposed upstream. DOLFIN fixes in jan/feature-424-eigen-backend.

[Chris Richardson]

@blechta - OK, let's make a note in the code, so this can be adopted when the fix is widespread.

[Jan Blechta]

So far Eigen guys seem to be ignorant of mathematics and being reluctant to merge the fix in.

[Prof Garth Wells]

Switched to bicgstab for now.

[Martin Sandve Alnæs]

Removing milestone: 1.6 (automated comment)