Getting start with amgcl on Visual Studio 2012

[ZhuonanLin]

Hi Denis,

Thanks for your work about this open source package. I am a student working on solve large dense linear system (i.e. Ax = b where A is dense and large). I wrote my own iterative complex gmres solver with preconditioner and I hope I can get a better performance with this amazing amgcl package. I am new to amgcl and based on my previous experience, to use this header-only package, I only need to let the compiler know the path to include the headers. And also to use amgcl, I also tried to build a boost library. However, when I tried to include headers in the "getting start" example in tutorial, it has several errors:

Code:

#include <amgcl/make_solver.hpp>
#include <amgcl/solver/bicgstab.hpp>
#include <amgcl/amg.hpp>
#include <amgcl/coarsening/smoothed_aggregation.hpp>
#include <amgcl/relaxation/spai0.hpp>
#include <amgcl/adapter/crs_tuple.hpp>

int main(){
return 0;
}

Errors:

My env: Visual Studio 2012 w/ msvc-11.0 Boost 1.69

I wonder what other steps I need or what I did wrong to make amgcl work on my system. I am sorry if this question is too simple since I am really new to these stuff. Any help will be appreciated.

[ddemidov]

First, I don't think amgcl is a good choice for a dense system solution. You should be much better off with a package, specifically created to deal with dense systems (i am not an expert here, but libraries like linpack or magma come to mind).

Second, regarding the compilation errors: it looks like the compiler lacks C++11 support, which amgcl requires. Can you upgrade to a newer version?

[ZhuonanLin]

Thanks! It seems to me that my compiler is not new enough to support some C++ 1x features. I wonder what is the minimum requirement of mscv to work with amgcl? Also for my dense problem. The usual method we use is to use part of the dense coefficient matrix A to do precondition (ie. A = A1 + A2 where A1 is sparse and A2 is the part to make A dense, and we use A1 to do precondition, like ILUt etc., this may not be strict, but piratically it helps a lot in our cases). I wonder if i use amgcl, whether it allows me to provide a function that calculate prod = Ax for built-in gmres/fgmres solver, and meanwhile provide an sparse matrix A1 for AMG preconditioner.

Thanks again for your kind help.

[ddemidov]

I don't develop with MSVC, so I am not really sure what is the minimum requirement. Appveyor build server running automated tests uses Visual Studio 14 2015 (https://ci.appveyor.com/project/ddemidov/amgcl/branch/master#L12), so that should be definitely enough.

[ddemidov]

I wonder if i use amgcl, whether it allows me to provide a function that calculate prod = Ax for built-in gmres/fgmres solver, and meanwhile provide an sparse matrix A1 for AMG preconditioner.

I think that should be possible using this signature:

https://github.com/ddemidov/amgcl/blob/d5990756522d4947f4b6d2efa8b3e138becbb71a/amgcl/make_solver.hpp#L128-L130

Here the make_solver class should be created for the sparse A1 part (this will be used to create the preconditioner), and you should pass your (custom) class for A, which will be used for matrix-vector products in the iterative solver, and for which you would need to provide a specialization of amgcl::backend::spmv_impl:

https://github.com/ddemidov/amgcl/blob/d5990756522d4947f4b6d2efa8b3e138becbb71a/amgcl/backend/interface.hpp#L168-L173

[ZhuonanLin]

Hi Denis,

I downloaded VS 2017 w/ msvc14.1 and boost 1.69. Now it's working. Seems the problem is that though msvc11.0 starts to support C++ 11, it only supports partial functionalities.

Let me take a look at how to work with amgcl in my case. Thanks a lot for your quick reply and kind help.

[ZhuonanLin]

Hi Denis,

Basically the package can work in my system. However, when I tried to use the fgmres solver, I just change the "getting start" example solver from bicgstab to fgmres, however it has a problem when compile.

my code

