ug4
Loading...
Searching...
No Matches
external_solvers.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2014-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__EXTERNAL_SOLVER_
34#define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__EXTERNAL_SOLVER_
35
36#include "common/common.h"
38
39
40#ifdef UG_PARALLEL
42#endif
43#include "common/progress.h"
45
48
49
50namespace ug{
51
53{
54public:
55 virtual bool init(const CPUAlgebra::matrix_type &A) = 0;
57 virtual const char* name() const = 0;
59};
60
61template <typename TAlgebra>
63 : public IMatrixOperatorInverse<typename TAlgebra::matrix_type,
64 typename TAlgebra::vector_type>,
65 public VectorDebugWritingObject<typename TAlgebra::vector_type>
66{
67 public:
68
69 virtual const char *double_name() const = 0;
70
71 virtual const char* name() const
72 {
73 return double_name();
74 }
75
76 // Algebra type
77 typedef TAlgebra algebra_type;
78
79 // Vector type
80 typedef typename TAlgebra::vector_type vector_type;
81
82 // Matrix type
83 typedef typename TAlgebra::matrix_type matrix_type;
84
86
87 public:
88 // Constructor
95
96 // Clone
97
99 {
100 UG_THROW("");
101 return SPNULL;
102 }
103
104
106 virtual bool supports_parallel() const {return false;}
107
109 void set_disable_preprocessing(bool bDisable) {m_bDisablePreprocessing = bDisable;}
110
112
113 virtual void double_init(const CPUAlgebra::matrix_type &mat) = 0;
114
116 {
117 // do not do a thing if preprocessing disabled
119 return;
120
121 if( A.num_rows() == 0 || A.num_cols() == 0) { m_size = 0; return; }
122 STATIC_ASSERT(matrix_type::rows_sorted, Matrix_has_to_have_sorted_rows);
123
125
126 #ifdef UG_PARALLEL
127 // add slave rows to master rows (in indices which this is possible for)
128 matrix_type A_tmp; A_tmp = A;
131 }else {
133
134 // set zero on slaves
135 std::vector<IndexLayout::Element> vIndex;
136 CollectUniqueElements(vIndex, A.layouts()->slave());
137 SetDirichletRow(A_tmp, vIndex);
138 }
139 // Even after this setting of Dirichlet rows, it is possible that there are
140 // zero rows on a proc because of the distribution:
141 // For example, if one has a horizontal grid interface between two SHADOW_RIM_COPY
142 // vertices, but the shadowing element for the hSlave side is vMaster (without being in
143 // any horizontal interface). Then, as the horizontal interface (on the shadowed level)
144 // is not part of the algebraic layouts, the hSlave is not converted into a Dirichlet row
145 // by the previous commands.
146 // As a first aid, we will simply convert any zero row on the current proc into a
147 // Dirichlet row.
148 // TODO: The corresponding rhs vector entry could be non-zero!
149 // It is definitely not set to zero by change_storage_type(PST_UNIQUE) as the index is not contained
150 // in the vector layouts either. Still, the defect assembling process might contain a vertex
151 // loop and assemble something that is not solution-dependent! What do we do then?
152 size_t nInd = A_tmp.num_rows();
153 for (size_t i = 0; i < nInd; ++i)
154 {
155 if (!A_tmp.num_connections(i))
156 A_tmp(i,i) = 1.0;
157 }
158
161 #else
162 const matrix_type& A_tmp = A;
163 #endif
164
166 m_blockSize = mat.num_rows()/A_tmp.num_rows();
167
168 if(m_size != 0)
169 double_init(mat);
170 }
171
174 {
175 // remember operator
176 m_spOperator = Op;
177
178 // init LU operator
180 // we're done
181 return true;
182 }
183
184
185 protected:
186
187 // Preprocess routine
189 {
190 matrix_type &A = *pOp;
192
193 return true;
194 }
195
197 {
198 for(size_t i=0, k=0; i<v.size(); i++)
199 {
200 for(size_t j=0; j<GetSize(v[i]); j++)
201 v_scalar[k++] = BlockRef(v[i],j);
202 }
203 }
204
206 {
207 for(size_t i=0, k=0; i<v.size(); i++)
208 {
209 for(size_t j=0; j<GetSize(v[i]); j++)
210 BlockRef(v[i],j) = v_scalar[k++];
211 }
212 }
213
215
216
217 // Stepping routine
218 virtual bool apply(vector_type& c, const vector_type& d)
219 {
220 if (m_size == 0)
221 {
222#ifdef UG_PARALLEL
223 c.set_storage_type(PST_CONSISTENT);
224#endif
225 return true;
226 }
229
230#ifdef UG_PARALLEL
231
234 // make defect unique
235 SmartPtr<vector_type> spDtmp = d.clone();
237 spDtmp->change_storage_type(PST_CONSISTENT);
238 } else{
239 spDtmp->change_storage_type(PST_UNIQUE);
240 }
241
242
243 get_vector(m_d, *spDtmp);
244#else
245 get_vector(m_d, d);
246#endif
247
248 m_c.set(0.0);
249
250
252
253 set_vector(m_c, c);
254
255#ifdef UG_PARALLEL
256 // correction must always be consistent (but is unique by construction)
257 // (or regarded as such in the case of consistent interfaces)
258 c.set_storage_type(PST_UNIQUE);
259 c.change_storage_type(PST_CONSISTENT);
260#endif
261 return true;
262 }
263
265 {
266 // solve u
267 if(!apply(u, f)) return false;
268
269 // update defect
270 if(!m_spOperator->matmul_minus(f, u))
271 {
272 UG_LOG("ERROR in 'LU::apply_return_defect': "
273 "Cannot apply matmul_minus.\n");
274 return false;
275 }
276
277 // we're done
278 return true;
279 }
280
281public:
282 using VectorDebugWritingObject<typename TAlgebra::vector_type>::vector_debug_writer;
283
285 {
286 if(vector_debug_writer().valid() == false) return -1;
287 vector_debug_writer()->update_positions();
288 return vector_debug_writer()->get_dim();
289 }
290
291 template<int Tdim>
292 bool get_positions(std::vector<MathVector<Tdim> > &coord)
293 {
294 UG_COND_THROW(vector_debug_writer().valid() == false, "no debug writer set.");
295 int dim = get_dim();
296 UG_COND_THROW(dim != Tdim, "wrong dimension");
297
298 return copy_pos<Tdim, Tdim>(coord, get_positions<Tdim>());
299 }
300
301 bool get_positions3(std::vector<MathVector<3> > &coord)
302 {
303 UG_COND_THROW(vector_debug_writer().valid() == false, "no debug writer set.");
304 int dim = get_dim();
305 switch(dim)
306 {
307 case 1:
308 return copy_pos(coord, get_positions<1>());
309 case 2:
310 return copy_pos(coord, get_positions<2>());
311 case 3:
312 return copy_pos(coord, get_positions<3>());
313 case -1:
314 return false;
315 }
316 }
317
318
319 template<int dim>
320 const std::vector<MathVector<dim> > &get_positions()
321 {
322 if(vector_debug_writer().valid())
323 {
324 vector_debug_writer()->update_positions();
325 return vector_debug_writer()->template get_positions<dim>();
326 }
327 else UG_THROW("no debug_writer!");
328 }
329 template<int dim1, int dim2>
330 bool copy_pos(std::vector<MathVector<dim1> > &dest, const std::vector<MathVector<dim2> > &src)
331 {
332 UG_COND_THROW(m_size == 0 || m_blockSize == 0, "not initialized");
333 UG_COND_THROW(dim1 < dim2, "loss of data");
334
335 dest.resize(m_size);
336 for(size_t i=0; i<src.size(); i++)
337 {
338 for(size_t k=0; k<m_blockSize; k++)
339 {
340 dest[i*m_blockSize + k]=0;
341 for(size_t j=0; j<dim2; j++)
342 dest[i*m_blockSize + k][j] = src[i][j];
343 }
344 }
345 return true;
346 }
347
348
349 protected:
350 // Postprocess routine
351 virtual bool postprocess() {return true;}
352
355 size_t m_size;
358
359
360};
361
362} // namespace ug
363
364#endif /* __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__EXTERNAL_SOLVER_ */
Definition smart_pointer.h:108
Definition external_solvers.h:66
CPUAlgebra::vector_type m_d
Definition external_solvers.h:354
CPUAlgebra::vector_type m_c
Definition external_solvers.h:354
virtual void double_init(const CPUAlgebra::matrix_type &mat)=0
virtual bool double_apply(CPUAlgebra::vector_type &c, const CPUAlgebra::vector_type &d)=0
void mat_preprocess(const matrix_type &A)
Definition external_solvers.h:115
virtual bool apply(vector_type &c, const vector_type &d)
applies the inverse operator, i.e. returns u = A^{-1} * f
Definition external_solvers.h:218
int get_dim()
Definition external_solvers.h:284
const std::vector< MathVector< dim > > & get_positions()
Definition external_solvers.h:320
virtual const char * double_name() const =0
bool get_positions(std::vector< MathVector< Tdim > > &coord)
Definition external_solvers.h:292
SmartPtr< ILinearIterator< vector_type > > clone()
clone
Definition external_solvers.h:98
SmartPtr< MatrixOperator< matrix_type, vector_type > > m_spOperator
Definition external_solvers.h:172
size_t m_blockSize
Definition external_solvers.h:356
virtual bool supports_parallel() const
returns if parallel solving is supported
Definition external_solvers.h:106
virtual const char * name() const
returns the name of the operator inverse
Definition external_solvers.h:71
bool copy_pos(std::vector< MathVector< dim1 > > &dest, const std::vector< MathVector< dim2 > > &src)
Definition external_solvers.h:330
TAlgebra::vector_type vector_type
Definition external_solvers.h:80
TAlgebra::matrix_type matrix_type
Definition external_solvers.h:83
bool get_positions3(std::vector< MathVector< 3 > > &coord)
Definition external_solvers.h:301
void get_vector(CPUAlgebra::vector_type &v_scalar, const vector_type &v)
Definition external_solvers.h:196
void enable_consistent_interfaces(bool enable)
Definition external_solvers.h:111
virtual bool postprocess()
Definition external_solvers.h:351
bool m_bUseConsistentInterfaces
Definition external_solvers.h:353
bool m_bDisablePreprocessing
Definition external_solvers.h:357
virtual bool preprocess(SmartPtr< MatrixOperator< matrix_type, vector_type > > pOp)
Definition external_solvers.h:188
void set_disable_preprocessing(bool bDisable)
disable preprocessing (if underlying matrix has not changed)
Definition external_solvers.h:109
virtual bool init(SmartPtr< MatrixOperator< matrix_type, vector_type > > Op)
Definition external_solvers.h:173
TAlgebra algebra_type
Definition external_solvers.h:77
virtual bool apply_return_defect(vector_type &u, vector_type &f)
applies the inverse operator and updates the defect
Definition external_solvers.h:264
size_t m_size
Definition external_solvers.h:355
void set_vector(CPUAlgebra::vector_type &v_scalar, vector_type &v)
Definition external_solvers.h:205
IExternalSolver()
Definition external_solvers.h:89
Definition external_solvers.h:53
virtual const char * name() const =0
virtual bool apply(CPUAlgebra::vector_type &c, const CPUAlgebra::vector_type &d)=0
virtual bool init(const CPUAlgebra::matrix_type &A)=0
virtual ~IExternalSolverImplementation()
Definition external_solvers.h:58
describes an inverse linear mapping X->Y based on a matrix
Definition matrix_operator_inverse.h:64
a mathematical Vector with N entries.
Definition math_vector.h:97
Definition matrix_operator.h:49
virtual M & get_matrix()
Definition matrix_operator.h:74
Wrapper for sequential matrices to handle them in parallel.
Definition parallel_matrix.h:65
void set_storage_type(uint type)
sets the storage type
Definition parallel_matrix.h:101
void set_layouts(ConstSmartPtr< AlgebraLayouts > layouts)
sets the algebra layouts
Definition parallel_matrix.h:97
void set(number w, ParallelStorageType type)
set all entries to value and the storage type
Definition parallel_vector_impl.h:223
void set_storage_type(uint type)
sets the storage type
Definition parallel_vector.h:104
Definition debug_writer.h:266
SmartPtr< IVectorDebugWriter< vector_type > > vector_debug_writer()
returns the debug writer
Definition debug_writer.h:286
void resize(size_t newSize, bool bCopyValues=true)
Definition vector.h:109
void MatAddSlaveRowsToMasterRowOverlap0(TMatrix &mat)
Generates a set of unique global algebra ids.
Definition parallelization_util.h:106
void MatMakeConsistentOverlap0(TMatrix &mat)
Generates a set of unique global algebra ids.
Definition parallelization_util.h:126
@ PST_CONSISTENT
Definition parallel_storage_type.h:68
@ PST_UNIQUE
Definition parallel_storage_type.h:70
@ PST_ADDITIVE
Definition parallel_storage_type.h:69
size_t GetDoubleSparseFromBlockSparse(TDoubleSparse &A, const TSparseMatrix &S)
Definition sparsematrix_util.h:1217
void SetDirichletRow(TSparseMatrix &A, size_t i, size_t alpha)
Definition sparsematrix_util.h:796
virtual void init()
const NullSmartPtr SPNULL
The equivalent to NULL for smart pointers.
Definition smart_pointer.h:90
#define UG_THROW(msg)
Definition error.h:57
#define STATIC_ASSERT(expr, msg)
Definition common.h:61
#define UG_LOG(msg)
Definition log.h:367
#define UG_COND_THROW(cond, msg)
UG_COND_THROW(cond, msg) : performs a UG_THROW(msg) if cond == true.
Definition error.h:61
the ug namespace
double & BlockRef(T &vec, size_t i)
Definition blocks.h:66
SmartPtr< AlgebraLayouts > CreateLocalAlgebraLayouts()
Definition parallelization_util.cpp:936
size_t GetSize(const T &t)