Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Martin Rupp
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__ILUT_SCALAR__
34 : #define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__ILUT_SCALAR__
35 :
36 : #include "common/util/smart_pointer.h"
37 : #include "lib_algebra/operator/interface/preconditioner.h"
38 : #ifdef UG_PARALLEL
39 : #include "lib_algebra/parallelization/parallelization.h"
40 : #endif
41 : #include "common/progress.h"
42 : #include "common/util/ostream_util.h"
43 :
44 : #include "lib_algebra/algebra_common/vector_util.h"
45 : #include "lib_algebra/operator/linear_solver/linear_solver.h"
46 : #include "lib_algebra/cpu_algebra_types.h"
47 : #include "ilut.h"
48 :
49 : namespace ug{
50 :
51 : template <typename TAlgebra>
52 : class ILUTScalarPreconditioner : public IPreconditioner<TAlgebra>
53 : {
54 : public:
55 : // Algebra type
56 : typedef TAlgebra algebra_type;
57 :
58 : // Vector type
59 : typedef typename TAlgebra::vector_type vector_type;
60 :
61 : // Matrix type
62 : typedef typename TAlgebra::matrix_type matrix_type;
63 :
64 : /// Matrix Operator type
65 : typedef typename IPreconditioner<TAlgebra>::matrix_operator_type matrix_operator_type;
66 :
67 : private:
68 : typedef typename matrix_type::value_type block_type;
69 : typedef IPreconditioner<TAlgebra> base_type;
70 :
71 : public:
72 : // Constructor
73 0 : ILUTScalarPreconditioner(double eps=1e-6) :
74 0 : m_eps(eps), m_info(false), m_show_progress(true), m_bSort(true)
75 0 : {};
76 :
77 : /// clone constructor
78 0 : ILUTScalarPreconditioner( const ILUTScalarPreconditioner<TAlgebra> &parent )
79 0 : : base_type(parent)
80 : {
81 0 : m_eps = parent.m_eps;
82 0 : set_info(parent.m_info);
83 0 : set_show_progress(parent.m_show_progress);
84 0 : set_sort(parent.m_bSort);
85 0 : }
86 :
87 : /// Clone
88 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
89 : {
90 0 : return make_sp(new ILUTScalarPreconditioner<algebra_type>(*this));
91 : }
92 :
93 : // Destructor
94 0 : virtual ~ILUTScalarPreconditioner()
95 : {
96 0 : };
97 :
98 : /// returns if parallel solving is supported
99 0 : virtual bool supports_parallel() const {return true;}
100 :
101 : /// sets threshold for incomplete LU factorisation (added 01122010ih)
102 0 : void set_threshold(number thresh)
103 : {
104 0 : m_eps = thresh;
105 0 : }
106 :
107 : /// sets storage information output to true or false
108 0 : void set_info(bool info)
109 : {
110 0 : m_info = info;
111 0 : }
112 :
113 : /// sets the indication of the progress to true or false
114 : void set_show_progress(bool b)
115 : {
116 0 : m_show_progress = b;
117 : }
118 :
119 0 : void set_sort(bool b)
120 : {
121 0 : m_bSort = b;
122 0 : }
123 :
124 : public:
125 0 : void preprocess(const matrix_type &M)
126 : {
127 : #ifdef UG_PARALLEL
128 : matrix_type A;
129 : A = M;
130 :
131 : MatAddSlaveRowsToMasterRowOverlap0(A);
132 :
133 : // set zero on slaves
134 : std::vector<IndexLayout::Element> vIndex;
135 : CollectUniqueElements(vIndex, M.layouts()->slave());
136 : SetDirichletRow(A, vIndex);
137 : #else
138 : const matrix_type &A = M;
139 : #endif
140 :
141 : STATIC_ASSERT(matrix_type::rows_sorted, Matrix_has_to_have_sorted_rows);
142 :
143 0 : ilut = make_sp(new ILUTPreconditioner<CPUAlgebra>(m_eps));
144 0 : ilut->set_info(m_info);
145 0 : ilut->set_show_progress(m_show_progress);
146 0 : ilut->set_sort(m_bSort);
147 :
148 0 : mo = make_sp(new MatrixOperator<CPUAlgebra::matrix_type, CPUAlgebra::vector_type>);
149 0 : CPUAlgebra::matrix_type &mat = mo->get_matrix();
150 :
151 : #ifdef UG_PARALLEL
152 : mat.set_storage_type(PST_ADDITIVE);
153 : mat.set_layouts(CreateLocalAlgebraLayouts());
154 : #endif
155 :
156 0 : m_size = GetDoubleSparseFromBlockSparse(mat, A);
157 :
158 0 : SmartPtr<StdConvCheck<CPUAlgebra::vector_type> > convCheck(new StdConvCheck<CPUAlgebra::vector_type>);
159 : convCheck->set_maximum_steps(100);
160 : convCheck->set_minimum_defect(1e-50);
161 : convCheck->set_reduction(1e-20);
162 : convCheck->set_verbose(false);
163 :
164 0 : linearSolver.set_preconditioner(ilut);
165 0 : linearSolver.set_convergence_check(convCheck);
166 0 : linearSolver.init(mo);
167 0 : }
168 : protected:
169 : // Name of preconditioner
170 0 : virtual const char* name() const {return "ILUTScalar";}
171 :
172 :
173 : // Preprocess routine
174 0 : virtual bool preprocess(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp)
175 : {
176 0 : matrix_type &A = *pOp;
177 0 : preprocess(A);
178 :
179 0 : return true;
180 : }
181 :
182 : void get_vector(CPUAlgebra::vector_type &v_scalar, const vector_type& v)
183 : {
184 0 : for(size_t i=0, k=0; i<v.size(); i++)
185 : {
186 0 : for(size_t j=0; j<GetSize(v[i]); j++)
187 0 : v_scalar[k++] = BlockRef(v[i],j);
188 : }
189 : }
190 :
191 : void set_vector(CPUAlgebra::vector_type &v_scalar, vector_type& v)
192 : {
193 0 : for(size_t i=0, k=0; i<v.size(); i++)
194 : {
195 0 : for(size_t j=0; j<GetSize(v[i]); j++)
196 0 : BlockRef(v[i],j) = v_scalar[k++];
197 : }
198 : }
199 :
200 : // Stepping routine
201 0 : virtual bool step(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp, vector_type& c, const vector_type& d)
202 : {
203 : #ifdef UG_PARALLEL
204 : SmartPtr<vector_type> spDtmp = d.clone();
205 : spDtmp->change_storage_type(PST_UNIQUE);
206 : bool b = apply_double(c, *spDtmp);
207 :
208 : c.set_storage_type(PST_ADDITIVE);
209 : c.change_storage_type(PST_CONSISTENT);
210 : return b;
211 : #else
212 0 : return apply_double(c, d);
213 : #endif
214 : return true;
215 : }
216 :
217 0 : bool apply_double(vector_type &c, const vector_type &d)
218 : {
219 0 : m_c.resize(m_size);
220 0 : m_d.resize(m_size);
221 :
222 : #ifdef UG_PARALLEL
223 : m_d.set_storage_type(PST_ADDITIVE);
224 : m_c.set_storage_type(PST_CONSISTENT);
225 : m_c.set_layouts(CreateLocalAlgebraLayouts());
226 : m_d.set_layouts(CreateLocalAlgebraLayouts());
227 : #endif
228 : get_vector(m_d, d);
229 : m_c.set(0.0);
230 :
231 : //ApplyLinearSolver(mo, m_u, m_b, linearSolver);
232 0 : ilut->apply(m_c, m_d);
233 : //ilut->apply(m_c, m_d);
234 :
235 : set_vector(m_c, c);
236 0 : return true;
237 : }
238 :
239 :
240 :
241 : public:
242 0 : bool solve(vector_type &c, const vector_type &d)
243 : {
244 0 : m_c.resize(m_size);
245 0 : m_d.resize(m_size);
246 :
247 : #ifdef UG_PARALLEL
248 : m_d.set_storage_type(PST_ADDITIVE);
249 : m_c.set_storage_type(PST_CONSISTENT);
250 : m_c.set_layouts(CreateLocalAlgebraLayouts());
251 : m_d.set_layouts(CreateLocalAlgebraLayouts());
252 : #endif
253 : get_vector(m_d, d);
254 : m_c.set(0.0);
255 :
256 : //ApplyLinearSolver(mo, m_u, m_b, linearSolver);
257 0 : linearSolver.apply_return_defect(m_c, m_d);
258 : //ilut->apply(m_c, m_d);
259 :
260 : set_vector(m_c, c);
261 0 : return true;
262 : }
263 :
264 : public:
265 0 : virtual std::string config_string() const
266 : {
267 0 : std::stringstream ss ; ss << "ILUTScalar(threshold = " << m_eps << ", sort = " << (m_bSort?"true":"false") << ")";
268 0 : if(m_eps == 0.0) ss << " = Sparse LU";
269 0 : return ss.str();
270 0 : }
271 :
272 : protected:
273 : // Postprocess routine
274 0 : virtual bool postprocess() {return true;}
275 :
276 : protected:
277 : SmartPtr<ILUTPreconditioner<CPUAlgebra> > ilut;
278 : SmartPtr<MatrixOperator<CPUAlgebra::matrix_type, CPUAlgebra::vector_type> > mo;
279 : LinearSolver<CPUAlgebra::vector_type> linearSolver;
280 : CPUAlgebra::vector_type m_c, m_d;
281 : double m_eps;
282 : bool m_info, m_show_progress, m_bSort;
283 : size_t m_size;
284 : };
285 :
286 : } // end namespace ug
287 :
288 : #endif
|
