Plugins
Loading...
Searching...
No Matches
pcr_ilut.h
Go to the documentation of this file.
1/*
2 * Copyright (c) 2013-2014: G-CSC, Goethe University Frankfurt
3 * Authors: Martin Rupp, Christian Wehner
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
34#ifndef __H__UG__PLUGINS__NAVIER_STOKES__INCOMPRESSIBLE__FVCR__PCRILUT__
35#define __H__UG__PLUGINS__NAVIER_STOKES__INCOMPRESSIBLE__FVCR__PCRILUT__
36
39#ifdef UG_PARALLEL
41#endif
42
43namespace ug{
44
46/* Threshold ILU method
47*
48* parallelization equivalent to ILU parallelization (see ilu.h)
49*
50* The matrix entries can be devided into four sectors:
51* velocity-velocity, velocity-pressure, pressure-velocity, pressure-pressure
52* Threshold tolerance depends on sector of new entry
53* (given by parameters m_eps_vv, m_eps_vp, m_eps_pv, m_eps_pp)
54* Type of entry a_ij is determined by diagonal entry of rows i and j
55* if a_ii=0 i is a pressure index, else a velocity index
56*
57*/
58template <typename TAlgebra>
59class PCRILUTPreconditioner : public IPreconditioner<TAlgebra>
60{
61 public:
62 // Algebra type
63 typedef TAlgebra algebra_type;
64
65 // Vector type
66 typedef typename TAlgebra::vector_type vector_type;
67
68 // Matrix type
69 typedef typename TAlgebra::matrix_type matrix_type;
70
73
74 private:
75 typedef typename matrix_type::value_type block_type;
76 using IPreconditioner<TAlgebra>::debug_writer;
77 using IPreconditioner<TAlgebra>::set_debug;
78
79 public:
80 // Constructor
81 PCRILUTPreconditioner(double eps=1e-6) :
82 m_eps(eps), m_eps_vv(eps), m_eps_vp(eps), m_eps_pv(eps), m_eps_pp(eps)
83 {
84 m_info = false;
85 };
86
87 PCRILUTPreconditioner(number threshvv,number thresh_vp_pv_pp)
88 {
89 set_threshold(threshvv,thresh_vp_pv_pp);
90 m_info = false;
91 }
92
93 PCRILUTPreconditioner(number threshvv,number threshvp,number threshpv,number threshpp)
94 {
95 set_threshold(threshvv,threshvp,threshpv,threshpp);
96 m_info = false;
97 };
98
99 PCRILUTPreconditioner(double eps,bool info) :
100 m_eps(eps), m_eps_vv(eps), m_eps_vp(eps), m_eps_pv(eps), m_eps_pp(eps), m_info(info)
101 {
102 }
103
104 PCRILUTPreconditioner(number threshvv,number thresh_vp_pv_pp,bool info)
105 {
106 set_threshold(threshvv,thresh_vp_pv_pp);
107 m_info = info;
108 }
109
110 PCRILUTPreconditioner(number threshvv,number threshvp,number threshpv,number threshpp,bool info)
111 {
112 set_threshold(threshvv,threshvp,threshpv,threshpp);
113 m_info = info;
114 }
115
116 // Clone
117
119 {
121 newInst->set_debug(debug_writer());
122 newInst->set_damp(this->damping());
123 newInst->set_threshold(this->m_eps_vv,this->m_eps_vp,this->m_eps_pv,this->m_eps_pp);
124 newInst->set_info(this->m_info);
125 return newInst;
126 }
127
128 // Destructor
130 {
131 };
132
134 virtual bool supports_parallel() const {return true;}
135
137 void set_threshold(number threshvv,number threshvp,number threshpv,number threshpp)
138 {
139 m_eps_vv = threshvv;
140 m_eps_vp = threshvp;
141 m_eps_pv = threshpv;
142 m_eps_pp = threshpp;
143 m_eps = std::min(std::min(m_eps_vv,m_eps_vp),std::min(m_eps_pv,m_eps_pp));
144 }
145
147 void set_threshold(number threshvv,number thresh_vp_pv_pp)
148 {
149 m_eps_vv = threshvv;
150 m_eps_vp = thresh_vp_pv_pp;
151 m_eps_pv = thresh_vp_pv_pp;
152 m_eps_pp = thresh_vp_pv_pp;
153 m_eps = std::min(m_eps_vv,m_eps_vp);
154 }
155
157 {
158 m_eps_vv = thresh;
159 m_eps_vp = thresh;
160 m_eps_pv = thresh;
161 m_eps_pp = thresh;
162 m_eps = thresh;
163 }
164
166 void set_info(bool info)
167 {
168 m_info = info;
169 }
170
171
172protected:
173 // Name of preconditioner
174 virtual const char* name() const {return "PCRILUTPreconditioner";}
175
176 // Preprocess routine
178 {
179 STATIC_ASSERT(matrix_type::rows_sorted, Matrix_has_to_have_sorted_rows);
180 matrix_type& mat = *pOp;
181
182#ifdef UG_PARALLEL
183 // copy original matrix
184 // commented-out function seems to be broken (A.Vogel)
185 // MakeConsistent(mat, m_ILUT);
186
187 // we use new instead (A.Vogel)
188 m_A = mat;
190
191 // set zero on slaves
192 std::vector<IndexLayout::Element> vIndex;
193 CollectUniqueElements(vIndex, m_A.layouts()->slave());
194 SetDirichletRow(m_A, vIndex);
195#else
196 // copy original matrix
197 m_A = mat;
198#endif
199 typedef typename matrix_type::connection connection;
200 m_L.resize_and_clear(m_A.num_rows(), m_A.num_cols());
201 m_U.resize_and_clear(m_A.num_rows(), m_A.num_cols());
202
203 // con is the current line of L/U
204 std::vector<typename matrix_type::connection> con;
205 con.reserve(300);
206 con.resize(0);
207
208 static const size_t velocity=0;
209 static const size_t pressure=1;
210
211 // init row 0 of U
212 for(typename matrix_type::row_iterator i_it = m_A.begin_row(0); i_it != m_A.end_row(0); ++i_it)
213 con.push_back(connection(i_it.index(), i_it.value()));
214 m_U.set_matrix_row(0, &con[0], con.size());
215
216 size_t totalentries=0;
217 size_t maxentries=0;
218
219 for(size_t i=1; i<m_A.num_rows(); i++)
220 {
221 size_t itype;
222 if (m_A(i,i)!=0) itype=velocity; else itype=pressure;
223
224 con.resize(0);
225 size_t u_part=0;
226
227 // get the row A(i, .) into con
228 double dmax=0;
229 m_remove_zeros = false;
230 if (m_remove_zeros==false){
231 for(typename matrix_type::row_iterator i_it = m_A.begin_row(i); i_it != m_A.end_row(i); ++i_it)
232 {
233 con.push_back(connection(i_it.index(), i_it.value()));
234 if(dmax < BlockNorm(i_it.value()))
235 dmax = BlockNorm(i_it.value());
236 }
237 } else {
238 typename matrix_type::row_iterator i_it = m_A.begin_row(i);
239 con.push_back(connection(i_it.index(), i_it.value()));
240 ++i_it;
241 for(; i_it != m_A.end_row(i); ++i_it)
242 {
243 if (i_it.value()==0) continue;
244 con.push_back(connection(i_it.index(), i_it.value()));
245 if(dmax < BlockNorm(i_it.value()))
246 dmax = BlockNorm(i_it.value());
247 }
248 }
249 // eliminate all entries A(i, k) with k<i with rows U(k, .) and k<i
250 for(size_t i_it = 0; i_it < con.size(); ++i_it)
251 {
252 size_t k = con[i_it].iIndex;
253 if(k >= i)
254 {
255 // safe where U begins / L ends in con
256 u_part = i_it;
257 break;
258 }
259 if(con[i_it].dValue == 0.0) continue;
260 UG_ASSERT(m_U.num_connections(k) != 0 && m_U.begin_row(k).index() == k, "");
261 block_type &ukk = m_U.begin_row(k).value();
262
263 // add row k to row i by A(i, .) -= U(k,.) A(i,k) / U(k,k)
264 // so that A(i,k) is zero.
265 // safe A(i,k)/U(k,k) in con, (later L(i,k) )
266 block_type &d = con[i_it].dValue = con[i_it].dValue / ukk;
267
268 typename matrix_type::row_iterator k_it = m_U.begin_row(k); // upper row iterator
269 ++k_it; // skip diag
270 size_t j = i_it+1;
271 while(k_it != m_U.end_row(k) && j < con.size())
272 {
273
274 // (since con and U[k] is sorted, we can do sth like a merge on the two lists)
275 if(k_it.index() == con[j].iIndex)
276 {
277 // match
278 con[j].dValue -= k_it.value() * d;
279 ++k_it; ++j;
280 }
281 else if(k_it.index() < con[j].iIndex)
282 {
283 // we have a value in U(k, (*k_it).iIndex), but not in A.
284 // check tolerance criteria
285
286 typename matrix_type::connection
287 c(k_it.index(), k_it.value() * d * -1.0);
288
289 if(BlockNorm(c.dValue) > dmax * m_eps){
290
291 size_t kitType;
292 if (m_A(k_it.index(),k_it.index())!=0) kitType=velocity; else kitType=pressure;
293
294 if ((itype==velocity)&&(kitType==velocity)){
295 if(BlockNorm(c.dValue) > dmax * m_eps_vv)
296 {
297 // insert sorted
298 con.insert(con.begin()+j, c);
299 ++j;
300 }
301 }
302 if ((itype==velocity)&&(kitType==pressure)){
303 if(BlockNorm(c.dValue) > dmax * m_eps_vp)
304 {
305 // insert sorted
306 con.insert(con.begin()+j, c);
307 ++j;
308 }
309 }
310 if ((itype==pressure)&&(kitType==velocity)){
311 if(BlockNorm(c.dValue) > dmax * m_eps_pv)
312 {
313 // insert sorted
314 con.insert(con.begin()+j, c);
315 ++j;
316 }
317 }
318 if ((itype==pressure)&&(kitType==pressure)){
319 if(BlockNorm(c.dValue) > dmax * m_eps_pp)
320 {
321 // insert sorted
322 con.insert(con.begin()+j, c);
323 ++j;
324 }
325 }
326 }
327 // else do some lumping
328 ++k_it;
329 }
330 else
331 {
332 // we have a value in A(k, con[j].iIndex), but not in U.
333 ++j;
334 }
335 }
336 // insert new connections after last connection of row i
337 if (k_it!=m_U.end_row(k)){
338 for (;k_it!=m_U.end_row(k);++k_it){
339 typename matrix_type::connection c(k_it.index(),-k_it.value()*d);
340 if(BlockNorm(c.dValue) > dmax * m_eps){
341 size_t kitType;
342 if (m_A(k_it.index(),k_it.index())!=0) kitType=velocity; else kitType=pressure;
343 if ((itype==velocity)&&(kitType==velocity)){
344 if(BlockNorm(c.dValue) > dmax * m_eps_vv)
345 {
346 con.push_back(c);
347 }
348 } else {
349 if ((itype==velocity)&&(kitType==pressure)){
350 if(BlockNorm(c.dValue) > dmax * m_eps_vp)
351 {
352 con.push_back(c);
353 }
354 } else {
355 if ((itype==pressure)&&(kitType==velocity)){
356 if(BlockNorm(c.dValue) > dmax * m_eps_pv)
357 {
358 con.push_back(c);
359 }
360 } else {
361 // if ((itype==pressure)&&(kitType==pressure)){
362 if(BlockNorm(c.dValue) > dmax * m_eps_pp)
363 {
364 con.push_back(c);
365 }
366 // }
367 }
368 }
369 }
370
371 };
372 }
373 };
374 }
375
376 totalentries+=con.size();
377 if(maxentries < con.size()) maxentries = con.size();
378
379 // safe L and U
380 m_L.set_matrix_row(i, &con[0], u_part);
381 m_U.set_matrix_row(i, &con[u_part], con.size()-u_part);
382
383 }
384
385 m_L.defragment();
386 m_U.defragment();
387 m_A.defragment();
388
389 if (m_info==true){
390 UG_LOG("\n PCRILUT storage information:\n");
391 UG_LOG(" A nr of connections: " << m_A.total_num_connections() << "\n");
392 UG_LOG(" L+U nr of connections: " << m_L.total_num_connections()+m_U.total_num_connections() << "\n");
393 UG_LOG(" Increase factor: " << (float)(m_L.total_num_connections() + m_U.total_num_connections() )/m_A.total_num_connections() << "\n");
394 }
395
396 return true;
397 }
398
399 // Stepping routine
401 {
402#ifdef UG_PARALLEL
403 SmartPtr<vector_type> spDtmp = d.clone();
404 spDtmp->change_storage_type(PST_UNIQUE);
405 // apply iterator: c = LU^{-1}*d (damp is not used)
406 // L
407 for(size_t i=0; i < m_L.num_rows(); i++)
408 {
409 // c[i] = d[i] - m_L[i]*c;
410 c[i] = (*spDtmp)[i];
411 for(typename matrix_type::row_iterator it = m_L.begin_row(i); it != m_L.end_row(i); ++it)
412 MatMultAdd(c[i], 1.0, c[i], -1.0, it.value(), c[it.index()] );
413 // lii = 1.0.
414 }
415 // U
416 //
417 // last row diagonal U entry might be close to zero with corresponding zero rhs
418 // when solving Navier Stokes system, therefore handle separately
419 {
420 size_t i=m_U.num_rows()-1;
421 typename matrix_type::row_iterator it = m_U.begin_row(i);
422 UG_ASSERT(it.index() == i, "");
423 block_type &uii = it.value();
424 typename vector_type::value_type s = c[i];
425 // check if close to zero
426 if (BlockNorm(uii)<m_small_lower){
427 // set correction to zero
428 c[i] = 0;
430 UG_LOG("Warning: zero entry in last row of U with corresponding non-zero rhs entry (" << BlockNorm(s) << ")\n");
431 }
432 } else {
433 // c[i] = s/uii;
434 InverseMatMult(c[i], 1.0, uii, s);
435 }
436 }
437 // handle all other rows
438 for(size_t i=m_U.num_rows()-2; ; i--)
439 {
440 typename matrix_type::row_iterator it = m_U.begin_row(i);
441 UG_ASSERT(it.index() == i, "");
442 block_type &uii = it.value();
443
444 typename vector_type::value_type s = c[i];
445 ++it; // skip diag
446 for(; it != m_U.end_row(i); ++it)
447 // s -= it.value() * c[it.index()];
448 MatMultAdd(s, 1.0, s, -1.0, it.value(), c[it.index()] );
449
450 // c[i] = s/uii;
451 InverseMatMult(c[i], 1.0, uii, s);
452
453 if(i==0) break;
454 }
455 // Correction is always consistent
456 c.set_storage_type(PST_ADDITIVE);
457 c.change_storage_type(PST_CONSISTENT);
458#else
459 // apply iterator: c = LU^{-1}*d (damp is not used)
460 // L
461 for(size_t i=0; i < m_L.num_rows(); i++)
462 {
463 // c[i] = d[i] - m_L[i]*c;
464 c[i] = d[i];
465 for(typename matrix_type::row_iterator it = m_L.begin_row(i); it != m_L.end_row(i); ++it)
466 MatMultAdd(c[i], 1.0, c[i], -1.0, it.value(), c[it.index()] );
467 // lii = 1.0.
468 }
469 // U
470 //
471 // last row diagonal U entry might be close to zero with corresponding zero rhs
472 // when solving Navier Stokes system, therefore handle separately
473 {
474 size_t i=m_U.num_rows()-1;
475 typename matrix_type::row_iterator it = m_U.begin_row(i);
476 UG_ASSERT(it.index() == i, "");
477 block_type &uii = it.value();
478 typename vector_type::value_type s = c[i];
479 // check if close to zero
480 if (BlockNorm(uii)<m_small_lower){
481 // set correction to zero
482 c[i] = 0;
484 UG_LOG("Warning: zero entry in last row of U with corresponding non-zero rhs entry (" << BlockNorm(s) << ")\n");
485 }
486 } else {
487 // c[i] = s/uii;
488 InverseMatMult(c[i], 1.0, uii, s);
489 }
490 }
491 // handle all other rows
492 for(size_t i=m_U.num_rows()-2; ; i--)
493 {
494 typename matrix_type::row_iterator it = m_U.begin_row(i);
495 UG_ASSERT(it.index() == i, "");
496 block_type &uii = it.value();
497
498 typename vector_type::value_type s = c[i];
499 ++it; // skip diag
500 for(; it != m_U.end_row(i); ++it)
501 // s -= it.value() * c[it.index()];
502 MatMultAdd(s, 1.0, s, -1.0, it.value(), c[it.index()] );
503
504 // c[i] = s/uii;
505 InverseMatMult(c[i], 1.0, uii, s);
506
507 if(i==0) break;
508 }
509#endif
510 return true;
511 }
512
513 // Postprocess routine
514 virtual bool postprocess() {return true;}
515
516protected:
525 bool m_info;
527 static const number m_small_lower;
528 static const number m_small_upper;
529};
530
531template <typename TAlgebra>
533template <typename TAlgebra>
535
536} // end namespace ug
537
538#endif
parameterString s
Definition Biogas.lua:2
virtual void set_debug(SmartPtr< IDebugWriter< algebra_type > > spDebugWriter)
SmartPtr< IDebugWriter< algebra_type > > debug_writer()
SmartPtr< IDamping< X, Y > > damping()
PCRILUTPreconditioner class.
Definition pcr_ilut.h:60
void set_threshold(number threshvv, number threshvp, number threshpv, number threshpp)
sets threshold for incomplete LU factorisation
Definition pcr_ilut.h:137
TAlgebra::vector_type vector_type
Definition pcr_ilut.h:66
bool m_info
Definition pcr_ilut.h:525
number m_eps
Definition pcr_ilut.h:520
number m_eps_vv
Definition pcr_ilut.h:521
bool m_remove_zeros
Definition pcr_ilut.h:526
PCRILUTPreconditioner(number threshvv, number threshvp, number threshpv, number threshpp, bool info)
Definition pcr_ilut.h:110
virtual bool preprocess(SmartPtr< MatrixOperator< matrix_type, vector_type > > pOp)
Definition pcr_ilut.h:177
PCRILUTPreconditioner(number threshvv, number thresh_vp_pv_pp)
Definition pcr_ilut.h:87
virtual bool postprocess()
Definition pcr_ilut.h:514
PCRILUTPreconditioner(double eps=1e-6)
Definition pcr_ilut.h:81
number m_eps_pp
Definition pcr_ilut.h:524
IPreconditioner< TAlgebra >::matrix_operator_type matrix_operator_type
Matrix Operator type.
Definition pcr_ilut.h:72
void set_threshold(number threshvv, number thresh_vp_pv_pp)
sets threshold for incomplete LU factorisation
Definition pcr_ilut.h:147
PCRILUTPreconditioner(double eps, bool info)
Definition pcr_ilut.h:99
TAlgebra::matrix_type matrix_type
Definition pcr_ilut.h:69
matrix_type m_U
Definition pcr_ilut.h:518
matrix_type m_A
Definition pcr_ilut.h:519
virtual const char * name() const
Definition pcr_ilut.h:174
void set_info(bool info)
sets storage information output to true or false
Definition pcr_ilut.h:166
static const number m_small_lower
Definition pcr_ilut.h:527
matrix_type m_L
Definition pcr_ilut.h:517
number m_eps_pv
Definition pcr_ilut.h:523
number m_eps_vp
Definition pcr_ilut.h:522
void set_threshold(number thresh)
Definition pcr_ilut.h:156
TAlgebra algebra_type
Definition pcr_ilut.h:63
PCRILUTPreconditioner(number threshvv, number threshvp, number threshpv, number threshpp)
Definition pcr_ilut.h:93
matrix_type::value_type block_type
Definition pcr_ilut.h:75
virtual bool supports_parallel() const
returns if parallel solving is supported
Definition pcr_ilut.h:134
SmartPtr< ILinearIterator< vector_type > > clone()
Definition pcr_ilut.h:118
virtual bool step(SmartPtr< MatrixOperator< matrix_type, vector_type > > pOp, vector_type &c, const vector_type &d)
Definition pcr_ilut.h:400
virtual ~PCRILUTPreconditioner()
Definition pcr_ilut.h:129
static const number m_small_upper
Definition pcr_ilut.h:528
PCRILUTPreconditioner(number threshvv, number thresh_vp_pv_pp, bool info)
Definition pcr_ilut.h:104
void MatAddSlaveRowsToMasterRowOverlap0(TMatrix &mat)
PST_CONSISTENT
PST_UNIQUE
PST_ADDITIVE
bool InverseMatMult(DenseVector< FixedArray1< double, 1 > > &dest, double beta1, const DenseMatrix< FixedArray2< double, 1, 1 > > &A1, const DenseVector< FixedArray1< double, 1 > > &w1)
void MatMultAdd(DenseVector< vector_t > &dest, const number &alpha1, const DenseVector< vector_t > &v1, const number &beta1, const DenseMatrix< matrix_t > &A1, const DenseVector< vector_t > &w1)
#define UG_ASSERT(expr, msg)
#define STATIC_ASSERT(expr, msg)
#define UG_LOG(msg)
double number
void SetDirichletRow(TMatrix &mat, const DoFIndex &ind)
number BlockNorm(const number &a)