#include <amgcl/make_solver.hpp>
#include <amgcl/solver/fgmres.hpp>
//#include <amgcl/solver/bicgstab.hpp>
#include <amgcl/amg.hpp>
#include <amgcl/coarsening/smoothed_aggregation.hpp>
#include <amgcl/relaxation/spai0.hpp>
#include <amgcl/adapter/crs_tuple.hpp>
#include <vector>

int main(){
    std::vector<int>    ptr, col;
    std::vector<double> val, rhs;
    int n = poisson(1000, ptr, col, val, rhs);
    typedef amgcl::backend::builtin<double> Backend;
    typedef amgcl::make_solver<
        // Use AMG as preconditioner:
        amgcl::amg<
        Backend,
        amgcl::coarsening::smoothed_aggregation,
        amgcl::relaxation::spai0
        >,
        // And BiCGStab as iterative solver:
        amgcl::solver::fgmres<Backend>
    > Solver;

    Solver solve(std::tie(n, ptr, col, val));
    std::vector<double> x(n, 0.0);
    int    iters;
    double error;
    std::tie(iters, error) = solve(rhs, x);
}

it gives the error:

Error C2079 'amgcl::multi_array<double,2>::strides' uses undefined class 'std::array<int,2>' test_amgcl d:\programs\amgcl\amgcl\util.hpp

at line 234 in util.hpp

                   private:
This Line->        std::array<int, N> strides;
                   std::vector<T>  buf;

[ddemidov]

Thanks for pointing that out. 817ce4e should resolve this.

[ZhuonanLin]

Thanks for your reply. Not quite figure not how to implement self-defined to the make_solver class. Is there any example related to this implementation? I think what I need is to use the sparse part for maker_solver class constructor, then when actually solve the problem, I need provide the matrix-vector product function.

The amgcl-related codes in my application is as following:

1. Include the following headers (from amgcl/examples/solver_complex.cpp), possibly redundant

   
   #include <iostream>
   #include <string>

include <boost/program_options.hpp>

include <boost/property_tree/ptree.hpp>

include <boost/property_tree/json_parser.hpp>

include <boost/range/iterator_range.hpp>

include <amgcl/backend/builtin.hpp>

include <amgcl/value_type/complex.hpp>

include <amgcl/solver/runtime.hpp>

include <amgcl/solver/fgmres.hpp>

include <amgcl/coarsening/runtime.hpp>

include <amgcl/relaxation/runtime.hpp>

include <amgcl/preconditioner/runtime.hpp>

include <amgcl/make_solver.hpp>

include <amgcl/amg.hpp>

include <amgcl/adapter/crs_tuple.hpp>

include <amgcl/io/mm.hpp>

include <amgcl/profiler.hpp>

2. maker_solver template
```C++
typedef double T
typedef amgcl::backend::builtin<std::complex<T>> Backend;
typedef amgcl::make_solver<
    // Use AMG as preconditioner:
    amgcl::amg<
    Backend,
    amgcl::coarsening::smoothed_aggregation,
    amgcl::relaxation::spai0
    >,
    // And BiCGStab as iterative solver:
    amgcl::solver::fgmres<Backend>
> Solver;

3. make_solver object, define the vector using my own sparse_matrix class, this possibly can be done using run-time interface in future

   
   // Make jac->IA, ja and jac_a_complex to vector
   std::vector<int> ja(jac->vtxpr, jac->vtxpr + sizeof jac->vtxpr / sizeof jac->vtxpr[0]);
   std::vector<int> ia(jac->Rowoffset, jac->Rowoffset + sizeof jac->Rowoffset / sizeof jac->Rowoffset[0]);
   std::vector<complex<T>> a(jac_a_complex, sizeof jac_a_complex / sizeof jac_a_complex[0]);
   // constructor with jac

Solver solve(std::tie(size, ia, ja, a));

4.  solve the system with Matrix-vector product function 
 In my previous code (where I write a inhouse complex fgmres solver based on  Y. Saad's book[2003] Algo. 6.9), I have a specific function ( lhs_gmres_solver_complex(w, v) ) that calculates w = Av in each gmres step. In the previous post, you mentioned use 
```C++
/** \note Used in spmv() */
template <class Alpha, class Matrix, class Vector1, class Beta, class Vector2, class Enable = void>
struct spmv_impl {
    typedef typename Matrix::SPMV_NOT_IMPLEMENTED type;
};

