Hi,

first thank you for the library. I was looking for the library that could simplify my calculation of FEM method in electromagnetics. Till now I have been working on my procedures for calculations of sparse matrices and solvers (BiCG with SSOR precondition). The problem is I wrote relatively good solver that could easily solve 1mio by 1mio size matrix. My methods were relative stable for my cases (type of numbers complex<double>), however I could never make the true parallization (SSE) of my methods (main problem was slow SSOR precondition). I have decided to use your library for my calculation and I have found few things:

1. I am not able to get good convergence for relative small example 2500x 2500 matrix. I can get result/convergence only with diagonal precondition, while ILU precondition produce the divergence (in small example 5 by 5 it works great). Could this be due to pivoting?

2. My own BiCG method I tried to implement in the Eigen was not working as it should work like in my software (it is true I havent looked in details of your library to check if I am using wrong function calls, I have been using the library only for a week). I have decided to modified a little bit the stabilized BICG method you have implemented (according to http://www.netlib.org/templates/templates.pdf) and now I am getting a better convergence, at least for one case I have achieved the convergence while with your version of solver I could not produce the result.

I have changed the lines in BiCGSTAB.h:
kt = precond.solve(t);
ks = precond.solve(s);
w = kt.dot(ks) / kt.squaredNorm();
to:
w = t.dot(s) / t.squaredNorm();

3. I am aware the linear solvers are still in development or unstable phase, I just want to get some comment or I can share my experience and maybe even provide the test files for checking convergence of the methods, ….
4. I have also “bad” feeling the calls for some procedures: matrix (I am using RowMajor) vector multiplication or other procedures may still have some bugs which are than translated to BiCG method.
5. I have also tried to produce the SSOR precondition with omega = 1, (or SGS precondition). For instance I have changed the line:
z = precond.solve(s);
to:
z1=mat.template triangularView<UnitLower>().solve(s);
z=mat.template triangularView<Upper>().solve(z1);
However, I am not getting any good results. I am calling the procedures in wrong way? Any comment? Maybe I dont understand completely the calls or I am not aware what is happening behind these calls? Or maybe again we have pivoting problems?


Thank you for any comment or useful hint that could help me understand the library, …
I am not a good programmer, especially I am not so familiar with templates and is hard for me to follow the code in the library.

Best regards
Andrej

Thanks for your feedbacks.

1-
All these algorithms are very news and this is still a work in progress. Typically, the ILU preconditioner is very naive and does not do any pivoting, thus explaining the divergence. It is currently only first step towards a "true" ILU preconditioner with thresholding and pivoting. Currently, it's probably not recommended to use it, as you noticed !

2-
Did you observe this behavior with the diagonal preconditioner or with the naive ILU? Because your change is mostly discarding some preconditioning steps. If that's only with ILU, then what happens on the number of iterations with the diagonal preconditioner with and without your change?

3-
You're very welcome. The sparse module is my top first priority for Eigen at the moment.

4-
I'd doubt that, we have quite extensive unit tests.

5-
I think that's rather:

z = mat.template triangularView<Lower>().solve(s);
z = mat.diagonal() * z;
z = mat.template triangularView<Upper>().solve(z);

Alternatively, you could pre scale the columns of the lower part by the diagonal coefficients and then use your simplified version. It is better to store mat.diagonal() into a dense vector once and for all because the extraction is quite costly.


hope that helps.

by curiosity, I tried a SGS preconditioner on random but easy matrices (nice diagonal) and it seems to work very very badly. The naive ILU works well for such easy cases.

Here is the code:

Code: Select all

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.

#ifndef EIGEN_SSORPRECONDITIONER_H
#define EIGEN_SSORPRECONDITIONER_H

template <typename _Scalar>
class SSORPreconditioner
{
    typedef _Scalar Scalar;
    typedef Matrix<Scalar,Dynamic,1> Vector;
    typedef typename Vector::Index Index;
    typedef SparseMatrix<Scalar,ColMajor> FactorType;

  public:
    typedef Matrix<Scalar,Dynamic,Dynamic> MatrixType;

    SSORPreconditioner() : m_isInitialized(false) {}

    template<typename MatrixType>
    SSORPreconditioner(const MatrixType& mat) : m_isInitialized(false)
    {
      compute(mat);
    }

    Index rows() const { return m_lu.rows(); }
    Index cols() const { return m_lu.cols(); }

    template<typename MatrixType>
    SSORPreconditioner& compute(const MatrixType& mat)
    {
      m_lu = mat;
      for(int j=0; j<m_lu.outerSize(); ++j)
      {
        Scalar d = 0;
        typename FactorType::InnerIterator it(m_lu,j);
        while(it && it.index()<j)
          ++it;
        if(it && it.index()==j)
        {
          d = 1/it.value();
          ++it;
        }
        for(; it; ++it)
        {
          eigen_assert(it.index()>j);
          it.valueRef() *= d;
        }
      }

      m_isInitialized = true;
      return *this;
    }

    template<typename Rhs, typename Dest>
    void _solve(const Rhs& b, Dest& x) const
    {
      x = m_lu.template triangularView<UnitLower>().solve(b);
      x = m_lu.template triangularView<Upper>().solve(x);
    }

    template<typename Rhs> inline const internal::solve_retval<SSORPreconditioner, Rhs>
    solve(const MatrixBase<Rhs>& b) const
    {
      eigen_assert(m_isInitialized && "SSORPreconditioner is not initialized.");
      eigen_assert(cols()==b.rows()
                && "SSORPreconditioner::solve(): invalid number of rows of the right hand side matrix b");
      return internal::solve_retval<SSORPreconditioner, Rhs>(*this, b.derived());
    }

  protected:
    FactorType m_lu;
    bool m_isInitialized;
};

namespace internal {

template<typename _MatrixType, typename Rhs>
struct solve_retval<SSORPreconditioner<_MatrixType>, Rhs>
  : solve_retval_base<SSORPreconditioner<_MatrixType>, Rhs>
{
  typedef SSORPreconditioner<_MatrixType> Dec;
  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)

  template<typename Dest> void evalTo(Dest& dst) const
  {
    dec()._solve(rhs(),dst);
  }
};

}

