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