ITL_OEMSolver.hh 3.75 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
#include "ITL_Preconditioner.h"

namespace AMDiS {
  template< typename Vector >
  struct PreconPair { 
    /// Pointer to the left preconditioner
    ITL_BasePreconditioner< Vector >* l;

    /// Pointer to the right preconditioner
    ITL_BasePreconditioner< Vector >* r;
    PreconPair():
12
13
      l(NULL), r(NULL)
    {}
14
15
16
17
18
19

  };

  struct BaseITL_runner { 
    OEMSolver& oem;
    BaseITL_runner(OEMSolver* oemPtr):
20
      oem(*oemPtr)
21
22
23
24
25
    {
      TEST_EXIT(oemPtr != NULL)("need a real oemsolver\n");
    }

    template< typename Vector , typename Matrix >
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
    void setPrecon(PreconPair< Vector >& pair, const Matrix& matrix)
    {
      /// Creator for the left preconditioner
      CreatorInterface< ITL_BasePreconditioner< Vector > >* leftCreator(NULL);

      /// Creator for the right preconditioner
      CreatorInterface< ITL_BasePreconditioner< Vector > >* rightCreator(NULL);

      std::string preconType("no");
      Parameters::get(oem.getName() + "->left precon", preconType);

      leftCreator = CreatorMap<ITL_BasePreconditioner< Vector > >::getCreator(preconType);
      TEST_EXIT(leftCreator != NULL)
	("there is no creator for the given left preconditioner");

      preconType = "no";
      Parameters::get(oem.getName() + "->right precon", preconType);
      rightCreator = CreatorMap<ITL_BasePreconditioner< Vector > >::getCreator(preconType);
      TEST_EXIT(rightCreator != NULL)
	("there is no creator for the given right preconditioner");

      pair.l = leftCreator->create(matrix);
      pair.r = rightCreator->create(matrix);
    }
50
51
  };

52
  template< typename ITLSolver, typename Matrix, typename Vector >
53
54
  struct ITL_OEMSolver_para_runner : BaseITL_runner {
    ITL_OEMSolver_para_runner(OEMSolver* oem_):
55
56
      BaseITL_runner(oem_),
      ell(1)
57
58
59
60
61
    {
      ell=oem.getMaxIterations();
      Parameters::get(oem.getName() + "->ell", ell);
    }

62
63
64
65
66
    void init(const Matrix& A)
    {
	BaseITL_runner::setPrecon(preconPair, A); 
    }

67
68
    int solve(const Matrix& A, Vector& x, Vector& b) 
    {
69
70
71
      FUNCNAME("ITL_OEMSolver_para_runner::solve()");
      TEST_EXIT(preconPair.l != NULL)("there is no left preconditioner");
      TEST_EXIT(preconPair.r != NULL)("there is no right preconditioner");
72
73
74

      mtl::dense_vector<typename Matrix::value_type> r(A * x - b); 
      itl::cyclic_iteration<typename Vector::value_type> iter(r, oem.getMaxIterations(), oem.getRelative(), 
75
							      oem.getTolerance(), oem.getPrint_cycle());
76

77
      int error = solver(A, x, b, *(preconPair.l), *(preconPair.r), iter, ell);
78
79
80
      oem.setIterations(iter.iterations());
      oem.setResidual(iter.resid());
      oem.setErrorCode(error);
81

82
83
      return error;
    }
84
  private:
85
86
    PreconPair< Vector > preconPair;
    ITLSolver solver;
87
88
89
    int ell;
  };

90
  template< typename ITLSolver, typename Matrix, typename Vector >
91
92
93
94
  struct ITL_OEMSolver_runner : BaseITL_runner {
   
   
    ITL_OEMSolver_runner(OEMSolver* oem):
95
      BaseITL_runner(oem)
96
97
    {}

98
99
100
101
102
    void init(const Matrix& A)
    {
	BaseITL_runner::setPrecon(preconPair, A);
    }
    int solve(const Matrix& A , Vector& x, Vector& b) 
103
    {
104
105
106
      FUNCNAME("ITL_OEMSolver_runner::solve()");
      TEST_EXIT(preconPair.l != NULL)("there is no left preconditioner");
      TEST_EXIT(preconPair.r != NULL)("there is no right preconditioner");
107
108
109

      mtl::dense_vector<typename Matrix::value_type> r(A * x - b); 
      itl::cyclic_iteration<typename Matrix::value_type> iter(r, oem.getMaxIterations(), oem.getRelative(), 
110
							      oem.getTolerance(), oem.getPrint_cycle());
111

112
      int error = solver(A, x, b, *(preconPair.l), *(preconPair.r), iter);
113
114
115
      oem.setErrorCode(error);
      oem.setIterations(iter.iterations());
      oem.setResidual(iter.resid());
116

117
118
      return error;
    }
119
120
121
    private:
    ITLSolver solver;
    PreconPair< Vector > preconPair;
122
123
  };
}