MTL4Solver.h 4.25 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
// ============================================================================
// ==                                                                        ==
// == 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"
Naumann, Andreas's avatar
Naumann, Andreas committed
25
#include "Timer.h"
26
#include <iostream>
Naumann, Andreas's avatar
Naumann, Andreas committed
27
28
29
30
31
32
#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
33
34
35

namespace AMDiS {

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

Naumann, Andreas's avatar
Naumann, Andreas committed
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
67
    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));
    }

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

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

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

Naumann, Andreas's avatar
Naumann, Andreas committed
87
88
      MTLVector mtl_b;
      init(mtl_b,mtl::traits::is_distributed<MTLVector>());
89
90
91
      set_to_zero(mtl_b);
      VecMap< Vector, Mapper> bVecMap(b, mapper);
      mtl_b << bVecMap;
Naumann, Andreas's avatar
Naumann, Andreas committed
92
93
   
      MSG("fill MTL4 matrix needed %5.f seconds\n", t.elapsed());
94

Naumann, Andreas's avatar
Naumann, Andreas committed
95
      t.reset();
96
      error = worker.solve(matrix ,mtl_x, mtl_b);
Naumann, Andreas's avatar
Naumann, Andreas committed
97
      MSG("solve MTL4 matrix needed %5.f seconds\n", t.elapsed());
98

Naumann, Andreas's avatar
Naumann, Andreas committed
99
      MTLVector r(mtl_b); 
100
101
102
103
      r -= matrix * mtl_x; 
      double residual = two_norm(r);
      MSG("MTL4Solver: ||b-Ax||= %e\n", residual);

104
105
106
      mtl_x >> xVecMap;
      return error;
    }
107

108
109
110
  public:
    MTL4Solver(std::string n):
      OEMSolver(n),
111
      matrix(0,0),
112
113
      worker(this)
    {}
114

115
    virtual int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
116
117
118
			    SystemVector& x, 
			    SystemVector& b,
			    VectorialMapper& m)
119
120
121
    {
      return solve(A,x,b,m);
    }
122
123
124
125
126

      
  };
}
#endif