MTL4Solver.h 4.43 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
#include <boost/mpl/bool.hpp>
#include <boost/numeric/mtl/utility/is_distributed.hpp>

#if defined(HAVE_PARALLEL_DOMAIN_AMDIS) && defined(HAVE_PARALLEL_MTL4)
31
#include <boost/numeric/mtl/par/distribution.hpp>
Naumann, Andreas's avatar
Naumann, Andreas committed
32
#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
    template< typename Mapper, bool IsDist>
    void init(Mapper& mapper, boost::mpl::bool_<IsDist>)
    {
45
      matrix.change_dim(mapper.getNumRows(), mapper.getNumCols());
Naumann, Andreas's avatar
Naumann, Andreas committed
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
    }

#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
    template< typename Matrix, typename Vector, typename Mapper >
69
70
71
    int solve(const Matrix& A, Vector& x, Vector& b, Mapper& mapper,
	      bool createMatrixData,
	      bool storeMatrixData) 
72
    {
73
74
      FUNCNAME("MTL4Solver::solve()");

Naumann, Andreas's avatar
Naumann, Andreas committed
75
      Timer t;
76
      if (createMatrixData) {
Naumann, Andreas's avatar
Naumann, Andreas committed
77
78
	init(mapper, mtl::traits::is_distributed<MTLMatrix>());

79
        set_to_zero(matrix);
80
        MatMap< const Matrix, Mapper > matMap(A, mapper);
81
82
83
        matrix << matMap;
	worker.init(matrix);
      }
84

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

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

Naumann, Andreas's avatar
Naumann, Andreas committed
99
      t.reset();
100
      error = worker.solve(matrix ,mtl_x, mtl_b);
Naumann, Andreas's avatar
Naumann, Andreas committed
101
      MSG("solve MTL4 matrix needed %.5f seconds\n", t.elapsed());
102

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

108
      mtl_x >> xVecMap;
109
110
111

      worker.exit();

112
113
      return error;
    }
114

115
116
117
  public:
    MTL4Solver(std::string n):
      OEMSolver(n),
118
      matrix(0,0),
119
120
      worker(this)
    {}
121

122
    virtual int solveSystem(const SolverMatrix<Matrix<DOFMatrix*> >& A,
123
124
			    SystemVector& x, 
			    SystemVector& b,
125
126
127
			    VectorialMapper& m,
			    bool createMatrixData,
			    bool storeMatrixData) 
128
    {
129
      return solve(A, x, b, m, createMatrixData, storeMatrixData);
130
    }
131
132
133
134
135

      
  };
}
#endif