MTL4Solver.h 4.06 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  TU Dresden                                                            ==
// ==                                                                        ==
// ==  Institut für Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
// ==                                                                        ==
// ============================================================================

#ifndef MTL4SOLVER_H
#define MTL4SOLVER_H

#include "OEMSolver.h"
#include "MatrixStreams.h"
#include <iostream>
Naumann, Andreas's avatar
Naumann, Andreas committed
26
27
28
29
30
31
#include <boost/mpl/bool.hpp>
#include <boost/numeric/mtl/utility/is_distributed.hpp>

#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) && defined(HAVE_PARALLEL_MTL4)
# include <boost/numeric/mtl/par/distribution.hpp>
#endif
32
33
34

namespace AMDiS {

35
36
  template< typename MTLMatrix, typename MTLVector, typename Worker >
  class MTL4Solver : public OEMSolver {    
37
    MTLMatrix matrix;
38
39
  protected:    
    Worker worker;
40

Naumann, Andreas's avatar
Naumann, Andreas committed
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
    template< typename Mapper, bool IsDist>
    void init(Mapper& mapper, boost::mpl::bool_<IsDist>)
    {
        matrix.change_dim(mapper.getNumRows(), mapper.getNumCols());
    }

#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) && defined(HAVE_PARALLEL_MTL4)
    void init(ParallelMapper& mapper, boost::mpl::true_)
    {
	mtl::par::block_distribution dist(mapper.getNumRows() /* , communicator */ );
	dist.setup_from_local_size(mapper.meshDistributor().getNumberRankDofs()*mapper.getNumComponents());
        matrix.init_distribution(dist, dist, mapper.getNumRows(), mapper.getNumRows());
    }
#endif

    template< bool IsDist>
    void init(MTLVector& v, boost::mpl::bool_<IsDist>)
    {
      v.change_dim(num_rows(matrix));
    }

    void init(MTLVector& v, boost::mpl::true_)
    {
      v.init_distribution(row_distribution(matrix), num_rows(matrix));
    }

67
68
69
    template< typename Matrix, typename Vector, typename Mapper >
    int solve(const Matrix& A, Vector& x, Vector& b, Mapper& mapper) 
    {
70
      if(num_rows(matrix) == 0 || !getMultipleRhs() ) {
Naumann, Andreas's avatar
Naumann, Andreas committed
71
72
	init(mapper, mtl::traits::is_distributed<MTLMatrix>());

73
74
75
76
77
        set_to_zero(matrix);
        MatMap< const Matrix, Mapper > matMap(A,mapper);
        matrix << matMap;
	worker.init(matrix);
      }
78

Naumann, Andreas's avatar
Naumann, Andreas committed
79
80
      MTLVector mtl_x;
      init(mtl_x,mtl::traits::is_distributed<MTLVector>());
81
82
83
      set_to_zero(mtl_x);
      VecMap< Vector, Mapper > xVecMap(x, mapper);
      mtl_x << xVecMap;
84

Naumann, Andreas's avatar
Naumann, Andreas committed
85
86
      MTLVector mtl_b;
      init(mtl_b,mtl::traits::is_distributed<MTLVector>());
87
88
89
      set_to_zero(mtl_b);
      VecMap< Vector, Mapper> bVecMap(b, mapper);
      mtl_b << bVecMap;
90

91
      error = worker.solve(matrix ,mtl_x, mtl_b);
92

Naumann, Andreas's avatar
Naumann, Andreas committed
93
      MTLVector r(mtl_b); 
94
95
96
97
      r -= matrix * mtl_x; 
      double residual = two_norm(r);
      MSG("MTL4Solver: ||b-Ax||= %e\n", residual);

98
99
100
      mtl_x >> xVecMap;
      return error;
    }
101

102
103
104
  public:
    MTL4Solver(std::string n):
      OEMSolver(n),
105
      matrix(0,0),
106
107
      worker(this)
    {}
108

109
    virtual int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
110
111
112
			    SystemVector& x, 
			    SystemVector& b,
			    VectorialMapper& m)
113
114
115
    {
      return solve(A,x,b,m);
    }
116
117
118
119
120

      
  };
}
#endif