Complex-valued problems

christoph-sohrmann commented 8 years ago

Hi.

I would like to solve a complex-valued problem without rewriting it into a real-valued one. Unfortunately, I wasn't successful in using std::complex instead of double. I get a bunch of errors, for a large part related to implicit conversion where you wrote "1/x" instead of "1./x". Also, some comparisons have to be rewritten as "fabs(eps)>fabs(x)". Both can be fixed straightforwardly.

However, I was wondering whether you see more substantial issues within the algorithms when plugging in a complex type?

Thanks, Chris

christoph-sohrmann commented 8 years ago

Update:

In fact, it seems std::complex DOES work as a value type. There are a few things which need to be changed, though: First, write all floating-point literals as "1." instead of "1". Second, rewrite all comparisons as "fabs(x)>fabs(y)". Third, enable the complex vector-type like so:

namespace amgcl {
namespace backend {
  template <>
  struct is_builtin_vector< std::vector<std::complex<double>> > : boost::true_type {};
} // namespace backend
} // namespace amgcl

Et voila, at first sight the results seem pretty promising!

ddemidov commented 8 years ago

That is great, could you make a PR for your changes? Thanks for testing this!

On Sat, Oct 17, 2015, 02:55 dondieselkopf notifications@github.com wrote:

Update:

In fact, it seems std::complex DOES work as a value type. There are a few things which need to be changed, though: First, write all floating-point literals as "1." etc. instead of "1". Second, rewrite all comparisons as "fabs(x)>fabs(y)". Third, enable the complex type like so:

namespace amgcl { namespace backend { template <> struct is_builtin_vector< std::vectorstd::complex > : boost::true_type {}; } // namespace backend } // namespace amgcl

Et voila, at first sight the results seem pretty promising!

— Reply to this email directly or view it on GitHub https://github.com/ddemidov/amgcl/issues/17#issuecomment-148867536.

Cheers, Denis

christoph-sohrmann commented 8 years ago

Will do that. I'm just struggling with the implementation of the test case. Just adding a complex backend to test_solver.cpp is probably not a good idea, since some algorithms seem not to be valid for complex input (e.g. Chebyshev smoother). Now I am trying to write a new testcase with the valid algorithms only and in the style of the existing one (i.e. using the runtime interface). However, the runtime interface seems to statically instantiate all possible algorithms with the complex backend, which I would actually like to avoid for the real-only-algorithms. Do you have a good idea how to go about that? Maybe I should avoid the runtime interface?

ddemidov commented 8 years ago

First, write all floating-point literals as "1." instead of "1"

I think this should be written as static_cast<value_type>(1). 1. is double and won't work when value_type is std::complex<float>.

Third, enable the complex vector-type like so

What I don't like about this is that <complex> would always have to be included. May be its better to move this out into some other header (e.g. <amgcl/backend/enable_complex.hpp>). Also, this should be

template <typename T>
   struct is_builtin_vector< std::vector<std::complex<T>> > : boost::true_type {};

some algorithms seem not to be valid for complex input (e.g. Chebyshev smoother).

If this is only about relaxation, then we can use this template to disable some of the smoothers. Either disable everything by default for complex types, then enable what works, or disable what does not work, whatever is easier. The runtime interface will not try to instantiate unsupported smoothers.

christoph-sohrmann commented 8 years ago

I've already put everything into a separate header file such as enable_complex.hpp, as you suggested. Thanks for the hint with the template relaxation_is_supported. Will try to use that to disable invalid relaxation algorithms. Although I am not sure yet whether it's only them which cause problems. We will see how the test are running through. Guess I need a few more days until sending the PR.

ddemidov commented 8 years ago

No hurry, and thanks for working on this! If you need any help, feel free to ping me.

ddemidov commented 8 years ago

Hi, do you have any news on this? Do you need any help?

christoph-sohrmann commented 8 years ago

Hi. Thanks for aksing. I am sorry for the silence on my part. I'm currently a bit busy otherwise, so I had to put this away for the time-being. Hope I can continue asap. I will keep you posted as soon as there are any news.

ddemidov commented 8 years ago

I may look into this myself then. Do you have a simple complex-valued matrix to test this with?

christoph-sohrmann commented 8 years ago

What I did was to copy your function "sample_problem", changed the template parameter name to "complex" and checked its type with "boost::is_complex::type". Then I have simply used each value for both, the real and the imaginary part of the matrix entry. This should be sufficient for testing purposes.

ddemidov commented 8 years ago

Allowing direct use of complex values through all amgcl components proved to be harder than it seemed initially. Instead, I implemented a complex matrix adapter. It uses the fact that complex-valued system

    Ax = b

may be written as

