ug4
Loading...
Searching...
No Matches
pilut.h
Go to the documentation of this file.
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__PILUT__
34#define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__PILUT__
35
38#ifdef UG_PARALLEL
40#endif
41#include "common/progress.h"
43namespace ug{
44
45template <typename TAlgebra>
46class PILUTPreconditioner : public IPreconditioner<TAlgebra>
47{
48 public:
49 // Algebra type
50 typedef TAlgebra algebra_type;
51
52 // Vector type
53 typedef typename TAlgebra::vector_type vector_type;
54
55 // Matrix type
56 typedef typename TAlgebra::matrix_type matrix_type;
57
60
61 private:
62 typedef typename matrix_type::value_type block_type;
63 using IPreconditioner<TAlgebra>::debug_writer;
64 using IPreconditioner<TAlgebra>::set_debug;
65
66 public:
67 // Constructor
68 PILUTPreconditioner(double eps=0) :
69 m_eps(eps), m_info(false), m_show_progress(true)
70 {};
71
72 // Clone
73
75 {
77 newInst->set_debug(debug_writer());
78 newInst->set_damp(this->damping());
79 newInst->set_info(m_info);
80 return newInst;
81 }
82
83 // Destructor
85 {
86 };
87
89 void set_threshold(number thresh)
90 {
91 m_eps = thresh;
92 }
93
95 void set_info(bool info)
96 {
97 m_info = info;
98 }
99
101 void set_show_progress(bool s)
102 {
104 }
105
106 protected:
107 // Name of preconditioner
108 virtual const char* name() const {return "PILUTPreconditioner";}
109
110 // Preprocess routine
112 {
113 matrix_type mat;
114 mat.set_as_copy_of(*pOp);
115 STATIC_ASSERT(matrix_type::rows_sorted, Matrix_has_to_have_sorted_rows);
116
117 // Prepare Inverse Matrix
118 matrix_type* A = &mat;
119 typedef typename matrix_type::connection connection;
120 m_L.resize_and_clear(A->num_rows(), A->num_cols());
121 m_U.resize_and_clear(A->num_rows(), A->num_cols());
122
123 // con is the current line of L/U
124 // i also tried using std::list here or a special custom vector-based linked list
125 // but vector is fastest, even with the insert operation.
126 std::vector<typename matrix_type::connection> con;
127 con.reserve(300);
128 con.resize(0);
129
130 // init row 0 of U
131 for(typename matrix_type::row_iterator i_it = A->begin_row(0); i_it != A->end_row(0); ++i_it)
132 con.push_back(connection(i_it.index(), i_it.value()));
133 m_U.set_matrix_row(0, &con[0], con.size());
134
135 size_t totalentries=0;
136 size_t maxentries=0;
137 Progress prog;
139 PROGRESS_START_WITH(prog, A->num_rows(),
140 "Using ILUT(" << m_eps << ") on " << A->num_rows() << " x " << A->num_rows() << " matrix...");
141
142 PROFILE_BEGIN(PILUT1);
143 for(size_t i=1; i<A->num_rows()/2; i++)
144 {
145 if(m_show_progress) {PROGRESS_UPDATE(prog, i);}
146 con.resize(0);
147 size_t u_part=0;
148
149 // get the row A(i, .) into con
150 double dmax=0;
151 for(typename matrix_type::row_iterator i_it = A->begin_row(i); i_it != A->end_row(i); ++i_it)
152 {
153 con.push_back(connection(i_it.index(), i_it.value()));
154 if(dmax < BlockNorm(i_it.value()))
155 dmax = BlockNorm(i_it.value());
156 }
157
158
159 u_part = eliminate_row(con, 0, i, dmax);
160 totalentries+=con.size();
161 if(maxentries < con.size()) maxentries = con.size();
162 // safe L and U
163 for(size_t i_it=0; i_it<con.size(); i_it++)
164 {
165 if(con[i_it].iIndex < i) m_L(i, con[i_it].iIndex) = con[i_it].dValue;
166 else m_U(i, con[i_it].iIndex) = con[i_it].dValue;
167 }
168 }
169 PROFILE_END();
170 PROFILE_BEGIN(PILUT2);
171 for(size_t i=A->num_rows()/2; i<A->num_rows(); i++)
172 {
173 con.resize(0);
174 size_t u_part=0;
175
176 // get the row A(i, .) into con
177 double dmax=0;
178 for(typename matrix_type::row_iterator i_it = A->begin_row(i); i_it != A->end_row(i); ++i_it)
179 {
180 con.push_back(connection(i_it.index(), i_it.value()));
181 if(dmax < BlockNorm(i_it.value()))
182 dmax = BlockNorm(i_it.value());
183 }
184
185
186 u_part = eliminate_row(con, 0, A->num_rows()/2, dmax);
187 totalentries+=con.size();
188 if(maxentries < con.size()) maxentries = con.size();
189 // safe L and U
190 A->set_matrix_row(i, &con[0], con.size());
191 }
192 con.clear();
193 size_t u_part;
194 for(typename matrix_type::row_iterator i_it = A->begin_row(A->num_rows()/2); i_it != A->end_row(A->num_rows()/2); ++i_it)
195 {
196 if(i_it.index() <= A->num_rows()/2)
197 u_part=con.size();
198 con.push_back(connection(i_it.index(), i_it.value()));
199 }
200 m_L.set_matrix_row(A->num_rows()/2, &con[0], u_part);
201 m_U.set_matrix_row(A->num_rows()/2, &con[u_part], con.size()-u_part);
202 PROFILE_END();
203 PROFILE_BEGIN(PILUT3);
204 for(size_t i=A->num_rows()/2+1; i<A->num_rows(); i++)
205 {
206 if(m_show_progress) {PROGRESS_UPDATE(prog, i);}
207 con.resize(0);
208 size_t u_part=0;
209
210 // get the row A(i, .) into con
211 double dmax=0;
212 for(typename matrix_type::row_iterator i_it = A->begin_row(i); i_it != A->end_row(i); ++i_it)
213 {
214 con.push_back(connection(i_it.index(), i_it.value()));
215 if(dmax < BlockNorm(i_it.value()))
216 dmax = BlockNorm(i_it.value());
217 }
218
219
220 u_part = eliminate_row(con, A->num_rows()/2, i, dmax);
221 totalentries+=con.size();
222 if(maxentries < con.size()) maxentries = con.size();
223 // safe L and U
224 for(size_t i_it=0; i_it<con.size(); i_it++)
225 {
226 if(con[i_it].iIndex < i) m_L(i, con[i_it].iIndex) = con[i_it].dValue;
227 else m_U(i, con[i_it].iIndex) = con[i_it].dValue;
228 }
229 }
230 PROFILE_END();
231 //m_L.print();
232 //m_U.print();
233
235
236 m_L.defragment();
237 m_U.defragment();
238
239 if (m_info==true){
240 UG_LOG("\n ILUT storage information:\n");
241 UG_LOG(" A nr of connections: " << A->total_num_connections() << "\n");
242 UG_LOG(" L+U nr of connections: " << m_L.total_num_connections()+m_U.total_num_connections() << "\n");
243 UG_LOG(" Increase factor: " << (float)(m_L.total_num_connections() + m_U.total_num_connections() )/A->total_num_connections() << "\n");
244 UG_LOG(reset_floats << "Total entries: " << totalentries << " (" << ((double)totalentries) / (A->num_rows()*A->num_rows()) << "% of dense)");
245 }
246
247 return true;
248 }
249
250 size_t eliminate_row(std::vector<typename matrix_type::connection> &con, size_t start, size_t stop, double dmax)
251 {
252 size_t u_part;
253 // eliminate all entries A(i, k) with k<i with rows U(k, .) and k<i
254 for(size_t i_it = 0; i_it < con.size(); ++i_it)
255 {
256 size_t k = con[i_it].iIndex;
257 if(k < start) continue;
258 if(k >= stop)
259 {
260 // safe where U begins / L ends in con
261 u_part = i_it;
262 break;
263 }
264 if(con[i_it].dValue == 0.0) continue;
265 UG_ASSERT(m_U.num_connections(k) != 0 && m_U.begin_row(k).index() == k, "");
266 block_type &ukk = m_U.begin_row(k).value();
267
268 // add row k to row i by A(i, .) -= U(k,.) A(i,k) / U(k,k)
269 // so that A(i,k) is zero.
270 // safe A(i,k)/U(k,k) in con, (later L(i,k) )
271 block_type &d = con[i_it].dValue = con[i_it].dValue / ukk;
272
273 typename matrix_type::row_iterator k_it = m_U.begin_row(k); // upper row iterator
274 ++k_it; // skip diag
275 size_t j = i_it+1;
276 while(k_it != m_U.end_row(k) && j < con.size())
277 {
278 // (since con and U[k] is sorted, we can do sth like a merge on the two lists)
279 if(k_it.index() == con[j].iIndex)
280 {
281 // match
282 con[j].dValue -= k_it.value() * d;
283 ++k_it; ++j;
284 }
285 else if(k_it.index() < con[j].iIndex)
286 {
287 // we have a value in U(k, (*k_it).iIndex), but not in A.
288 // check tolerance criteria
289
290 typename matrix_type::connection
291 c(k_it.index(), k_it.value() * d * -1.0);
292
293 if(BlockNorm(c.dValue) > dmax * m_eps)
294 {
295 // insert sorted
296 con.insert(con.begin()+j, c);
297 ++j; // don't do this when using iterators
298 }
299 // else do some lumping
300 ++k_it;
301 }
302 else
303 {
304 // we have a value in A(k, con[j].iIndex), but not in U.
305 ++j;
306 }
307 }
308 // insert new connections after last connection of row i
309 if (k_it!=m_U.end_row(k)){
310 for (;k_it!=m_U.end_row(k);++k_it){
311 typename matrix_type::connection c(k_it.index(),-k_it.value()*d);
312 if(BlockNorm(c.dValue) > dmax * m_eps)
313 {
314 con.push_back(c);
315 }
316 }
317 };
318 }
319 return u_part;
320
321 }
322
323 // Stepping routine
325 {
326 // apply iterator: c = LU^{-1}*d (damp is not used)
327 // L
328 for(size_t i=0; i < m_L.num_rows(); i++)
329 {
330 // c[i] = d[i] - m_L[i]*c;
331 c[i] = d[i];
332 for(typename matrix_type::row_iterator it = m_L.begin_row(i); it != m_L.end_row(i); ++it)
333 MatMultAdd(c[i], 1.0, c[i], -1.0, it.value(), c[it.index()] );
334 // lii = 1.0.
335 }
336 // U
337 //
338 // last row diagonal U entry might be close to zero with corresponding zero rhs
339 // when solving Navier Stokes system, therefore handle separately
340 {
341 size_t i=m_U.num_rows()-1;
342 typename matrix_type::row_iterator it = m_U.begin_row(i);
343 UG_ASSERT(it.index() == i, "");
344 block_type &uii = it.value();
345 typename vector_type::value_type s = c[i];
346 // check if close to zero
347 if (BlockNorm(uii)<m_small){
348 // set correction to zero
349 c[i] = 0;
350 if (BlockNorm(s)>m_small){
351 UG_LOG("Warning: near-zero diagonal entry in last row of U with corresponding non-near-zero rhs entry (" << BlockNorm(s) << ")\n");
352 }
353 } else {
354 // c[i] = s/uii;
355 InverseMatMult(c[i], 1.0, uii, s);
356 }
357 }
358 // handle all other rows
359 for(size_t i=m_U.num_rows()-2; ; i--)
360 {
361 typename matrix_type::row_iterator it = m_U.begin_row(i);
362 UG_ASSERT(it.index() == i, "");
363 block_type &uii = it.value();
364
365 typename vector_type::value_type s = c[i];
366 ++it; // skip diag
367 for(; it != m_U.end_row(i); ++it)
368 // s -= it.value() * c[it.index()];
369 MatMultAdd(s, 1.0, s, -1.0, it.value(), c[it.index()] );
370
371 // c[i] = s/uii;
372 InverseMatMult(c[i], 1.0, uii, s);
373
374 if(i==0) break;
375 }
376
377#ifdef UG_PARALLEL
378 // Correction is always consistent
379 // todo: We set here correction to consistent, but it is not. Think about how to use ilu in parallel.
380 c.set_storage_type(PST_CONSISTENT);
381#endif
382 return true;
383 }
384
385 // Postprocess routine
386 virtual bool postprocess() {return true;}
387
388 protected:
391 double m_eps;
392 bool m_info;
394 static const number m_small;
395};
396
397// define constant
398template <typename TAlgebra>
400
401} // end namespace ug
402
403#endif
parameterString s
Definition smart_pointer.h:108
virtual void set_debug(SmartPtr< IDebugWriter< algebra_type > > spDebugWriter)
set debug writer
Definition debug_writer.h:384
SmartPtr< IDebugWriter< algebra_type > > debug_writer()
returns the debug writer
Definition debug_writer.h:391
SmartPtr< IDamping< TAlgebra::vector_type, TAlgebra::vector_type > > damping()
returns the scaling
Definition linear_iterator.h:173
describes a linear iterator that is based on a matrix operator
Definition preconditioner.h:103
Definition pilut.h:47
virtual const char * name() const
returns the name of iterator
Definition pilut.h:108
matrix_type m_U
Definition pilut.h:390
double m_eps
Definition pilut.h:391
matrix_type m_L
Definition pilut.h:389
virtual bool preprocess(SmartPtr< MatrixOperator< matrix_type, vector_type > > pOp)
initializes the preconditioner
Definition pilut.h:111
TAlgebra::matrix_type matrix_type
Definition pilut.h:56
void set_info(bool info)
sets storage information output to true or false
Definition pilut.h:95
void set_threshold(number thresh)
sets threshold for incomplete LU factorisation (added 01122010ih)
Definition pilut.h:89
virtual ~PILUTPreconditioner()
Definition pilut.h:84
PILUTPreconditioner(double eps=0)
Definition pilut.h:68
TAlgebra::vector_type vector_type
Definition pilut.h:53
size_t eliminate_row(std::vector< typename matrix_type::connection > &con, size_t start, size_t stop, double dmax)
Definition pilut.h:250
TAlgebra algebra_type
Definition pilut.h:50
bool m_info
Definition pilut.h:392
bool m_show_progress
Definition pilut.h:393
void set_show_progress(bool s)
set whether the progress meter should be shown
Definition pilut.h:101
SmartPtr< ILinearIterator< vector_type > > clone()
clone
Definition pilut.h:74
static const number m_small
Definition pilut.h:394
virtual bool step(SmartPtr< MatrixOperator< matrix_type, vector_type > > pOp, vector_type &c, const vector_type &d)
computes a new correction c = B*d
Definition pilut.h:324
matrix_type::value_type block_type
Definition pilut.h:62
IPreconditioner< TAlgebra >::matrix_operator_type matrix_operator_type
Matrix Operator type.
Definition pilut.h:59
virtual bool postprocess()
cleans the operator
Definition pilut.h:386
Definition progress.h:50
@ PST_CONSISTENT
Definition parallel_storage_type.h:68
static std::ios_base & reset_floats(std::ios_base &o)
used to unset flags set by std::scientific or std::fixed.
Definition ostream_util.h:52
#define UG_ASSERT(expr, msg)
Definition assert.h:70
#define STATIC_ASSERT(expr, msg)
Definition common.h:61
#define UG_LOG(msg)
Definition log.h:367
double number
Definition types.h:124
the ug namespace
bool MatMultAdd(vector_t &dest, const number &alpha1, const vector_t &v1, const number &beta1, const matrix_t &A1, const vector_t &w1)
calculates dest = alpha1*v1 + beta1 * A1 *w1;
Definition operations_mat.h:68
bool InverseMatMult(number &dest, const double &beta, const TMat &mat, const TVec &vec)
you can implement this function with GetInverse and MatMult
double BlockNorm(const TYPE &v)
Definition blocks.h:57
#define PROFILE_BEGIN(name)
Definition profiler.h:254
#define PROFILE_END()
Definition profiler.h:256
#define PROGRESS_FINISH(progVarName)
Definition progress.h:118
#define PROGRESS_START_WITH(progVarName, dSize, msg)
Definition progress.h:114
#define PROGRESS_UPDATE(progVarName, d)
Definition progress.h:117