to carry out matrix-vector product wrapper. However I only found the function object

int    iters;
double error;
std::tie(iters, error) = solve(A, b_in, x_in);

which takes in A for solver. I wonder how should I use the self-defined lhs function to replace the Matrix-class type A above.

Thanks in advance~

[ddemidov]

Include the following headers (from amgcl/examples/solver_complex.cpp), possibly redundant

These are indeed redundant. You seem to use the compile-time interface (specifying amg components like smoothed_aggregation or spia0 explicitly in template parameters), so you can get rid of the boost dependencies. You will need to remove any runtime.hpp headers though, and replace with the ones you use (e.g. <amgcl/coarsening/smoothed_aggregation.hpp>).

solve the system with Matrix-vector product function In my previous code (where I write a inhouse complex fgmres solver based on Y. Saad's book[2003] Algo. 6.9), I have a specific function ( lhs_gmres_solver_complex(w, v) ) that calculates w = Av in each gmres step.

You can introduce a dummy struct:

struct my_matrix {};

(or use your existing type that carries your dense matrix data).

After that, you will need to specify the amgcl::backend::spmv_impl<> template. An example can be found here. You should put into amgcl::backend namespace (somewhere before your main function is a good place):

namespace amgcl {
namespace backend {

template <class Alpha, class Vector1, class Beta, class Vector2>
struct spmv_impl<Alpha, /* this is your struct: */my_matrix, Vector1, Beta, Vector2> {
  static void apply(Alpha alpha, const my_matrix &A, const Vector1 &x, Beta beta, Vector2 &y)
  {
    // This should implement operation y = beta * y + alpha * A * x
    // you can probably adjust implementation of your `lhs_gmres_solver_complex()` function,
    // but right now you would need a temporary vector to do this:
    std::vector<std::complex<double>> t(y.size());
    lhs_gmres_solver_complex(t, x);
    for (int i = 0; i < y.size(); ++i) {
      y[i] = beta * y[i] + alpha * t[i];
    }
  }

};

} // namespace backend
} // namespace amgcl

You will also need to specialize amgcl::backend::residual_impl<> for your type (see here for an exampe). It should be equivalent to r = rhs - A * x.

Then you could call the solver as

my_matrix A;
std::tie(iters, error) = solve(A, b_in, x_in);

This way the solver will call your implementation of spmv and residual.

[ZhuonanLin]

Thanks so much.

Do I also need to define a

amgcl::backend:clear_impl<> 

function? (I didn't find an similar example like in amgcl/matrix_ops.hpp ) Seems I have compiler errors about this.

Also in the example you write in your last kind respond, is it missing my_matrix in template of spmv_impl

Your example above:

namespace amgcl {
namespace backend {

template <class Alpha, class Vector1, class Beta, class Vector2>
struct spmv_impl<Alpha, /* this is your struct: */my_matrix, Vector1, Beta, Vector2>

...

} // namespace backend
} // namespace amgcl

I wonder if it should be

namespace amgcl {
namespace backend {

template <class Alpha, 
        class my_matrix, 
        class Vector1, class Beta, class Vector2>
struct spmv_impl<Alpha, /* this is your struct: */my_matrix, Vector1, Beta, Vector2>

...

} // namespace backend
} // namespace amgcl

[ddemidov]

Do I also need to define a amgcl::backend:clear_impl<> function?

No, that is for clearing vectors, not matrices. It should work as is with vectors supported by the builtin backend. If you get compile errors, can you post a minimal reproducible example, or at least the compiler output?

Also in the example you write in your last kind respond, is it missing my_matrix in template of spmv_impl

The example is correct as is. You are specializing this template for the case of my_matrix class. See also https://en.cppreference.com/w/cpp/language/template_specialization.

[ZhuonanLin]

Hi Denis,

Thanks so much for reply. The following is the specific way and errors I had in my project:

I have a class named solver, which includes how I calculate the matrix-vector product (public lhs_gmres_solver_complex(Ax, x);), and to use this public function, it needs other attribute and method from the solver object.

To use amgcl as I desire, what I do is to write a function in my solver class, named gmres_amgcl_complex, contains the preconditioner and solver from amgcl:

<solver.cpp>

template<typename T>
int solver::gmres_amgcl_complex(
    int &size_in,
    complex<T> * b_in,
    complex<T> *x_in,
) 
{
        typedef amgcl::backend::builtin<std::complex<T>> Backend;
    typedef amgcl::make_solver<
        // Use AMG as preconditioner:
        amgcl::amg<
        Backend,
        amgcl::coarsening::smoothed_aggregation,
        amgcl::relaxation::spai0
        >,
        // And BiCGStab as iterative solver:
        amgcl::solver::fgmres<Backend>
    > Solver;

        // jac is my sparse matrix object for preconditioner
    //  jac->vtxpr, jac->RowOffset and jac_a_complex are crs matrix arrays from my sprase matrix class
        // Here I just convert them to vector in a naive way for amgcl 
    std::vector<int> ja(jac->vtxpr, jac->vtxpr + sizeof(jac->vtxpr) / sizeof(jac->vtxpr[0]));
    std::vector<int> ia(jac->RowOffset, jac->RowOffset + sizeof(jac->RowOffset) / sizeof(jac->RowOffset[0]));
    std::vector<complex<T>> a(jac_a_complex, jac_a_complex + sizeof(jac_a_complex) / sizeof(jac_a_complex[0]));
    // constructor with jac
    Solver solve(std::tie(size_in, ia, ja, a));

    int    iters;
        double error;
    std::tie(iters, error) = solve(A, b_in, x_in);
}

And then to solve real systems, in my main.cpp, I first define the amgcl::backend::spmv_impl<> and amgcl::backend::residul_impl<> at the beginning of my main.cpp:

<main.cpp>
#include <amgcl/make_solver.hpp>
#include <amgcl/solver/fgmres.hpp>
#include <amgcl/value_type/complex.hpp>
#include <amgcl/amg.hpp>
#include <amgcl/coarsening/smoothed_aggregation.hpp>
#include <amgcl/relaxation/spai0.hpp>
#include <amgcl/adapter/crs_tuple.hpp>

namespace amgcl {
    namespace backend {

        template <class Alpha, class Vector1, class Beta, class Vector2>
        struct spmv_impl<Alpha, /* this is your struct: */LLLG_Eigen_solver_transverse, Vector1, Beta, Vector2> {
            static void apply(Alpha alpha, const LLLG_Eigen_solver_transverse &A, const Vector1 &x, Beta beta, Vector2 &y)
            {
                double omega = A.get_oemga();
                typedef double T;
                complex<T>* v = &x[0];
                complex<T>* lhs = (complex<T>*)malloc(y.size() * sizeof(complex<T>));
                // This should implement operation y = beta * y + alpha * A * x
                // you can probably adjust implementation of your `lhs_gmres_solver_complex()` function,
                // but right now you would need a temporary vector to do this:
                //std::vector<std::complex<T>> t(y.size());
                A.lhs_gmres_solver_complex<T>(lhs, v, omega);
                for (int i = 0; i < y.size(); ++i) {
                    y[i] = beta * y[i] + alpha * lhs[i];
                }
                delete[]lhs;
            }
        };
        template <class LLLG_Eigen_solver_transverse, class Vector1, class Vector2, class Vector3>
        struct residual_impl<
            LLLG_Eigen_solver_transverse, Vector1, Vector2, Vector3
        >
        {
            static void apply(
                Vector1 const &rhs,
                LLLG_Eigen_solver_transverse  const &A,
                Vector2 const &x,
                Vector3       &res
            )
            {
                double omega = A.get_oemga();
                typedef double T;
                complex<T> *x_in = &x[0];
                complex<T>* lhs_complex = (complex<T>*)malloc(rhs.size() * sizeof(complex<T>));
                A.lhs_gmres_solver_complex<T>(lhs_complex, x, omega);
                for (int p = 0; p < y.size(); p++)
                {
                    res[p] = rhs[p] - lhs_complex[p];
                }
                delete[]lhs_complex;
            }
        };
    } // namespace backend
} //namespace amgcl

And Finally I have the computation codes in after the codes above in my int main() function:

int main(){
       solver *A = new solver(input_files); // a new solver object A for computing
       int size_in = A->get_size();
       complex<double> * b_in = (complex<double>*)malloc(size_in *sizeof( complex<double>));
       complex<double> * x_in = (complex<double>*)malloc(size_in *sizeof( complex<double>));
       // and initialization b_in with actual rhs
       A->gmres_amgcl_complex(A->get_size(), b_in, x_in);

      return
}

I have simplified my main function above and in actual one it calls several methods in A to set up the linear system.

And when build the project (compile the codes), I have the following errors:

I wonder if you can kindly help me take a look at when you have time. I really appreciate your help and time!

[ddemidov]

could you please either attach the complete source and header files here (you can drag and drop the files onto the comment field), or upload them to https://gist.github.com and post a link here?

[ddemidov]

But, at the first glance, the compiler mentions my_matrix type, but I can see in the code above that you have specialized the spmv_impl for LLLG_Eigen_solver_transverse. Could that be a problem? Also, I forgot to mention earlier that you either need to put typedef std::complex<T> value_type inside your matrix class, or specialize value_type template (example).

Also, visual studio should have a 'build output' window, which should have more complete compiler output in text format, that you can copy an paste here instead of the screenshots.

[ZhuonanLin]

Hi Denis,

I just make a minimal project that can reproduce the same errors I have in my project. There are only 3 files:

my_solver.h is the header of my solver class my_solver.cpp is the source file of my sovler class amgcl_w_my_class.cpp is my main source file

I have uploaded these file to: https://gist.github.com/ZhuonanLin/477ed73022e53d6b8fbb32075850ad7c

The compile errors given by my machine are: errors.txt

[ddemidov]

There was a bunch of unrelated syntax errors, so I just reimplemented the whole example from scratch. This does compile for me with

g++ -o amgcl_my_problem amgcl_my_problem.cpp -I ~/work/amgcl

Hope this helps:

#include <iostream>

#include <vector>
#include <complex>

#include <amgcl/make_solver.hpp>
#include <amgcl/solver/fgmres.hpp>
#include <amgcl/value_type/complex.hpp>
#include <amgcl/amg.hpp>
#include <amgcl/coarsening/smoothed_aggregation.hpp>
#include <amgcl/relaxation/spai0.hpp>
#include <amgcl/adapter/crs_tuple.hpp>

/* This is just a helper struct that holds your specific problem details,
 * and knows how to do the full spmv operations.
 * There is no need to incapsulate amgcl solver here.
 */
template <typename T>
struct my_problem {
    typedef std::complex<T> value_type;

    int    size;
    double omega;

    my_problem(int size, double omega) : size(size), omega(omega) {}

    int lhs_gmres_complex(
            std::complex<T>* w, const std::complex<T>* v, double omega) const
    {
        for(int i = 0; i < size; ++i) { 
            w[i] = 2. * v[i] - omega * i / 2;
        }

        return 0;
    }
};

/* Now we tell amgcl that my_problem may be treated as a system matrix.
 * That is, we teach amgcl how to perform spmv and residual operations with my_problem.
 */
namespace amgcl { namespace backend {

template <class Alpha, class T, class Vector1, class Beta, class Vector2>
struct spmv_impl<Alpha, my_problem<T>, Vector1, Beta, Vector2> {
    static void apply(Alpha alpha, const my_problem<T> &A, const Vector1 &x, Beta beta, Vector2 &y)
    {
        std::vector<std::complex<T>> t(y.size());

        A.lhs_gmres_complex(&t[0], &x[0], A.omega);
        for (int i = 0; i < y.size(); ++i) {
            y[i] = beta * y[i] + alpha * t[i];
        }
    }
};

template <class T, class Vector1, class Vector2, class Vector3>
struct residual_impl<my_problem<T>, Vector1, Vector2, Vector3> {
    static void apply(
            const Vector1 &rhs,
            const my_problem<T> &A,
            const Vector2 &x,
            Vector3 &res
            )
    {
        A.lhs_gmres_complex(&res[0], &x[0], A.omega);
        for (int i = 0; i < rhs.size(); ++i)
        {
            res[i] = rhs[i] - res[i];
        }
    }
};

} }

int main() {
    int N = 100;

    my_problem<double> test(N, 20.);

    std::vector<std::complex<double>> b(N);
    std::vector<std::complex<double>> x(N);

    // Assemble the sparse part of the matrix:
    std::vector<int> ja(N);
    std::vector<int> ia(N+1);
    std::vector<std::complex<double>> a(N, std::complex<double>(1, 0));

    ia[0] = 0;
    for (int i = 0; i < N; ++i)
    {
        ja[i] = i;
        ia[i + 1] = i;
    }

    // Setup amgcl for the sparse component of the problem:
    typedef amgcl::backend::builtin< std::complex<double> > Backend;
    typedef amgcl::make_solver<
        amgcl::amg<
            Backend,
            amgcl::coarsening::smoothed_aggregation,
            amgcl::relaxation::spai0
            >,
        amgcl::solver::fgmres<Backend>
    > Solver;

    Solver solve(std::tie(N, ia, ja, a));

    int    iters;
    double error;

    // Solve the system for the matrix defined by my_problem (using custom
    // matrix-less spmv and residual operations), and preconditioned by AMG:
    std::tie(iters, error) = solve(test, b, x);

    std::cout << iters << " " << error << std::endl;
}

[ZhuonanLin]

Hi Denis,

Thanks for your help. I was working on something else recently and today I re-work on this problem. With your kind help, I have managed to compiled the amgcl codes with my class and the my codes (with amgcl) give consistent results with my benchmark test.

One other question is that when I tried to solve some complicated system with large number of unknowns ( another benchmark test in my project), the iteration number of amgcl is not as small as I expected. I wonder what configuration I can improve AMG settings to can a better performance (e.g. less iteration number). I read a previous issue post discussing about this and I wonder if there is any rule of thumb to follow, or if there is some literature I can read. I never use AMG before ( usually I use various preconditioners from ILU family) and any suggestion is appreciated!

[ddemidov]

What kind of a problem you are solving (heat transfer, elastisity, navier-stokes)? What is the expected number of iterations, and what did you get with AMG?

You could try to replace amgcl::relaxation::spai0 with amgcl::relaxation::ilu0 if ILU is working better for your problems.

If you are solving a non-scalar problem (with multiple physical unknowns per grid node, e.g. three components of displacement in 3D elasticity problem, or 3 components of velocity plus pressure in Navier-Stokes problem), you can try to either set precond.coarsening.aggr.block_size=<m> where m is the number of unknowns per grid point, or use backend::builtin<amgcl::static_matrix<double, m, m>> instead of backend::builtin<double>.

[ZhuonanLin]

Hi Denis,

Thanks for your reply. My research focus is about micro-magnetism, which tends to solve dynamics of magnetization (generally 3D vector) in material, with the Landau-Lifshitz-Gilbert (LLG) equation as governing equation. The specific problem I am working on is to solve the system with linearized LLG in frequency domain with weak excitation, in finite element regime, which mathematically equivalent to solve large linear system (in our case, the dimension can be up to ~1M ). Based on our previous experience, I am using only the sparse part (A1) of the coefficient matrix (A) as the precondition matrix.

In my codes, I wrote a complex GMRES solver myself based on Prof. Y. Saad 's book, and the preconditioner used as benchmark is actually a direct sparse solver (dss). The dss solves A1 * (Ax) = A1 * b in each step, and this method should give us a lower bound of iteration number that using A1 as preconditioner only can help. However, solving a sparse system as each step will be computationally expensive, so that I am trying to use this AMGCL for help.

After successfully merged amgcl work with my codes, I run some cases with the amgcl setting from the the sample (amgcl::coarsening::smoothed_aggregation, amgcl::relaxation::spai0, and typedef amgcl::backend::builtin< std::complex<double> > Backend;)

1. In first sanity-check test (where i have analytical solution), the system has ~1K nodes (so that 3K unknowns), with same tolerance (1e-8), dss and amgcl can solve the problem in 7 and 9 iterations respectively, which indicates they have similar performance, in term of iteration;

2. A benchmark test, the system has ~30k nodes with ~90K unknowns (In this case actually I have some problems with specific value of omega, the frequency of excitation source, in my solver + dss preconditioner, I will discuss the problem later). For those in-problematic omega (e.g. f = 1Gz), my codes can solve the problem in ~9 iterations, while amgcl with the setting above can NOT solve it in 100 iterations.

I was trying your suggestion in last reply, however, I am not sure if setting typedef amgcl::backend::builtin<amgcl::static_matrix<double, m, m>> can help, since in our data structure, we use m_x1 ... m_xN, m_y1 ... m_yN, m_z1 ... m_zN instead of m_x1, m_y1, m_z1 ... m_xN, m_yN, m_zN. And also in my complex solver, I haven't tried ilu family preconditioner, so I will try to use ilu0 instead of spai0 but I am not sure it can guarantee a better performance.

I also have a general problem in my codes (complex gmres + dss) and I wonder if you have similar previous experience. In my coefficient matrix A, the parameter omega only exists in diagonal terms in imaginary part. At beginning, I thought that a larger omega will also lead to a better convergence, since it help improve the "diagonal dominant" property of A, which is essential in ilu preconditioner family. However, in my actual codes, I observe the following phenomenon. For the benchmark test 2 above, with f = 1 GHz, it can solve with ~9 iteration, however, with f = 25~230 GHz, the convergence is bad and can NOT solve the problem within 100 iterations, then increasing f > 240 GHz, the solver can again solver the problem in only several iterations. We have nearly ruled out the physics reason (e.g. specific omega hits the resonance frequency), and I wonder if you have experience on similar complex GMRES problem. Mathematically or numerically speaking, is there any "competing phenomenon" between real part and imaginary part on diagonal elements? Or is there any way that can help this situation (~30K nodes with >100 iterations is not acceptable in our case)? I have been stuck by this problem for nearly a month and I do really appreciate for any kind suggestion from you.

Thanks and regards, Zhuonan

[ddemidov]

In first sanity-check test (where i have analytical solution), the system has ~1K nodes (so that 3K unknowns), with same tolerance (1e-8), dss and amgcl can solve the problem in 7 and 9 iterations respectively, which indicates they have similar performance, in term of iteration;

By default, amgcl uses sparse direct solver for any grids smaller than 3000 unknowns, so it seems you just reproduced your original approach with dss here.

I was trying your suggestion in last reply, however, I am not sure if setting typedef amgcl::backend::builtin<amgcl::static_matrix<double, m, m>> can help, since in our data structure, we use m_x1 ... m_xN, m_y1 ... m_yN, m_z1 ... m_zN instead of m_x1, m_y1, m_z1 ... m_xN, m_yN, m_zN.

Indeed, blockwise approach (either with static_matrix as value type, or with blockwise aggregation) requires node-wise ordering of the variables. You could in principle try reordering adapter to try this out:

https://github.com/ddemidov/amgcl/blob/master/amgcl/adapter/reorder.hpp

And also in my complex solver, I haven't tried ilu family preconditioner, so I will try to use ilu0 instead of spai0 but I am not sure it can guarantee a better performance.

ilu0 in general has better convergence than spai0 (spai0 is basically a diagonal preconditioner), so it is worth trying.

I also have a general problem in my codes (complex gmres + dss) and I wonder if you have similar previous experience.

Sorry, I have no experience here.

[ZhuonanLin]

Hi Denis,

Thanks for the response. Double checking my code recently...

Just a quick question about the package. I know that the amgcl, together with the vexcl, have a parallel version on GPU with CUDA. I wonder if the gmres/fgmres solver is also paralleled in this situation?

[ddemidov]

I know that the amgcl, together with the vexcl, have a parallel version on GPU with CUDA.

And you can also use cuda backend directly (without vexcl).

I wonder if the gmres/fgmres solver is also paralleled in this situation?

Yes.

[ZhuonanLin]

Hi Denis,

Although I'm still struggling with my convergence problem, I am moving my codes on the basis of the amgcl package. I have several questions I wonder if you can kindly clarify.

As a recap, in my linear problem Ax = b, I have a dense matrix A = A1 + A2, where A1 is sparse and A2 is dense. I use A1 as a left preconditioner matrix for my problem

(1). In my code, I use a direct sparse solver to "solve" the preconditioner matrix, so that in each iteration of fgmres, it first calculates x' = A * x_temp, then slove A1 * x' = (inv(A1) * b) with a direct sparse solver. From the previous post, to make the scheme above as an available option, I want to add a flag in skyline_lu.hpp:

        static size_t coarse_enough() {
            //return 3000 / math::static_rows<value_type>::value;
        if(use_dss){
                return true;
            }
            else{
            return false;
           }
        }

Also I want to make sure that if by doing this, together with the setup we discussed before, I am actually using the scheme as I describe above.

2. In my code I have a flag to print out the error of the current step, in order to monitor the trend of convergence. I wonder if I want to so a similar thing with amgcl, where I should add the output line.

I really appreciate your kind help!

[ddemidov]

From the previous post, to make the scheme above as an available option, I want to add a flag in skyline_lu.hpp

If I understand correctly, you want to have an option to always use direct sparse solver. First, the skyline_lu solver does not scale very well: it will be more and more expensive to setup after something like 10000 unknowns. If you are fine with that, then you could just set precond.coarse_enough option to something large (may be INT_MAX). This way you won't need to modify the code.

I wonder if I want to so a similar thing with amgcl, where I should add the output line.

It depends on the iterative solver you use, but each of those have a convergence check. You could add an output somewhere around that check. For example, fgmres has it here:

https://github.com/ddemidov/amgcl/blob/930993509a9f626e9ee3a13cfbb77431fb77b1d7/amgcl/solver/fgmres.hpp#L202-L205

[ZhuonanLin]

Hi Denis,

Thanks for your reply. RecentlyI am trying to integrate amgcl in my class, and I am not sure if I can do this. For some reason, I really hope that I can and the amgcl in my calss, so that I can have a method in my class (e.g. gmres_amgcl(complex<double> x, complex<double> rhs) ) and inside this function, I use amgcl to solve the linear system. So currently I have my_class.h, my_class.cpp and a main.cpp. In the main.cpp, the only job is to have a new my_class object, and call my_class->gmres_amgcl(complex<double> x, complex<double> rhs to solve the linear system, by passing in the initial guess and the rhs.

I have uploaded my codes with has compile problems, actually in a previous post where I uploaded similar codes, I also tried to do in this way, but I had several other problems, so that when you kindly helped my re-write the code, you just write the amgcl solver in the main function which takes in my_solver object in the solver and I am currently working in this way. I really appreciate your kind help

The link to my codes: https://gist.github.com/ZhuonanLin/9fdad2b5cc9a6766161523cc3db5ae7a

Compiler error I have: error.txt

[ddemidov]

Try to replace this with *this in this line:

https://gist.github.com/ZhuonanLin/9fdad2b5cc9a6766161523cc3db5ae7a#file-my_class-cpp-L95