Commit 249d1edc authored by Oliver Sander's avatar Oliver Sander

[!320] [bugfix] complex values in FCGSolvers

Merge branch 'issue/iterative_solver_complex' into 'master'

ref:core/dune-istl

### Summary

Fix complex<double> field type in RestartedFCGSolver and CompleteFCGSolver

### Details

Passing vector type with field_type=complex<double> results in compiler errors
due to assignment of complex to double. This is fixed by changing types to
field_type instead of real_type similar to CGSolver and PCGSolver.

A test with explicit template instantiation of complex type is added to catch
this error.

See merge request [!320]

  [!320]: gitlab.dune-project.org/core/dune-istl/merge_requests/320
parents eed5aea9 6c1f8c10
Pipeline #20930 passed with stage
in 4 minutes and 50 seconds
......@@ -1507,7 +1507,7 @@ private:
virtual void apply (X& x, X& b, InverseOperatorResult& res)
{
using rAlloc = ReboundAllocatorType<X,real_type>;
using rAlloc = ReboundAllocatorType<X,field_type>;
res.clear();
Iteration iteration(*this,res);
_prec->pre(x,b); // prepare preconditioner
......@@ -1516,7 +1516,7 @@ private:
//arrays for interim values:
std::vector<X> d(_mmax+1, x); // array for directions
std::vector<X> Ad(_mmax+1, x); // array for Ad[i]
std::vector<real_type,rAlloc> ddotAd(_mmax+1,0); // array for <d[i],Ad[i]>
std::vector<field_type,rAlloc> ddotAd(_mmax+1,0); // array for <d[i],Ad[i]>
X w(x);
real_type def = _sp->norm(b); // compute norm
......@@ -1526,7 +1526,7 @@ private:
}
// some local variables
real_type alpha;
field_type alpha;
// the loop
int i=1;
......@@ -1566,7 +1566,7 @@ private:
private:
//This function is called every iteration to orthogonalize against the last search directions
virtual void orthogonalizations(const int& i_bounded,const std::vector<X>& Ad, const X& w, const std::vector<real_type,ReboundAllocatorType<X,real_type>>& ddotAd,std::vector<X>& d) {
virtual void orthogonalizations(const int& i_bounded,const std::vector<X>& Ad, const X& w, const std::vector<field_type,ReboundAllocatorType<X,field_type>>& ddotAd,std::vector<X>& d) {
// The RestartedFCGSolver uses only values with lower array index;
for (int k = 0; k < i_bounded; k++) {
d[i_bounded].axpy(-_sp->dot(Ad[k], w) / ddotAd[k], d[k]); // d[i] -= <<Ad[k],w>/<d[k],Ad[k]>>d[k]
......@@ -1574,7 +1574,7 @@ private:
}
// This function is called every mmax iterations to handle limited array sizes.
virtual void cycle(std::vector<X>& Ad,std::vector<X>& d,std::vector<real_type,ReboundAllocatorType<X,real_type> >& ddotAd,int& i_bounded) {
virtual void cycle(std::vector<X>& Ad,std::vector<X>& d,std::vector<field_type,ReboundAllocatorType<X,field_type> >& ddotAd,int& i_bounded) {
// Reset loop index and exchange the first and last arrays
i_bounded = 1;
std::swap(Ad[0], Ad[_mmax]);
......@@ -1622,7 +1622,7 @@ private:
private:
// This function is called every iteration to orthogonalize against the last search directions.
virtual void orthogonalizations(const int& i_bounded,const std::vector<X>& Ad, const X& w, const std::vector<real_type,ReboundAllocatorType<X,real_type>>& ddotAd,std::vector<X>& d) override {
virtual void orthogonalizations(const int& i_bounded,const std::vector<X>& Ad, const X& w, const std::vector<field_type,ReboundAllocatorType<X,field_type>>& ddotAd,std::vector<X>& d) override {
// This FCGSolver uses values with higher array indexes too, if existent.
for (int k = 0; k < _k_limit; k++) {
if(i_bounded!=k)
......@@ -1635,7 +1635,7 @@ private:
};
// This function is called every mmax iterations to handle limited array sizes.
virtual void cycle(std::vector<X>& Ad,std::vector<X>& d,std::vector<real_type,ReboundAllocatorType<X,real_type> >& ddotAd,int& i_bounded) override {
virtual void cycle(std::vector<X>& Ad,std::vector<X>& d,std::vector<field_type,ReboundAllocatorType<X,field_type> >& ddotAd,int& i_bounded) override {
// Only the loop index i_bounded return to 0, if it reached mmax.
i_bounded = 0;
// Now all arrays are filled and the loop in void orthogonalizations can use the whole arrays.
......
......@@ -11,24 +11,41 @@
#include <dune/istl/solvers.hh>
#include "laplacian.hh"
#include <complex>
#include <iterator>
namespace Dune
{
using Vector = BlockVector<FieldVector<double,1>>;
using Vec1 = BlockVector<FieldVector<double,1>>;
using Vec2 = BlockVector<FieldVector<std::complex<double>,1>>;
// explicit template instantiation of all iterative solvers
template class InverseOperator<Vector,Vector>;
template class LoopSolver<Vector>;
template class GradientSolver<Vector>;
template class CGSolver<Vector>;
template class BiCGSTABSolver<Vector>;
template class MINRESSolver<Vector>;
template class RestartedGMResSolver<Vector>;
template class RestartedFlexibleGMResSolver<Vector>;
template class GeneralizedPCGSolver<Vector>;
template class RestartedFCGSolver<Vector>;
template class CompleteFCGSolver<Vector>;
// field_type = double
template class InverseOperator<Vec1,Vec1>;
template class LoopSolver<Vec1>;
template class GradientSolver<Vec1>;
template class CGSolver<Vec1>;
template class BiCGSTABSolver<Vec1>;
template class MINRESSolver<Vec1>;
template class RestartedGMResSolver<Vec1>;
template class RestartedFlexibleGMResSolver<Vec1>;
template class GeneralizedPCGSolver<Vec1>;
template class RestartedFCGSolver<Vec1>;
template class CompleteFCGSolver<Vec1>;
// field_type = complex<double>
template class InverseOperator<Vec2,Vec2>;
template class LoopSolver<Vec2>;
template class GradientSolver<Vec2>;
template class CGSolver<Vec2>;
template class BiCGSTABSolver<Vec2>;
template class MINRESSolver<Vec2>;
template class RestartedGMResSolver<Vec2>;
template class RestartedFlexibleGMResSolver<Vec2>;
template class GeneralizedPCGSolver<Vec2>;
template class RestartedFCGSolver<Vec2>;
template class CompleteFCGSolver<Vec2>;
} // end namespace Dune
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment