Line data Source code
1 :
2 : #ifndef __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__ILUT__
3 : #define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__ILUT__
4 :
5 : #include "common/util/smart_pointer.h"
6 : #include "lib_algebra/operator/interface/preconditioner.h"
7 : #ifdef UG_PARALLEL
8 : #include "lib_algebra/parallelization/parallelization.h"
9 : #endif
10 : #include "common/progress.h"
11 : #include "common/util/ostream_util.h"
12 : #include "common/profiler/profiler.h"
13 :
14 : #include "lib_algebra/algebra_common/vector_util.h"
15 :
16 : #include "lib_algebra/ordering_strategies/algorithms/IOrderingAlgorithm.h"
17 : #include "lib_algebra/ordering_strategies/algorithms/native_cuthill_mckee.h" // for backward compatibility
18 :
19 : #include "lib_algebra/algebra_common/permutation_util.h"
20 :
21 : namespace ug{
22 :
23 : template <typename TAlgebra>
24 : class ILUTPreconditioner : public IPreconditioner<TAlgebra>
25 : {
26 : public:
27 : // Algebra type
28 : typedef TAlgebra algebra_type;
29 :
30 : // Vector type
31 : typedef typename TAlgebra::vector_type vector_type;
32 : typedef typename vector_type::value_type vector_value;
33 :
34 : // Matrix type
35 : typedef typename TAlgebra::matrix_type matrix_type;
36 :
37 : typedef typename matrix_type::row_iterator matrix_row_iterator;
38 : typedef typename matrix_type::const_row_iterator const_matrix_row_iterator;
39 : typedef typename matrix_type::connection matrix_connection;
40 :
41 : /// Ordering type
42 : typedef std::vector<size_t> ordering_container_type;
43 : typedef IOrderingAlgorithm<TAlgebra, ordering_container_type> ordering_algo_type;
44 :
45 : /// Matrix Operator type
46 : typedef typename IPreconditioner<TAlgebra>::matrix_operator_type matrix_operator_type;
47 : using IPreconditioner<TAlgebra>::set_debug;
48 :
49 : protected:
50 : typedef typename matrix_type::value_type block_type;
51 :
52 : using IPreconditioner<TAlgebra>::debug_writer;
53 : using IPreconditioner<TAlgebra>::write_debug;
54 :
55 : private:
56 : typedef IPreconditioner<TAlgebra> base_type;
57 :
58 : public:
59 : /// Constructor
60 0 : ILUTPreconditioner(double eps=1e-6)
61 0 : : m_eps(eps), m_info(false), m_show_progress(true), m_bSortIsIdentity(false)
62 : {
63 : //default was set true
64 0 : m_spOrderingAlgo = make_sp(new NativeCuthillMcKeeOrdering<TAlgebra, ordering_container_type>());
65 0 : };
66 :
67 : /// clone constructor
68 0 : ILUTPreconditioner(const ILUTPreconditioner<TAlgebra> &parent)
69 0 : : base_type(parent), m_spOrderingAlgo(parent.m_spOrderingAlgo)
70 : {
71 0 : m_eps = parent.m_eps;
72 0 : set_info(parent.m_info);
73 0 : m_bSortIsIdentity = parent.m_bSortIsIdentity;
74 0 : }
75 :
76 : /// Clone
77 0 : virtual SmartPtr<ILinearIterator<vector_type> > clone()
78 : {
79 0 : return make_sp(new ILUTPreconditioner<algebra_type>(*this));
80 : }
81 :
82 :
83 : /// Destructor
84 0 : virtual ~ILUTPreconditioner()
85 0 : {};
86 :
87 : /// returns if parallel solving is supported
88 0 : virtual bool supports_parallel() const {return true;}
89 :
90 : /// sets threshold for incomplete LU factorisation (added 01122010ih)
91 0 : void set_threshold(number thresh)
92 : {
93 0 : m_eps = thresh;
94 0 : }
95 :
96 : /// sets storage information output to true or false
97 0 : void set_info(bool info)
98 : {
99 0 : m_info = info;
100 0 : }
101 :
102 : /// set whether the progress meter should be shown
103 : void set_show_progress(bool s)
104 : {
105 0 : m_show_progress = s;
106 : }
107 :
108 0 : virtual std::string config_string() const
109 : {
110 0 : std::stringstream ss;
111 0 : ss << "ILUT(threshold = " << m_eps << ", sort = " << (m_spOrderingAlgo.valid()?"true":"false") << ")";
112 0 : if(m_eps == 0.0) ss << " = Sparse LU";
113 0 : return ss.str();
114 0 : }
115 :
116 : /// sets an ordering algorithm
117 0 : void set_ordering_algorithm(SmartPtr<ordering_algo_type> ordering_algo){
118 0 : m_spOrderingAlgo = ordering_algo;
119 0 : }
120 :
121 : /// set cuthill-mckee sort on/off
122 0 : void set_sort(bool b)
123 : {
124 0 : if(b){
125 0 : m_spOrderingAlgo = make_sp(new NativeCuthillMcKeeOrdering<TAlgebra, ordering_container_type>());
126 : }
127 : else{
128 0 : m_spOrderingAlgo = SPNULL;
129 : }
130 :
131 : UG_LOG("\nILUT: please use 'set_ordering_algorithm(..)' in the future\n");
132 0 : }
133 :
134 :
135 : protected:
136 : // Name of preconditioner
137 0 : virtual const char* name() const {return "ILUT";}
138 :
139 : protected:
140 0 : virtual bool init(SmartPtr<ILinearOperator<vector_type> > J,
141 : const vector_type& u)
142 : {
143 : // cast to matrix based operator
144 0 : SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp =
145 : J.template cast_dynamic<MatrixOperator<matrix_type, vector_type> >();
146 :
147 : // Check that matrix if of correct type
148 0 : if(pOp.invalid())
149 0 : UG_THROW(name() << "::init': Passed Operator is "
150 : "not based on matrix. This Preconditioner can only "
151 : "handle matrix-based operators.");
152 :
153 0 : m_u = &u;
154 :
155 : // forward request to matrix based implementation
156 0 : return base_type::init(pOp);
157 : }
158 :
159 0 : bool init(SmartPtr<ILinearOperator<vector_type> > L)
160 : {
161 : // cast to matrix based operator
162 0 : SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp =
163 : L.template cast_dynamic<MatrixOperator<matrix_type, vector_type> >();
164 :
165 : // Check that matrix if of correct type
166 0 : if(pOp.invalid())
167 0 : UG_THROW(name() << "::init': Passed Operator is "
168 : "not based on matrix. This Preconditioner can only "
169 : "handle matrix-based operators.");
170 :
171 0 : m_u = NULL;
172 :
173 : // forward request to matrix based implementation
174 0 : return base_type::init(pOp);
175 : }
176 :
177 : bool init(SmartPtr<MatrixOperator<matrix_type, vector_type> > Op)
178 : {
179 : m_u = NULL;
180 :
181 : return base_type::init(Op);
182 : }
183 :
184 : // Preprocess routine
185 0 : virtual bool preprocess(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp)
186 : {
187 0 : return preprocess_mat(*pOp);
188 : }
189 :
190 : public:
191 0 : virtual bool preprocess_mat(matrix_type &mat)
192 : {
193 : #ifdef UG_PARALLEL
194 : matrix_type m2;
195 : m2 = mat;
196 :
197 : MatAddSlaveRowsToMasterRowOverlap0(m2);
198 :
199 : // set zero on slaves
200 : std::vector<IndexLayout::Element> vIndex;
201 : CollectUniqueElements(vIndex, mat.layouts()->slave());
202 : SetDirichletRow(m2, vIndex);
203 :
204 : // Even after this setting of Dirichlet rows, it is possible that there are
205 : // zero rows on a proc because of the distribution:
206 : // For example, if one has a horizontal grid interface between two SHADOW_RIM_COPY
207 : // vertices, but the shadowing element for the hSlave side is vMaster (without being in
208 : // any horizontal interface). Then, as the horizontal interface (on the shadowed level)
209 : // is not part of the algebraic layouts, the hSlave is not converted into a Dirichlet row
210 : // by the previous commands.
211 : // As a first aid, we will simply convert any zero row on the current proc into a
212 : // Dirichlet row.
213 : // TODO: The corresponding rhs vector entry could be non-zero!
214 : // It is definitely not set to zero by change_storage_type(PST_UNIQUE) as the index is not contained
215 : // in the vector layouts either. Still, the defect assembling process might contain a vertex
216 : // loop and assemble something that is not solution-dependent! What do we do then?
217 : size_t nInd = m2.num_rows();
218 : size_t cnt = 0;
219 : for (size_t i = 0; i < nInd; ++i)
220 : {
221 : if (!m2.num_connections(i))
222 : {
223 : m2(i,i) = 1.0;
224 : ++cnt;
225 : }
226 : }
227 : #ifndef NDEBUG
228 : if (cnt) UG_LOG_ALL_PROCS("Converted "<<cnt<<" zero rows into Dirichlet rows.\n");
229 : #endif
230 :
231 : return preprocess_mat2(m2);
232 : #else
233 0 : return preprocess_mat2(mat);
234 : #endif
235 : }
236 :
237 0 : virtual bool preprocess_mat2(matrix_type &mat)
238 : {
239 : PROFILE_BEGIN_GROUP(ILUT_preprocess, "ilut algebra");
240 : //matrix_type &mat = *pOp;
241 : STATIC_ASSERT(matrix_type::rows_sorted, Matrix_has_to_have_sorted_rows);
242 0 : write_debug(mat, "ILUT_PreprocessIn");
243 :
244 : matrix_type* A;
245 0 : matrix_type permA;
246 :
247 0 : if(m_spOrderingAlgo.valid()){
248 0 : if(m_u){
249 0 : m_spOrderingAlgo->init(&mat, *m_u);
250 : }
251 : else{
252 0 : m_spOrderingAlgo->init(&mat);
253 : }
254 : }
255 :
256 0 : if(m_spOrderingAlgo.valid())
257 : {
258 0 : m_spOrderingAlgo->compute();
259 0 : m_ordering = m_spOrderingAlgo->ordering();
260 :
261 0 : m_bSortIsIdentity = GetInversePermutation(m_ordering, m_old_ordering);
262 :
263 0 : if(!m_bSortIsIdentity){
264 0 : SetMatrixAsPermutation(permA, mat, m_ordering);
265 : A = &permA;
266 : }
267 : else{
268 : A = &mat;
269 : }
270 : }
271 : else
272 : {
273 : A = &mat;
274 : }
275 :
276 0 : m_L.resize_and_clear(A->num_rows(), A->num_cols());
277 0 : m_U.resize_and_clear(A->num_rows(), A->num_cols());
278 :
279 : size_t totalentries=0;
280 : size_t maxentries=0;
281 :
282 0 : if((A->num_rows() > 0) && (A->num_cols() > 0)){
283 : // con is the current line of L/U
284 : // i also tried using std::list here or a special custom vector-based linked list
285 : // but vector is fastest, even with the insert operation.
286 : std::vector<matrix_connection> con;
287 0 : con.reserve(300);
288 0 : con.resize(0);
289 :
290 : // init row 0 of U
291 0 : for(matrix_row_iterator i_it = A->begin_row(0); i_it != A->end_row(0); ++i_it)
292 0 : con.push_back(matrix_connection(i_it.index(), i_it.value()));
293 0 : m_U.set_matrix_row(0, &con[0], con.size());
294 :
295 0 : Progress prog;
296 0 : if(m_show_progress)
297 0 : PROGRESS_START_WITH(prog, A->num_rows(),
298 : "Using ILUT(" << m_eps << ") on " << A->num_rows() << " x " << A->num_rows() << " matrix...");
299 :
300 0 : for(size_t i=1; i<A->num_rows(); i++)
301 : {
302 0 : if(m_show_progress) {PROGRESS_UPDATE(prog, i);}
303 0 : con.resize(0);
304 : size_t u_part=0;
305 :
306 0 : UG_COND_THROW(A->num_connections(i) == 0, "row " << i << " has no connections");
307 :
308 : // get the row A(i, .) into con
309 : double dmax=0;
310 0 : for(matrix_row_iterator i_it = A->begin_row(i); i_it != A->end_row(i); ++i_it)
311 : {
312 0 : con.push_back(matrix_connection(i_it.index(), i_it.value()));
313 0 : if(dmax < BlockNorm(i_it.value()))
314 0 : dmax = BlockNorm(i_it.value());
315 : }
316 :
317 : // eliminate all entries A(i, k) with k<i with rows U(k, .) and k<i
318 0 : for(size_t i_it = 0; i_it < con.size(); ++i_it)
319 : {
320 0 : size_t k = con[i_it].iIndex;
321 0 : if(k >= i)
322 : {
323 : // safe where U begins / L ends in con
324 : u_part = i_it;
325 0 : break;
326 : }
327 0 : if(con[i_it].dValue == 0.0) continue;
328 0 : UG_COND_THROW(!(m_U.num_connections(k) != 0 && m_U.begin_row(k).index() == k), "");
329 : block_type &ukk = m_U.begin_row(k).value();
330 :
331 : // add row k to row i by A(i, .) -= U(k,.) A(i,k) / U(k,k)
332 : // so that A(i,k) is zero.
333 : // safe A(i,k)/U(k,k) in con, (later L(i,k) )
334 0 : con[i_it].dValue = con[i_it].dValue / ukk;
335 : block_type d = con[i_it].dValue;
336 0 : UG_COND_THROW(!BlockMatrixFiniteAndNotTooBig(d, 1e40), "i = " << i << " " << d);
337 :
338 : matrix_row_iterator k_it = m_U.begin_row(k); // upper row iterator
339 : ++k_it; // skip diag
340 0 : size_t j = i_it+1;
341 0 : while(k_it != m_U.end_row(k) && j < con.size())
342 : {
343 : // (since con and U[k] is sorted, we can do sth like a merge on the two lists)
344 0 : if(k_it.index() == con[j].iIndex)
345 : {
346 : // match
347 0 : con[j].dValue -= k_it.value() * d;
348 0 : ++k_it; ++j;
349 : }
350 0 : else if(k_it.index() < con[j].iIndex)
351 : {
352 : // we have a value in U(k, (*k_it).iIndex), but not in A.
353 : // check tolerance criteria
354 :
355 0 : matrix_connection c(k_it.index(), k_it.value() * d * -1.0);
356 :
357 0 : UG_COND_THROW(!BlockMatrixFiniteAndNotTooBig(c.dValue, 1e40), "i = " << i << " " << c.dValue);
358 0 : if(BlockNorm(c.dValue) > dmax * m_eps)
359 : {
360 : // insert sorted
361 0 : con.insert(con.begin()+j, c);
362 0 : ++j; // don't do this when using iterators
363 : }
364 : // else do some lumping
365 : ++k_it;
366 : }
367 : else
368 : {
369 : // we have a value in A(k, con[j].iIndex), but not in U.
370 0 : ++j;
371 : }
372 : }
373 : // insert new connections after last connection of row i
374 0 : if (k_it!=m_U.end_row(k)){
375 0 : for (;k_it!=m_U.end_row(k);++k_it){
376 0 : matrix_connection c(k_it.index(),-k_it.value()*d);
377 0 : if(BlockNorm(c.dValue) > dmax * m_eps)
378 : {
379 0 : con.push_back(c);
380 : }
381 : }
382 : };
383 : }
384 :
385 :
386 0 : totalentries+=con.size();
387 : if(maxentries < con.size()) maxentries = con.size();
388 :
389 : // safe L and U
390 0 : m_L.set_matrix_row(i, &con[0], u_part);
391 0 : m_U.set_matrix_row(i, &con[u_part], con.size()-u_part);
392 :
393 : }
394 0 : if(m_show_progress) {PROGRESS_FINISH(prog);}
395 :
396 0 : m_L.defragment();
397 0 : m_U.defragment();
398 0 : }
399 :
400 0 : if (m_info==true)
401 : {
402 0 : m_L.print("L");
403 0 : m_U.print("U");
404 : UG_LOG("\n ILUT storage information:\n");
405 : UG_LOG(" A is " << A->num_rows() << " x " << A->num_cols() << " matrix.\n");
406 : UG_LOG(" A nr of connections: " << A->total_num_connections() << "\n");
407 0 : UG_LOG(" L+U nr of connections: " << m_L.total_num_connections()+m_U.total_num_connections() << "\n");
408 0 : UG_LOG(" Increase factor: " << (float)(m_L.total_num_connections() + m_U.total_num_connections() )/A->total_num_connections() << "\n");
409 0 : UG_LOG(reset_floats << " Total entries: " << totalentries << " (" << ((double)totalentries) / (A->num_rows()*A->num_rows()) << "% of dense)\n");
410 0 : if(m_spOrderingAlgo.valid())
411 : {
412 0 : UG_LOG(" Using " << m_spOrderingAlgo->name() << "\n");
413 : }
414 : else
415 : {
416 : UG_LOG("Not using sort.");
417 : }
418 : }
419 :
420 0 : return true;
421 0 : }
422 :
423 : // Stepping routine
424 0 : virtual bool step(SmartPtr<MatrixOperator<matrix_type, vector_type> > pOp, vector_type& c, const vector_type& d)
425 : {
426 : #ifdef UG_PARALLEL
427 : SmartPtr<vector_type> spDtmp = d.clone();
428 : spDtmp->change_storage_type(PST_UNIQUE);
429 : bool b = solve(c, *spDtmp);
430 :
431 : c.set_storage_type(PST_ADDITIVE);
432 : c.change_storage_type(PST_CONSISTENT);
433 : return b;
434 : #else
435 0 : return solve(c, d);
436 : #endif
437 : return true;
438 : }
439 :
440 0 : virtual bool solve(vector_type& c, const vector_type& d)
441 : {
442 0 : if(m_spOrderingAlgo.invalid()|| m_bSortIsIdentity)
443 : {
444 0 : if(applyLU(c, d) == false) return false;
445 : }
446 : else
447 : {
448 : PROFILE_BEGIN_GROUP(ILUT_StepWithReorder, "ilut algebra");
449 :
450 0 : SetVectorAsPermutation(c, d, m_ordering);
451 0 : c2.resize(c.size());
452 0 : if(applyLU(c2, c) == false) return false;
453 0 : SetVectorAsPermutation(c, c2, m_old_ordering);
454 : }
455 : return true;
456 : }
457 :
458 :
459 0 : virtual bool applyLU(vector_type& c, const vector_type& d)
460 : {
461 : PROFILE_BEGIN_GROUP(ILUT_step, "ilut algebra");
462 : // apply iterator: c = LU^{-1}*d (damp is not used)
463 : // L
464 0 : for(size_t i=0; i < m_L.num_rows(); i++)
465 : {
466 : // c[i] = d[i] - m_L[i]*c;
467 0 : c[i] = d[i];
468 0 : for(matrix_row_iterator it = m_L.begin_row(i); it != m_L.end_row(i); ++it)
469 0 : MatMultAdd(c[i], 1.0, c[i], -1.0, it.value(), c[it.index()] );
470 : // lii = 1.0.
471 : }
472 :
473 : // U
474 : //
475 : // last row diagonal U entry might be close to zero with corresponding zero rhs
476 : // when solving Navier Stokes system, therefore handle separately
477 0 : if(m_U.num_rows() > 0)
478 : {
479 0 : size_t i=m_U.num_rows()-1;
480 0 : matrix_row_iterator it = m_U.begin_row(i);
481 : UG_ASSERT(it != m_U.end_row(i), i);
482 : UG_ASSERT(it.index() == i, i);
483 : block_type &uii = it.value();
484 0 : vector_value s = c[i];
485 : // check if diag part is significantly smaller than rhs
486 : // This may happen when matrix is indefinite with one eigenvalue
487 : // zero. In that case, the factorization on the last row is
488 : // nearly zero due to round-off errors. In order to allow ill-
489 : // scaled matrices (i.e. small matrix entries row-wise) this
490 : // is compared to the rhs, that is small in this case as well.
491 : if (false && BlockNorm(uii) <= m_small * m_eps * BlockNorm(s)){
492 : UG_LOG("ILUT("<<m_eps<<") Warning: Near-zero diagonal entry "
493 : "with norm "<<BlockNorm(uii)<<" in last row of U "
494 : " with corresponding non-near-zero rhs with norm "
495 : << BlockNorm(s) << ". Setting rhs to zero.\n");
496 : // set correction to zero
497 : c[i] = 0;
498 : } else {
499 : // c[i] = s/uii;
500 : InverseMatMult(c[i], 1.0, uii, s);
501 : }
502 : }
503 :
504 : // handle all other rows
505 0 : if(m_U.num_rows() > 1){
506 0 : for(size_t i=m_U.num_rows()-2; ; i--)
507 : {
508 0 : matrix_row_iterator it = m_U.begin_row(i);
509 : UG_ASSERT(it != m_U.end_row(i), i);
510 : UG_ASSERT(it.index() == i, i);
511 : block_type &uii = it.value();
512 :
513 0 : vector_value s = c[i];
514 : ++it; // skip diag
515 0 : for(; it != m_U.end_row(i); ++it){
516 : // s -= it.value() * c[it.index()];
517 0 : MatMultAdd(s, 1.0, s, -1.0, it.value(), c[it.index()] );
518 : }
519 : // c[i] = s/uii;
520 : InverseMatMult(c[i], 1.0, uii, s);
521 :
522 0 : if(i==0) break;
523 : }
524 : }
525 0 : return true;
526 : }
527 :
528 0 : virtual bool multi_apply(std::vector<vector_type> &vc, const std::vector<vector_type> &vd)
529 : {
530 : PROFILE_BEGIN_GROUP(ILUT_step, "ilut algebra");
531 : // apply iterator: c = LU^{-1}*d (damp is not used)
532 : // L
533 0 : for(size_t i=0; i < m_L.num_rows(); i++)
534 : {
535 :
536 0 : for(size_t e=0; e<vc.size(); e++)
537 : {
538 : vector_type &c = vc[e];
539 : const vector_type &d = vd[e];
540 : // c[i] = d[i] - m_L[i]*c;
541 0 : c[i] = d[i];
542 0 : for(matrix_row_iterator it = m_L.begin_row(i); it != m_L.end_row(i); ++it)
543 0 : MatMultAdd(c[i], 1.0, c[i], -1.0, it.value(), c[it.index()] );
544 : // lii = 1.0.
545 : }
546 : }
547 :
548 : // U
549 : //
550 : // last row diagonal U entry might be close to zero with corresponding zero rhs
551 : // when solving Navier Stokes system, therefore handle separately
552 0 : if(m_U.num_rows() > 0)
553 : {
554 0 : for(size_t e=0; e<vc.size(); e++)
555 : {
556 : vector_type &c = vc[e];
557 0 : size_t i=m_U.num_rows()-1;
558 0 : matrix_row_iterator it = m_U.begin_row(i);
559 : UG_ASSERT(it != m_U.end_row(i), i);
560 : UG_ASSERT(it.index() == i, i);
561 : block_type &uii = it.value();
562 0 : vector_value s = c[i];
563 : // check if diag part is significantly smaller than rhs
564 : // This may happen when matrix is indefinite with one eigenvalue
565 : // zero. In that case, the factorization on the last row is
566 : // nearly zero due to round-off errors. In order to allow ill-
567 : // scaled matrices (i.e. small matrix entries row-wise) this
568 : // is compared to the rhs, that is small in this case as well.
569 : if (false && BlockNorm(uii) <= m_small * m_eps * BlockNorm(s)){
570 : UG_LOG("ILUT("<<m_eps<<") Warning: Near-zero diagonal entry "
571 : "with norm "<<BlockNorm(uii)<<" in last row of U "
572 : " with corresponding non-near-zero rhs with norm "
573 : << BlockNorm(s) << ". Setting rhs to zero.\n");
574 : // set correction to zero
575 : c[i] = 0;
576 : } else {
577 : // c[i] = s/uii;
578 : InverseMatMult(c[i], 1.0, uii, s);
579 : }
580 : }
581 : }
582 :
583 : // handle all other rows
584 0 : if(m_U.num_rows() > 1){
585 0 : for(size_t i=m_U.num_rows()-2; ; i--)
586 : {
587 0 : for(size_t e=0; e<vc.size(); e++)
588 : {
589 : vector_type &c = vc[e];
590 0 : matrix_row_iterator it = m_U.begin_row(i);
591 : UG_ASSERT(it != m_U.end_row(i), i);
592 : UG_ASSERT(it.index() == i, i);
593 : block_type &uii = it.value();
594 :
595 0 : vector_value s = c[i];
596 : ++it; // skip diag
597 0 : for(; it != m_U.end_row(i); ++it){
598 : // s -= it.value() * c[it.index()];
599 0 : MatMultAdd(s, 1.0, s, -1.0, it.value(), c[it.index()] );
600 : }
601 : // c[i] = s/uii;
602 : InverseMatMult(c[i], 1.0, uii, s);
603 :
604 0 : if(i==0) break;
605 : }
606 : }
607 : }
608 0 : return true;
609 : }
610 :
611 : // Postprocess routine
612 0 : virtual bool postprocess() {return true;}
613 :
614 : protected:
615 : vector_type c2;
616 : matrix_type m_L;
617 : matrix_type m_U;
618 : double m_eps;
619 : bool m_info;
620 : bool m_show_progress;
621 : static const number m_small;
622 : std::vector<size_t> newIndex, oldIndex;
623 :
624 : /// for ordering algorithms
625 : SmartPtr<ordering_algo_type> m_spOrderingAlgo;
626 : ordering_container_type m_ordering, m_old_ordering;
627 :
628 : bool m_bSortIsIdentity;
629 :
630 : const vector_type* m_u;
631 : };
632 :
633 : // define constant
634 : template <typename TAlgebra>
635 : const number ILUTPreconditioner<TAlgebra>::m_small = 1e-8;
636 :
637 : } // end namespace ug
638 :
639 : #endif
|