    [real(A) -imag(A)] [real(x)] = [real(b)]
    [imag(A)  real(A)] [imag(x)]   [imag(b)]

With this, one can transparently use complex-valued matrix to setup scalar valued amg. Ranges for rhs and solution may be used in-place. See tests/test_complex.cpp for example of this.

christoph-sohrmann commented 8 years ago

Those equivalent real formulations (ERF) have a huge draw back, which is the loss of numerical stability. By such a transformation the matrix eigenvalues change in a highly unfavourable manner. For the problems I work on at the moment (electrodynamics) this is not a feasible solution. I was working with ERF but changed to your solver exactly for the reason of stability. Some complex problems can only be solved in the complex domain. With the changes I applied to your code, I have a running version. However, I would need a bit more time to send a pull request, if you are still interested.

http://www.ajur.uni.edu/v1no3/Munankarmy.pdf http://www.math.polytechnique.fr/~tringali/documents/numelec_2006_Angiulli-Tringali-DiMassa.pdf

ddemidov commented 8 years ago

Thanks for pointing that out, I was not aware of this. Another problem is that adapted matrix takes twice as much space, so yes, I am still interested. I tried to do something in complex-values branch. Feel free to look at that, but I am not too proud of what's there.

By the way, this parameter tells aggregation to take block-structure of the matrix into account. It reduced number of iterations by 2-3 times on my simple tests. I am curious, how much would (or would not) that help with your problems.

christoph-sohrmann commented 8 years ago

Taking the block-structure into account should indeed be the right way to cope with it. Nice! I was seeing great improvements with such a strategy when I tried it a while ago. However, it still doesn't compare with solving in complex space directly. You can maybe test it by using a matrix with a singular real part. This probably makes the ERF matrix singular and the ERF strategy may fail.

christoph-sohrmann commented 8 years ago

Hi, I actually realised that there is a quite a performance drop when using complex type instead of real. Therefore I would like to test performance of your ERF solution on my problems. Could you give me a quick hint how to use the ublas adapter together with the complex adapter?

ddemidov commented 8 years ago

Could you give me a quick hint how to use the ublas adapter together with the complex adapter?

Adapters may be combined, so something like this should work (not tested):

#include <amgcl/adapter/ublas.hpp>
#include <amgcl/adapter/complex.hpp>

...

Solver solve( amgcl::adapter::complex_matrix(amgcl::backend::map(ublas_complex_matrix)) )

It the test_complex.cpp the boost tuple adapter is combined with complex adapter in a similar way.

christoph-sohrmann commented 8 years ago

I tried it but I guess I did something wrong. Will try again. Thanks!

ddemidov commented 8 years ago

May be there is a limitation in ublas adapter. Could you post the error message or, even better, a minimal example? Another option is to get pointers to the underlying CRS arrays from ublas matrix and wrap those into boost tuple.

ddemidov commented 8 years ago

Another option is to get pointers to the underlying CRS arrays from ublas matrix and wrap those into boost tuple.

On the second thought, that is exactly what's happening inside ublas adapter, so it should not be necessary.

christoph-sohrmann commented 8 years ago

The problem was that Solver needs the matrix size as its first parameter, so I was a bit puzzled where the matrix goes. Anyway, seems to work now. However, now I am struggeling with using amgcl::adapter::complex_range on a ublas vector...

ddemidov commented 8 years ago

The problem was that Solver needs the matrix size as its first parameter

Boost tuple adapter needs that. Solver just takes the matrix.

However, now I am struggeling with using amgcl::adapter::complex_range on a ublas vector

This example does compile for me: https://gist.github.com/98ecf0d9ad9e52795868

christoph-sohrmann commented 8 years ago

Hi, here is an update on my latest activities. I have tested the ERF approach on my own data and unfortunately it appears insufficient for my problem type. Therefore I am still looking into putting std::complex directly into the solver - which seems more tricky than I expected. At least the following things need to be considered when using a complex value type:

All epsilons, residuals and norm-functions should remain real-valued.
This implies an additional value type, holding the complex inner real type.
The inner product needs to be extended for complex-valued vectors.
It is not immediately obvious which of the implemented algorithms trivially extent to the complex case. Therefore all algorithms should be disabled by default until they have been reviewed and shown to work properly on a complex type. As a very crude workaround, operations such as min() and max() could be overloaded to use std::abs().
OpenMP seems not work with complex types, especially #pragma omp parallel reduction(+:sum) complains about it.

Especially the question about which algorithms work on complex may be tricky. Emperically GMRES with Gauss-Seidel seems to work pretty good (which is known). However, there are papers discussing how AMG extends to complex:

http://ta.twi.tudelft.nl/mf/users/oosterle/oosterlee/camg.pdf

It seems there are a few things that need to be adjusted in the code (more than adding an enable_real.hpp). Do you think it is worth the effort?

ddemidov commented 8 years ago

The issue is related to #1, which should help in case of coupled systems of equations (an additional problem there is that matrix values and solution/rhs values may be of different types, like 4x4 matrices and 4x1 vectors). I made a couple of unsuccessful attempts (see nonscalar* branches) to address that, which failed mostly because I stumbled on an algorithm where I did not know how to proceed. I think I'll try again when I have some time.

Do you think it is worth the effort?

I think it is. Even if the convergence rates stay the same, we may get some speedup due to higher ratio of arithmetic operations to memory reads (which, in turn, could result in better openmp scalability). And the convergence rates should improve for many problems.

ddemidov / amgcl

Complex-valued problems #17