general_purpose_block_preconditioners.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26
27#include "oomph_definitions.h"
28#include "preconditioner.h"
30#include "matrices.h"
32
33namespace oomph
34{
35 //============================================================================
36 /// setup for the block diagonal preconditioner
37 //============================================================================
38 template<typename MATRIX>
40 {
41 // clean the memory
42 this->clean_up_memory();
43
44 // Note: Subsidiary block preconditioners don't really use their
45 // mesh pointers (since the lookup schemes inferred from them are
46 // set up by their masters). Generally we insist that (for uniformity)
47 // a mesh should be set, but sometimes subsidiary preconditioners are
48 // set up on the fly and we can't sensibly set meshes, so we're
49 // forgiving...)
50 if (this->is_master_block_preconditioner())
51 {
52#ifdef PARANOID
53 if (this->gp_nmesh() == 0)
54 {
55 std::ostringstream err_msg;
56 err_msg << "There are no meshes set.\n"
57 << "Did you remember to call add_mesh(...)?";
58 throw OomphLibError(
59 err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
60 }
61#endif
62
63 // Set all meshes if this is master block preconditioner
64 this->gp_preconditioner_set_all_meshes();
65 }
66
67 // If we're meant to build silently
68 if (this->Silent_preconditioner_setup == true)
69 {
70 // Store the output stream pointer
71 this->Stream_pt = oomph_info.stream_pt();
72
73 // Now set the oomph_info stream pointer to the null stream to
74 // disable all possible output
76 }
77
78 // Set up the block look up schemes
79 this->gp_preconditioner_block_setup();
80
81 // number of types of blocks.
82 unsigned nblock_types = this->nblock_types();
83
84 // Create any subsidiary preconditioners needed
85 this->fill_in_subsidiary_preconditioners(nblock_types);
86
87 // The total time for extracting all the blocks from the "global" matrix
88 double t_extraction_total = 0.0;
89
90 // The total time for setting up the subsidiary preconditioners
91 double t_subsidiary_setup_total = 0.0;
92
93 // If using two level parallelisation then we need to use a
94 // PrecondtionerArray which requires very different setup. ??ds possibly
95 // it should have it's own class?
96 if (Use_two_level_parallelisation)
97 {
98 // Get the blocks. We have to use new because you can't have containers
99 // of matrices because the copy constructor is buggy (so we create a
100 // container of pointers instead). ??ds
101 Vector<CRDoubleMatrix*> block_diagonal_matrix_pt(nblock_types, 0);
102 for (unsigned i = 0; i < nblock_types; i++)
103 {
104 // Allocate space for the new matrix
105 block_diagonal_matrix_pt[i] = new CRDoubleMatrix;
106
107 // Get the start time
108 double t_extract_start = TimingHelpers::timer();
109
110 // Extract the i-th block
111 this->get_block(
112 i, get_other_diag_ds(i, nblock_types), *block_diagonal_matrix_pt[i]);
113
114 // Get the end time
115 double t_extract_end = TimingHelpers::timer();
116
117 // Update the timing total
118 t_extraction_total += (t_extract_end - t_extract_start);
119 }
120
121 // Construct the preconditioner array
122 Preconditioner_array_pt = new PreconditionerArray;
123
124 // Get the start time
125 double t_subsidiary_setup_start = TimingHelpers::timer();
126
127 // Set up the preconditioners
128 Preconditioner_array_pt->setup_preconditioners(
129 block_diagonal_matrix_pt,
130 this->Subsidiary_preconditioner_pt,
131 this->comm_pt());
132
133 // Get the end time
134 double t_subsidiary_setup_end = TimingHelpers::timer();
135
136 // Update the timing total
137 t_subsidiary_setup_total +=
138 (t_subsidiary_setup_end - t_subsidiary_setup_start);
139
140 // and delete the blocks
141 for (unsigned i = 0; i < nblock_types; i++)
142 {
143 delete block_diagonal_matrix_pt[i];
144 block_diagonal_matrix_pt[i] = 0;
145 }
146
147 // Preconditioner array is weird, it calls delete on all the
148 // preconditioners you give it and requires new ones each time!
149 this->Subsidiary_preconditioner_pt.clear();
150 }
151 // Otherwise just set up each block's preconditioner in order
152 else
153 {
154 for (unsigned i = 0; i < nblock_types; i++)
155 {
156 // Get the block
157 unsigned j = get_other_diag_ds(i, nblock_types);
158
159 // Get the start time
160 double t_extract_start = TimingHelpers::timer();
161
162 // Get the block
163 CRDoubleMatrix block = this->get_block(i, j);
164
165 // Get the end time
166 double t_extract_end = TimingHelpers::timer();
167
168 // Update the timing total
169 t_extraction_total += (t_extract_end - t_extract_start);
170
171 // Get the start time
172 double t_subsidiary_setup_start = TimingHelpers::timer();
173
174 // Set the (subsidiary) preconditioner up for this block
175 this->Subsidiary_preconditioner_pt[i]->setup(&block);
176
177 // Get the end time
178 double t_subsidiary_setup_end = TimingHelpers::timer();
179
180 // Tell the user
181 oomph_info << "Took "
182 << t_subsidiary_setup_end - t_subsidiary_setup_start
183 << "s to setup." << std::endl;
184
185 // Update the timing total
186 t_subsidiary_setup_total +=
187 (t_subsidiary_setup_end - t_subsidiary_setup_start);
188 }
189 }
190
191 // If we're meant to build silently, reassign the oomph stream pointer
192 if (this->Silent_preconditioner_setup == true)
193 {
194 // Store the output stream pointer
195 oomph_info.stream_pt() = this->Stream_pt;
196
197 // Reset our own stream pointer
198 this->Stream_pt = 0;
199 }
200
201 // Tell the user
202 oomph_info << "Total block extraction time [sec]: " << t_extraction_total
203 << "\nTotal subsidiary preconditioner setup time [sec]: "
204 << t_subsidiary_setup_total << std::endl;
205 }
206
207 //=============================================================================
208 /// Preconditioner solve for the block diagonal preconditioner
209 //=============================================================================
210 template<typename MATRIX>
212 const DoubleVector& r, DoubleVector& z)
213 {
214 // Cache umber of block types
215 unsigned n_block = this->nblock_types();
216
217 // Get the right hand side vector (b in Ax = b) in block form.
218 Vector<DoubleVector> block_r;
219 this->get_block_vectors(r, block_r);
220
221 // Make sure z vector is built
222 if (!z.built())
223 {
224 z.build(this->distribution_pt(), 0.0);
225 }
226
227 // Vector of vectors for storage of solution in block form.
228 Vector<DoubleVector> block_z(n_block);
229
230 if (Use_two_level_parallelisation)
231 {
232 Preconditioner_array_pt->solve_preconditioners(block_r, block_z);
233 }
234 else
235 {
236 // solve each diagonal block
237 for (unsigned i = 0; i < n_block; i++)
238 {
239 double t_start = 0.0;
240 if (Doc_time_during_preconditioner_solve)
241 {
242 t_start = TimingHelpers::timer();
243 }
244
245 // this->Subsidiary_preconditioner_pt[i]->preconditioner_solve(block_r[i],
246 // block_z[i]);
247
248
249 // See if the i-th subsidiary preconditioner is a block preconditioner
250 if (dynamic_cast<BlockPreconditioner<CRDoubleMatrix>*>(
251 this->Subsidiary_preconditioner_pt[i]) == 0)
252 {
253 this->Subsidiary_preconditioner_pt[i]->preconditioner_solve(
254 block_r[i], block_z[i]);
255 }
256 // If the subsidiary preconditioner is a block preconditioner
257 else
258 {
259 // This is pretty inefficient but there is no other way to do this
260 // with block sub preconditioners as far as I can tell because they
261 // demand to have the full r and z vectors...
262 DoubleVector block_z_with_size_of_full_z(r.distribution_pt());
263 DoubleVector r_updated(r.distribution_pt());
264
265 // Construct the new r vector (with the update given by back subs
266 // below).
267 this->return_block_vectors(block_r, r_updated);
268
269 // Solve (blocking is handled inside the block preconditioner).
270 this->Subsidiary_preconditioner_pt[i]->preconditioner_solve(
271 r_updated, block_z_with_size_of_full_z);
272
273 // Extract this block's z vector because we need it for the next step
274 // (and throw away block_z_with_size_of_full_z).
275 this->get_block_vector(i, block_z_with_size_of_full_z, block_z[i]);
276 }
277
278
279 if (Doc_time_during_preconditioner_solve)
280 {
281 oomph_info << "Time for application of " << i
282 << "-th block preconditioner: "
283 << TimingHelpers::timer() - t_start << std::endl;
284 }
285 }
286 }
287
288 // copy solution in block vectors block_z back to z
289 this->return_block_vectors(block_z, z);
290 }
291
292
293 //============================================================================
294 /// setup for the block triangular preconditioner
295 //============================================================================
296 template<typename MATRIX>
298 {
299 // clean the memory
300 this->clean_up_memory();
301
302 // Subsidiary preconditioners don't really need the meshes
303 if (this->is_master_block_preconditioner())
304 {
305#ifdef PARANOID
306 if (this->gp_nmesh() == 0)
307 {
308 std::ostringstream err_msg;
309 err_msg << "There are no meshes set.\n"
310 << "Did you remember to call add_mesh(...)?";
311 throw OomphLibError(
312 err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
313 }
314#endif
315
316 // Set all meshes if this is master block preconditioner
317 this->gp_preconditioner_set_all_meshes();
318 }
319
320 // If we're meant to build silently
321 if (this->Silent_preconditioner_setup == true)
322 {
323 // Store the output stream pointer
324 this->Stream_pt = oomph_info.stream_pt();
325
326 // Now set the oomph_info stream pointer to the null stream to
327 // disable all possible output
329 }
330
331 // Set up the block look up schemes
332 this->gp_preconditioner_block_setup();
333
334 // number of block types
335 unsigned nblock_types = this->nblock_types();
336
337 // storage for the off diagonal matrix vector products
338 Off_diagonal_matrix_vector_products.resize(nblock_types, nblock_types, 0);
339
340 // Fill in any null subsidiary preconditioners
341 this->fill_in_subsidiary_preconditioners(nblock_types);
342
343 // The total time for extracting all the blocks from the "global" matrix
344 double t_extraction_total = 0.0;
345
346 // The total time for setting up the subsidiary preconditioners
347 double t_subsidiary_setup_total = 0.0;
348
349 // The total time for setting up the matrix-vector products
350 double t_mvp_setup_total = 0.0;
351
352 // build the preconditioners and matrix vector products
353 for (unsigned i = 0; i < nblock_types; i++)
354 {
355 // Get the block and set up the preconditioner.
356 {
357 // Get the start time
358 double t_extract_start = TimingHelpers::timer();
359
360 // Grab the i-th diagonal block
361 CRDoubleMatrix block_matrix = this->get_block(i, i);
362
363 // Get the end time
364 double t_extract_end = TimingHelpers::timer();
365
366 // Update the timing total
367 t_extraction_total += (t_extract_end - t_extract_start);
368
369 // Get the start time
370 double t_subsidiary_setup_start = TimingHelpers::timer();
371
372 // Set up the i-th subsidiary preconditioner with this block
373 this->Subsidiary_preconditioner_pt[i]->setup(&block_matrix);
374
375 // Get the end time
376 double t_subsidiary_setup_end = TimingHelpers::timer();
377
378 // Update the timing total
379 t_subsidiary_setup_total +=
380 (t_subsidiary_setup_end - t_subsidiary_setup_start);
381 }
382
383 // next setup the off diagonal mat vec operators
384 unsigned l, u;
385 if (Upper_triangular)
386 {
387 l = i + 1;
388 u = nblock_types;
389 }
390 else
391 {
392 l = 0;
393 u = i;
394 }
395
396 for (unsigned j = l; j < u; j++)
397 {
398 // Get the start time
399 double t_extract_start = TimingHelpers::timer();
400
401 // Get the block
402 CRDoubleMatrix block_matrix = this->get_block(i, j);
403
404 // Get the end time
405 double t_extract_end = TimingHelpers::timer();
406
407 // Update the timing total
408 t_extraction_total += (t_extract_end - t_extract_start);
409
410 // Get the start time
411 double t_mvp_start = TimingHelpers::timer();
412
413 // Copy the block into a "multiplier" class. If trilinos is being
414 // used this should also be faster than oomph-lib's multiphys.
415 Off_diagonal_matrix_vector_products(i, j) = new MatrixVectorProduct();
416
417 // Set the damn thing up
418 this->setup_matrix_vector_product(
419 Off_diagonal_matrix_vector_products(i, j), &block_matrix, j);
420
421 // Get the end time
422 double t_mvp_end = TimingHelpers::timer();
423
424 // Update the timing total
425 t_mvp_setup_total += (t_mvp_end - t_mvp_start);
426 }
427 }
428
429 // Tell the user
430 oomph_info << "Total block extraction time [sec]: " << t_extraction_total
431 << "\nTotal subsidiary preconditioner setup time [sec]: "
432 << t_subsidiary_setup_total
433 << "\nTotal matrix-vector product setup time [sec]: "
434 << t_mvp_setup_total << std::endl;
435
436 // If we're meant to build silently, reassign the oomph stream pointer
437 if (this->Silent_preconditioner_setup == true)
438 {
439 // Store the output stream pointer
440 oomph_info.stream_pt() = this->Stream_pt;
441
442 // Reset our own stream pointer
443 this->Stream_pt = 0;
444 }
445 } // End of setup
446
447 //=============================================================================
448 /// Preconditioner solve for the block triangular preconditioner
449 //=============================================================================
450 template<typename MATRIX>
452 const DoubleVector& r, DoubleVector& z)
453 {
454 // Cache number of block types
455 unsigned n_block = this->nblock_types();
456
457 //
458 int start, end, step;
459 if (Upper_triangular)
460 {
461 start = n_block - 1;
462 end = -1;
463 step = -1;
464 }
465 else
466 {
467 start = 0;
468 end = n_block;
469 step = 1;
470 }
471
472 // vector of vectors for each section of residual vector
473 Vector<DoubleVector> block_r;
474
475 // rearrange the vector r into the vector of block vectors block_r
476 this->get_block_vectors(r, block_r);
477
478 // vector of vectors for the solution block vectors
479 Vector<DoubleVector> block_z(n_block);
480
481 //
482 for (int i = start; i != end; i += step)
483 {
484 // ??ds ugly, fix this?
485 if (dynamic_cast<BlockPreconditioner<CRDoubleMatrix>*>(
486 this->Subsidiary_preconditioner_pt[i]) == 0)
487 {
488 // solve on the block
489 this->Subsidiary_preconditioner_pt[i]->preconditioner_solve(block_r[i],
490 block_z[i]);
491 }
492 else
493 {
494 // This is pretty inefficient but there is no other way to do this with
495 // block sub preconditioners as far as I can tell because they demand
496 // to have the full r and z vectors...
497
498 DoubleVector block_z_with_size_of_full_z(r.distribution_pt()),
499 r_updated(r.distribution_pt());
500
501 // Construct the new r vector (with the update given by back subs
502 // below).
503 this->return_block_vectors(block_r, r_updated);
504
505 // Solve (blocking is handled inside the block preconditioner).
506 this->Subsidiary_preconditioner_pt[i]->preconditioner_solve(
507 r_updated, block_z_with_size_of_full_z);
508
509 // Extract this block's z vector because we need it for the next step
510 // (and throw away block_z_with_size_of_full_z).
511 this->get_block_vector(i, block_z_with_size_of_full_z, block_z[i]);
512 }
513
514 // substitute
515 for (int j = i + step; j != end; j += step)
516 {
517 DoubleVector temp;
518 Off_diagonal_matrix_vector_products(j, i)->multiply(block_z[i], temp);
519 block_r[j] -= temp;
520 }
521 }
522
523 // copy solution in block vectors block_r back to z
524 this->return_block_vectors(block_z, z);
525 }
526
527
528 //=============================================================================
529 /// Setup for the exact block preconditioner
530 //=============================================================================
531 template<typename MATRIX>
533 {
534 // If this is a master block preconditioner,
535 // we give the mesh pointers to the BlockPreconditioner base class.
536 // Subsidiary preconditioners don't need meshes...
537 if (this->is_master_block_preconditioner())
538 {
539#ifdef PARANOID
540 if (this->gp_nmesh() == 0)
541 {
542 std::ostringstream err_msg;
543 err_msg << "There are no meshes set.\n"
544 << "Did you remember to call add_mesh(...)?";
545 throw OomphLibError(
546 err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
547 }
548#endif
549
550 // Set all meshes if this is master block preconditioner
551 this->gp_preconditioner_set_all_meshes();
552 }
553
554 // Set up the block look up schemes
555 this->gp_preconditioner_block_setup();
556
557 // get the number of DOF types
558 unsigned nblock_types = this->nblock_types();
559
560 // Build the preconditioner matrix
561 VectorMatrix<BlockSelector> required_blocks(nblock_types, nblock_types);
562
563 // boolean indicating if we want the block, stored for readability.
564 const bool want_block = true;
565 for (unsigned b_i = 0; b_i < nblock_types; b_i++)
566 {
567 for (unsigned b_j = 0; b_j < nblock_types; b_j++)
568 {
569 required_blocks[b_i][b_j].select_block(b_i, b_j, want_block);
570 }
571 }
572
573 // Get the concatenated blocks
574 CRDoubleMatrix exact_block_matrix =
575 this->get_concatenated_block(required_blocks);
576
577 // Setup the preconditioner.
578 this->fill_in_subsidiary_preconditioners(1); // Only need one preconditioner
579 preconditioner_pt()->setup(&exact_block_matrix);
580 }
581
582 //=============================================================================
583 /// Preconditioner solve for the exact block preconditioner
584 //=============================================================================
585 template<typename MATRIX>
587 const DoubleVector& r, DoubleVector& z)
588 {
589 // get the block ordered components of the r vector for this preconditioner
590 DoubleVector block_order_r;
591 this->get_block_ordered_preconditioner_vector(r, block_order_r);
592
593 // vector for solution
594 DoubleVector block_order_z;
595
596 // apply the preconditioner
597 preconditioner_pt()->preconditioner_solve(block_order_r, block_order_z);
598
599 // copy solution back to z vector
600 this->return_block_ordered_preconditioner_vector(block_order_z, z);
601 }
602
608
609} // namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
Block "anti-diagonal" preconditioner, i.e. same as block diagonal but along the other diagonal of the...
Block diagonal preconditioner. By default SuperLU is used to solve the subsidiary systems,...
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to r.
virtual void setup()
Setup the preconditioner.
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to r.
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:888
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
bool built() const
Preconditioner that doesn't actually do any preconditioning, it just allows access to the Jacobian bl...
//////////////////////////////////////////////////////////////////////////// ////////////////////////...
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to r.
Matrix vector product helper class - primarily a wrapper to Trilinos's Epetra matrix vector product m...
std::ostream *& stream_pt()
Access function for the stream pointer.
An OomphLibError object which should be thrown when an run-time error is encountered....
PreconditionerArray - NOTE - first implementation, a number of assumptions / simplifications were mad...
VectorMatrix is a generalised, STL-map-based, matrix based on a Vector of Vectors.
Definition: vector_matrix.h:79
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
void start(const unsigned &i)
(Re-)start i-th timer
void clean_up_memory()
Clean up function that deletes anything dynamically allocated in this namespace.
double timer()
returns the time in seconds after some point in past
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Nullstream oomph_nullstream
///////////////////////////////////////////////////////////////////////// ///////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...