#endif // EIGEN_SSORPRECONDITIONER_H



ggael wrote:2-
Did you observe this behavior with the diagonal preconditioner or with the naive ILU? Because your change is mostly discarding some preconditioning steps. If that's only with ILU, then what happens on the number of iterations with the diagonal preconditioner with and without your change?

I did two tests from my FEM method (matrices RowMajor, type complex<double>)

"My modification of your BiCGstab method" BICGstab version described in the templates pdf

Case A: 3900 iterations
Case B: 3857 iterations

your way:
Case A: 3601 iterations
Case B: 6783 iterations


It looks like both are OK. We could check which method is better when we get better preconditioner.

I have tried some preconditions in my library (not eigen) and best was SSOR or SGS (symetric Gauss Siedl), I never got good results only with Gauss Siedl.

ggael wrote:by curiosity, I tried a SGS preconditioner on random but easy matrices (nice diagonal) and it seems to work very very badly.

If I am right this code is Gauss-Seidel iteration without backward iteration? I never got good results with GS its worse than Jacobi, however with SSOR I always get better results or convergence.

ggael wrote:z = mat.diagonal() * z;

this one is not working for me, strange, cant compile.

Thank you for your answers.

Regards Andrej

Hi one more thing, in the Bicgstab code you have one redundant line, since z and ks are the same!

Code: Select all

   
    z = precond.solve(s);
    t.noalias() = mat * z;

    kt = precond.solve(t);
    ks = precond.solve(s);


from the http://en.wikipedia.org/wiki/Biconjugat ... zed_method we have once K(preconditioner) and once K1 which is probably or it should be left preconditioner?

regards Andrej