general_purpose_space_time_block_preconditioner.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// Oomph-lib headers
28
29// Header file
31
32/// /////////////////////////////////////////////////////////////////////////
33/// /////////////////////////////////////////////////////////////////////////
34/// /////////////////////////////////////////////////////////////////////////
35
36namespace oomph
37{
38 //============================================================================
39 /// Set up the exact preconditioner
40 //============================================================================
41 template<typename MATRIX>
43 {
44 // Clean the memory
45 this->clean_up_memory();
46
47 // Subsidiary preconditioners don't really need the meshes
48 if (this->is_master_block_preconditioner())
49 {
50#ifdef PARANOID
51 if (this->gp_nmesh() == 0)
52 {
53 std::ostringstream err_msg;
54 err_msg << "There are no meshes set.\n"
55 << "Did you remember to call add_mesh(...)?";
56 throw OomphLibError(
57 err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
58 }
59#endif
60
61 // Set all meshes if this is master block preconditioner
62 this->gp_preconditioner_set_all_meshes();
63 }
64
65 // Initialise the memory usage variable
66 Memory_usage_in_bytes = 0.0;
67
68 // If we're meant to build silently
69 if (this->Silent_preconditioner_setup == true)
70 {
71 // Store the output stream pointer
72 this->Stream_pt = oomph_info.stream_pt();
73
74 // Now set the oomph_info stream pointer to the null stream to
75 // disable all possible output
77 }
78
79 // Set up the block look up schemes
80 this->gp_preconditioner_block_setup();
81
82 // Number of block types
83 unsigned n_block_types = this->nblock_types();
84
85 // Fill in any null subsidiary preconditioners
86 this->fill_in_subsidiary_preconditioners(n_block_types);
87
88 // The total time for setting up the subsidiary preconditioners
89 double t_subsidiary_setup_total = 0.0;
90
91 // Stores the blocks we want to extract
92 VectorMatrix<BlockSelector> required_blocks(n_block_types, n_block_types);
93
94 // Boolean indicating if we want the block, stored for readability.
95 const bool want_block = true;
96
97 // Loop over the block rows
98 for (unsigned i = 0; i < n_block_types; i++)
99 {
100 // Loop over the block columns
101 for (unsigned j = 0; j < n_block_types; j++)
102 {
103 // Indicate that we want this block
104 required_blocks[i][j].select_block(i, j, want_block);
105 }
106 } // for (unsigned i=0;i<n_block_types;i++)
107
108 // Get the start time
109 double t_extract_start = TimingHelpers::timer();
110
111 // Get the concatenated blocks
112 CRDoubleMatrix full_matrix = this->get_concatenated_block(required_blocks);
113
114 // The total time for extracting all the blocks from the "global" matrix
115 double t_extraction_total = (TimingHelpers::timer() - t_extract_start);
116
117 // Get the start time
118 double t_subsidiary_setup_start = TimingHelpers::timer();
119
120 // It's a block preconditioner so pass it a pointer to the global matrix!
121 this->Subsidiary_preconditioner_pt[0]->setup(&full_matrix);
122
123 // Update the timing total
124 t_subsidiary_setup_total +=
125 (TimingHelpers::timer() - t_subsidiary_setup_start);
126
127 // Remember that the preconditioner has been set up
128 Preconditioner_has_been_setup = true;
129
130 // If we're meant to build silently, reassign the oomph stream pointer
131 if (this->Silent_preconditioner_setup == true)
132 {
133 // Store the output stream pointer
134 oomph_info.stream_pt() = this->Stream_pt;
135
136 // Reset our own stream pointer
137 this->Stream_pt = 0;
138 }
139
140 // Tell the user
141 oomph_info << "Total block extraction time [sec]: " << t_extraction_total
142 << "\nTotal subsidiary preconditioner setup time [sec]: "
143 << t_subsidiary_setup_total << std::endl;
144
145 // Do we need to doc. the memory statistics?
146 // NOTE: We're going to assume that either:
147 // (1) GMRES is used as a subsidiary block preconditioner (preconditioned
148 // by the Navier-Stokes subsidiary block preconditioner), or;
149 // (2) SuperLU is used as the subsidiary preconditioner.
150 if (Compute_memory_statistics)
151 {
152 // How many rows are there in the global Jacobian?
153 unsigned n_row = this->matrix_pt()->nrow();
154
155 // How many nonzeros are there in the global Jacobian?
156 unsigned n_nnz = this->matrix_pt()->nnz();
157
158 // Add in the subsidiary preconditioners contribution
159 double total_memory_usage_for_setup_phase =
160 dynamic_cast<SuperLUPreconditioner*>(
161 this->Subsidiary_preconditioner_pt[0])
162 ->get_total_memory_needed_for_superlu();
163
164 // Add in the global Jacobian contribution
165 total_memory_usage_for_setup_phase +=
166 ((2 * ((n_row + 1) * sizeof(int))) +
167 (n_nnz * (sizeof(int) + sizeof(double))));
168
169 // How much memory have we used in the subsidiary preconditioners?
170 oomph_info << "\nTotal amount of memory being used after setup (MB): "
171 << total_memory_usage_for_setup_phase / 1.0e+06 << "\n"
172 << std::endl;
173 } // if (Compute_memory_statistics)
174 } // End of setup
175
176
177 //=============================================================================
178 /// Preconditioner solve for the exact preconditioner
179 //=============================================================================
180 template<typename MATRIX>
182 const DoubleVector& r, DoubleVector& z)
183 {
184 // Vector of vectors for each section of residual vector
185 Vector<DoubleVector> block_r;
186
187 // Rearrange the vector r into the vector of block vectors block_r
188 this->get_block_vectors(r, block_r);
189
190 // Allocate space for the properly rearranged RHS vector
191 DoubleVector rhs_reordered;
192
193 // Concatenate the DoubleVectors together
194 DoubleVectorHelpers::concatenate(block_r, rhs_reordered);
195
196 // Allocate space for the rearranged solution vector
197 DoubleVector z_reordered;
198
199 // Solve the whole system exactly
200 this->Subsidiary_preconditioner_pt[0]->preconditioner_solve(rhs_reordered,
201 z_reordered);
202
203 // Vector of vectors for each section of solution vector (copy the RHS
204 // block vector to get the sizing and distributions right)
205 Vector<DoubleVector> block_z(block_r);
206
207 // Split the solution vector into blocks
208 DoubleVectorHelpers::split(z_reordered, block_z);
209
210 // Copy the solution from the block vector block_z back into z
211 this->return_block_vectors(block_z, z);
212 } // End of preconditioner_solve
213
214 //============================================================================
215 /// Set up the block triangular preconditioner
216 //============================================================================
217 template<typename MATRIX>
219 {
220 // Clean the memory
221 this->clean_up_memory();
222
223 // Subsidiary preconditioners don't really need the meshes
224 if (this->is_master_block_preconditioner())
225 {
226#ifdef PARANOID
227 if (this->gp_nmesh() == 0)
228 {
229 std::ostringstream err_msg;
230 err_msg << "There are no meshes set.\n"
231 << "Did you remember to call add_mesh(...)?";
232 throw OomphLibError(
233 err_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
234 }
235#endif
236
237 // Set all meshes if this is master block preconditioner
238 this->gp_preconditioner_set_all_meshes();
239 }
240
241 // Initialise the memory usage variable
242 Memory_usage_in_bytes = 0.0;
243
244 // If we're meant to build silently
245 if (this->Silent_preconditioner_setup == true)
246 {
247 // Store the output stream pointer
248 this->Stream_pt = oomph_info.stream_pt();
249
250 // Now set the oomph_info stream pointer to the null stream to
251 // disable all possible output
253 }
254
255 // Set up the block look up schemes
256 this->gp_preconditioner_block_setup();
257
258 // Number of block types
259 unsigned n_block_types = this->nblock_types();
260
261 // Storage for the off diagonal matrix vector products
262 Off_diagonal_matrix_vector_products.resize(n_block_types, n_block_types, 0);
263
264 // Fill in any null subsidiary preconditioners
265 this->fill_in_subsidiary_preconditioners(n_block_types);
266
267 // The total time for extracting all the blocks from the "global" matrix
268 double t_extraction_total = 0.0;
269
270 // The total time for setting up the subsidiary preconditioners
271 double t_subsidiary_setup_total = 0.0;
272
273 // The total time for setting up the matrix-vector products
274 double t_mvp_setup_total = 0.0;
275
276 // Build the preconditioners and matrix vector products
277 for (unsigned i = 0; i < n_block_types; i++)
278 {
279 // See if the i-th subsidiary preconditioner is a block preconditioner
280 if (dynamic_cast<BlockPreconditioner<CRDoubleMatrix>*>(
281 this->Subsidiary_preconditioner_pt[i]) != 0)
282 {
283 // If we need to compute the memory statistics
284 if (Compute_memory_statistics)
285 {
286 // If we're dealing with a GMRES block preconditioner
287 if (dynamic_cast<GMRESBlockPreconditioner*>(
288 this->Subsidiary_preconditioner_pt[i]))
289 {
290 // Enable the doc-ing of memory usage in the GMRES block
291 // preconditioner
292 dynamic_cast<GMRESBlockPreconditioner*>(
293 this->Subsidiary_preconditioner_pt[i])
294 ->enable_doc_memory_usage();
295 }
296 } // if (Compute_memory_statistics)
297
298 // Get the start time
299 double t_subsidiary_setup_start = TimingHelpers::timer();
300
301 // It's a block preconditioner so pass it a pointer to the global
302 // matrix!
303 this->Subsidiary_preconditioner_pt[i]->setup(this->matrix_pt());
304
305 // Get the end time
306 double t_subsidiary_setup_end = TimingHelpers::timer();
307
308 // Update the timing total
309 t_subsidiary_setup_total +=
310 (t_subsidiary_setup_end - t_subsidiary_setup_start);
311 }
312 // It's just a preconditioner so pass it a deep copy of the block!
313 else
314 {
315 // Get the start time
316 double t_extract_start = TimingHelpers::timer();
317
318 // Grab the i-th diagonal block
319 CRDoubleMatrix block_matrix = this->get_block(i, i);
320
321 // Get the end time
322 double t_extract_end = TimingHelpers::timer();
323
324 // Update the timing total
325 t_extraction_total += (t_extract_end - t_extract_start);
326
327 // Get the start time
328 double t_subsidiary_setup_start = TimingHelpers::timer();
329
330 // Set up the i-th subsidiary preconditioner with this block
331 this->Subsidiary_preconditioner_pt[i]->setup(&block_matrix);
332
333 // Get the end time
334 double t_subsidiary_setup_end = TimingHelpers::timer();
335
336 // Update the timing total
337 t_subsidiary_setup_total +=
338 (t_subsidiary_setup_end - t_subsidiary_setup_start);
339 }
340
341 // The index of the lower column bound
342 unsigned l;
343
344 // The index of the upper column bound
345 unsigned u;
346
347 // If we're using the upper triangular form of the preconditioner
348 if (Upper_triangular)
349 {
350 // The lower bound column bound
351 l = i + 1;
352
353 // If the block bandwidth has been provided
354 if (Block_bandwidth >= 0)
355 {
356 // The upper bound is the last nonzero block column
357 u = std::min(n_block_types, l + Block_bandwidth);
358 }
359 // Default case; loop over all the block columns
360 else
361 {
362 // The upper bound is the last column
363 u = n_block_types;
364 }
365 }
366 // If we're using the lower triangular form of the preconditioner
367 else
368 {
369 // The upper bound column bound
370 u = i;
371
372 // If the block bandwidth has been provided
373 if (Block_bandwidth >= 0)
374 {
375 // The lower bound is the last nonzero block column. The explicit
376 // casting is a bit over the top but it's better to be safe than
377 // sorry...
378 l = std::max(0, int(int(u) - Block_bandwidth));
379 }
380 // Default case; loop over all the block columns
381 else
382 {
383 // The lower bound is the first column
384 l = 0;
385 }
386 } // for (unsigned i=0;i<nblock_types;i++)
387
388 // Loop over the chosen block columns
389 for (unsigned j = l; j < u; j++)
390 {
391 // Get the start time
392 double t_extract_start = TimingHelpers::timer();
393
394 // Get the (i,j)-th block matrix
395 CRDoubleMatrix block_matrix = this->get_block(i, j);
396
397 // Do we need to doc. the memory statistics?
398 if (Compute_memory_statistics)
399 {
400 // How many rows are there in this block?
401 unsigned n_row = block_matrix.nrow();
402
403 // How many nonzeros are there in this block?
404 unsigned n_nnz = block_matrix.nnz();
405
406 // Compute the memory usage. The only memory overhead here (for the
407 // setup phase) is the storage of the sub-diagonal blocks for the MVPs
408 Memory_usage_in_bytes += ((2 * ((n_row + 1) * sizeof(int))) +
409 (n_nnz * (sizeof(int) + sizeof(double))));
410 }
411
412 // Get the end time
413 double t_extract_end = TimingHelpers::timer();
414
415 // Update the timing total
416 t_extraction_total += (t_extract_end - t_extract_start);
417
418 // Get the start time
419 double t_mvp_start = TimingHelpers::timer();
420
421 // Copy the block into a "multiplier" class. If trilinos is being
422 // used this should also be faster than oomph-lib's multiphys.
423 Off_diagonal_matrix_vector_products(i, j) = new MatrixVectorProduct();
424
425 // Set the damn thing up
426 this->setup_matrix_vector_product(
427 Off_diagonal_matrix_vector_products(i, j), &block_matrix, j);
428
429 // Get the end time
430 double t_mvp_end = TimingHelpers::timer();
431
432 // Update the timing total
433 t_mvp_setup_total += (t_mvp_end - t_mvp_start);
434 }
435 } // for (unsigned i=0;i<n_block_types;i++)
436
437 // Remember that the preconditioner has been set up
438 Preconditioner_has_been_setup = true;
439
440 // If we're meant to build silently, reassign the oomph stream pointer
441 if (this->Silent_preconditioner_setup == true)
442 {
443 // Store the output stream pointer
444 oomph_info.stream_pt() = this->Stream_pt;
445
446 // Reset our own stream pointer
447 this->Stream_pt = 0;
448 }
449
450 // Tell the user
451 oomph_info << "Total block extraction time [sec]: " << t_extraction_total
452 << "\nTotal subsidiary preconditioner setup time [sec]: "
453 << t_subsidiary_setup_total
454 << "\nTotal matrix-vector product setup time [sec]: "
455 << t_mvp_setup_total << std::endl;
456
457 // Do we need to doc. the memory statistics?
458 // NOTE: We're going to assume that either:
459 // (1) GMRES is used as a subsidiary block preconditioner (preconditioned
460 // by the Navier-Stokes subsidiary block preconditioner), or;
461 // (2) SuperLU is used as the subsidiary preconditioner.
462 if (Compute_memory_statistics)
463 {
464 // Allocate space for the total memory usage
465 double total_memory_usage_for_setup_phase = 0.0;
466
467 // How many rows are there in the global Jacobian?
468 unsigned n_row = this->matrix_pt()->nrow();
469
470 // How many nonzeros are there in the global Jacobian?
471 unsigned n_nnz = this->matrix_pt()->nnz();
472
473 // Storage for the memory usage of the global system matrix.
474 // NOTE: This calculation is done by taking into account the storage
475 // scheme for a compressed row matrix
476 double memory_usage_for_storage_of_global_jacobian_in_bytes =
477 ((2 * ((n_row + 1) * sizeof(int))) +
478 (n_nnz * (sizeof(int) + sizeof(double))));
479
480 // Allocate storage for the memory usage of the subsidiary preconditioner
481 double memory_usage_for_subsidiary_preconditioner_in_bytes = 0.0;
482
483 // Allocate storage for the memory usage of the Navier-Stokes subsidiary
484 // block preconditioner (if it's used, that is)
485 double memory_usage_for_nssbp_preconditioner_in_bytes = 0.0;
486
487 // Boolean to indicate whether we've issue a warning about not being
488 // able to compute memory statistics. If it doesn't work with one
489 // block, chances are it won't work with any of them, so don't be a
490 // pain in the gluteus maximus and just spit out the warning once
491 bool have_issued_warning_about_memory_calculations = false;
492
493 // Build the preconditioners and matrix vector products
494 for (unsigned i = 0; i < n_block_types; i++)
495 {
496 // See if the i-th subsidiary preconditioner the GMRES block
497 // preconditioner
498 if (dynamic_cast<GMRESBlockPreconditioner*>(
499 this->Subsidiary_preconditioner_pt[i]) != 0)
500 {
501 // Upcast the subsidiary preconditioner pointer to a GMRES block
502 // preconditioner pointer
503 GMRESBlockPreconditioner* gmres_block_prec_pt =
504 dynamic_cast<GMRESBlockPreconditioner*>(
505 this->Subsidiary_preconditioner_pt[i]);
506
507 // Update the subsidiary preconditioner memory usage variable
508 memory_usage_for_subsidiary_preconditioner_in_bytes +=
509 gmres_block_prec_pt->get_memory_usage_in_bytes();
510
511 // Add on the memory usage for the NS subsidiary block preconditioner
512 memory_usage_for_nssbp_preconditioner_in_bytes +=
513 gmres_block_prec_pt->navier_stokes_subsidiary_preconditioner_pt()
515 }
516 // This block is preconditioned by the SuperLU preconditioner
517 else if (dynamic_cast<SuperLUPreconditioner*>(
518 this->Subsidiary_preconditioner_pt[i]) != 0)
519 {
520 // Update the subsidiary preconditioner memory usage variable
521 memory_usage_for_subsidiary_preconditioner_in_bytes +=
522 dynamic_cast<SuperLUPreconditioner*>(
523 this->Subsidiary_preconditioner_pt[i])
524 ->get_total_memory_needed_for_superlu();
525 }
526 // Don't know what to do otherwise!
527 else
528 {
529 // Have we already complained once?
530 if (!have_issued_warning_about_memory_calculations)
531 {
532 // Remember that we've issued a warning now
533 have_issued_warning_about_memory_calculations = true;
534
535 // Allocate storage for an output stream
536 std::ostringstream warning_message_stream;
537
538 // Create a warning message
539 warning_message_stream
540 << "Can't compute the memory statistics for the " << i
541 << "-th diagonal block in the banded\nblock "
542 << "triangular preconditioner so I'm just going "
543 << "to skip it." << std::endl;
544
545 // Give the user a warning
546 OomphLibWarning(warning_message_stream.str(),
547 OOMPH_CURRENT_FUNCTION,
548 OOMPH_EXCEPTION_LOCATION);
549 }
550 }
551 } // for (unsigned i=0;i<nblock_types;i++)
552
553 // Create some space
554 oomph_info << std::endl;
555
556 // If we've got a nonzero contribution from the NSSBP
557 if (memory_usage_for_nssbp_preconditioner_in_bytes > 0.0)
558 {
559 // Doc. the memory usage for the NSSBP
560 oomph_info << "Amount of memory used by Navier-Stokes subsidiary block "
561 << "preconditioner (MB): "
562 << memory_usage_for_nssbp_preconditioner_in_bytes / 1.0e+06
563 << std::endl;
564
565 // Add in the NSSBP contribution
566 total_memory_usage_for_setup_phase +=
567 memory_usage_for_nssbp_preconditioner_in_bytes;
568 }
569
570 // Add in the subsidiary preconditioners contribution
571 total_memory_usage_for_setup_phase +=
572 memory_usage_for_subsidiary_preconditioner_in_bytes;
573
574 // Add in the banded block triangular preconditioner contribution
575 total_memory_usage_for_setup_phase += get_memory_usage_in_bytes();
576
577 // Add in the global Jacobian contribution
578 total_memory_usage_for_setup_phase +=
579 memory_usage_for_storage_of_global_jacobian_in_bytes;
580
581 // How much memory have we used in the subsidiary preconditioners?
583 << "Amount of memory used by the subsidiary preconditioners "
584 << "(MB): "
585 << memory_usage_for_subsidiary_preconditioner_in_bytes / 1.0e+06
586 << "\nAmount of memory used by banded block triangular "
587 << "preconditioner (MB): " << get_memory_usage_in_bytes() / 1.0e+06
588 << "\nAmount of memory used to store (global) Jacobian (MB): "
589 << memory_usage_for_storage_of_global_jacobian_in_bytes / 1.0e+06
590 << "\nTotal amount of memory being used after setup (MB): "
591 << total_memory_usage_for_setup_phase / 1.0e+06 << "\n"
592 << std::endl;
593 } // if (Compute_memory_statistics)
594 } // End of setup
595
596 //=============================================================================
597 /// Preconditioner solve for the block triangular preconditioner
598 //=============================================================================
599 template<typename MATRIX>
601 const DoubleVector& r, DoubleVector& z)
602 {
603 // Cache number of block types
604 unsigned n_block = this->nblock_types();
605
606 // The start index
607 int start;
608
609 // The end index
610 int end;
611
612 // The max. bandwidth of the block matrix
613 int bandwidth_limit;
614
615 // The steps to take (=+1/-1) which is dependent on whether we're using
616 // the lower or upper triangular form of the preconditioner
617 int step;
618
619 // If we're using this as an block upper triangular preconditioner
620 if (Upper_triangular)
621 {
622 start = n_block - 1;
623 end = -1;
624 step = -1;
625 }
626 else
627 {
628 start = 0;
629 end = n_block;
630 step = 1;
631 }
632
633 // Storage for the iteration count. The entries will be updated as we loop
634 // over the blocks. Here -1 is used to indicate that a block was not solved
635 // by GMRES so we ignore it
636 Vector<int> iteration_count(n_block, -1);
637
638 // Storage for the number of blocks solved with GMRES
639 unsigned n_block_solved_with_gmres = 0;
640
641 // Vector of vectors for each section of residual vector
642 Vector<DoubleVector> block_r;
643
644 // Rearrange the vector r into the vector of block vectors block_r
645 this->get_block_vectors(r, block_r);
646
647 // Vector of vectors for the solution block vectors
648 Vector<DoubleVector> block_z(n_block);
649
650 // Loop over the block rows
651 for (int i = start; i != end; i += step)
652 {
653 // If the i-th subsidiary preconditioner is a block preconditioner
654 if (dynamic_cast<BlockPreconditioner<CRDoubleMatrix>*>(
655 this->Subsidiary_preconditioner_pt[i]) != 0)
656 {
657 // This is pretty inefficient but there is no other way to do this with
658 // block sub preconditioners as far as I can tell because they demand
659 // to have the full r and z vectors...
660 DoubleVector block_z_with_size_of_full_z(r.distribution_pt());
661 DoubleVector r_updated(r.distribution_pt());
662
663 // Construct the new r vector (with the update given by back subs
664 // below).
665 this->return_block_vectors(block_r, r_updated);
666
667 // Solve (blocking is handled inside the block preconditioner).
668 this->Subsidiary_preconditioner_pt[i]->preconditioner_solve(
669 r_updated, block_z_with_size_of_full_z);
670
671 // Extract this block's z vector because we need it for the next step
672 // (and throw away block_z_with_size_of_full_z).
673 this->get_block_vector(i, block_z_with_size_of_full_z, block_z[i]);
674 }
675 // If the subsidiary preconditioner is just a regular preconditioner
676 else
677 {
678 // Solve on the block
679 this->Subsidiary_preconditioner_pt[i]->preconditioner_solve(block_r[i],
680 block_z[i]);
681 }
682
683 // See if the i-th subsidiary preconditioner the GMRES block
684 // preconditioner
685 if (dynamic_cast<GMRESBlockPreconditioner*>(
686 this->Subsidiary_preconditioner_pt[i]) != 0)
687 {
688 // Update the iteration count
689 iteration_count[i] = dynamic_cast<GMRESBlockPreconditioner*>(
690 this->Subsidiary_preconditioner_pt[i])
691 ->iterations();
692
693 // Increment the counter
694 n_block_solved_with_gmres++;
695 }
696
697 // If the block bandwidth has been provided
698 if (Block_bandwidth >= 0)
699 {
700 // If we're using this as an block upper triangular preconditioner
701 if (Upper_triangular)
702 {
703 bandwidth_limit = std::max(end - i, step + step * Block_bandwidth);
704 }
705 else
706 {
707 bandwidth_limit = std::min(end - i, step + step * Block_bandwidth);
708 }
709 }
710 // Default case; loop over all the block rows
711 else
712 {
713 bandwidth_limit = end - i;
714 }
715
716 // Substitute (over all the rows)
717 for (int j = i + step; j != i + bandwidth_limit; j += step)
718 {
719 // Allocate space for the matrix-vector product (MVP)
720 DoubleVector temp;
721
722 // Calculate the MVP
723 Off_diagonal_matrix_vector_products(j, i)->multiply(block_z[i], temp);
724
725 // Now update the RHS vector
726 block_r[j] -= temp;
727 }
728 } // for (int i=start;i!=end;i+=step)
729
730 // Copy the solution from the block vector block_z back into z
731 this->return_block_vectors(block_z, z);
732
733 // Only compute the iteration count statistics if we have some data
734 if (n_block_solved_with_gmres > 0)
735 {
736 // Storage for the average iteration count
737 double n_iter_mean = 0.0;
738
739 // Storage for the variance of the iteration count
740 double n_iter_var = 0.0;
741
742 // Loop over the entries of the iteration count vector
743 for (unsigned i = 0; i < n_block; i++)
744 {
745 // If we used GMRES on this block
746 if (dynamic_cast<GMRESBlockPreconditioner*>(
747 this->Subsidiary_preconditioner_pt[i]) != 0)
748 {
749 // Sanity check: make sure there's a non-negative iteration count
750 if (iteration_count[i] < 0)
751 {
752 // Create an output stream
753 std::ostringstream error_message_stream;
754
755 // Create an error message
756 error_message_stream
757 << "Iteration count should be non-negative but "
758 << "you have an iteration\ncount of " << iteration_count[i] << "!"
759 << std::endl;
760
761 // Throw an error
762 throw OomphLibError(error_message_stream.str(),
763 OOMPH_CURRENT_FUNCTION,
764 OOMPH_EXCEPTION_LOCATION);
765 }
766 } // if (dynamic_cast<GMRESBlockPreconditioner*>...
767 } // for (unsigned i=0;i<n_block;i++)
768
769 // Loop over the entries of the iteration count vector
770 for (unsigned i = 0; i < n_block; i++)
771 {
772 // If we used GMRES on this block
773 if (dynamic_cast<GMRESBlockPreconditioner*>(
774 this->Subsidiary_preconditioner_pt[i]) != 0)
775 {
776 // Update the mean
777 n_iter_mean +=
778 (double(iteration_count[i]) / double(n_block_solved_with_gmres));
779 } // if (dynamic_cast<GMRESBlockPreconditioner*>...
780 } // for (unsigned i=0;i<n_block;i++)
781
782 // Loop over the entries of the iteration count vector (again)
783 // NOTE: We don't need to do the error checks again, we did them on the
784 // previous run through
785 for (unsigned i = 0; i < n_block; i++)
786 {
787 // If we used GMRES on this block
788 if (dynamic_cast<GMRESBlockPreconditioner*>(
789 this->Subsidiary_preconditioner_pt[i]) != 0)
790 {
791 // Update the variance
792 n_iter_var +=
793 (std::pow(double(iteration_count[i]) - n_iter_mean, 2.0) /
794 double(n_block_solved_with_gmres));
795 } // if (dynamic_cast<GMRESBlockPreconditioner*>...
796 } // for (unsigned i=0;i<n_block;i++)
797
798 // Doc. the statistics
799 oomph_info << "\nNumber of subsidiary blocks solved with GMRES: "
800 << n_block_solved_with_gmres
801 << "\nAverage subsidiary preconditioner iteration count: "
802 << n_iter_mean
803 << "\nStandard deviation of subsidiary preconditioner "
804 << "iteration count: " << std::sqrt(n_iter_var) << std::endl;
805 } // if (n_block_solved_with_gmres>0)
806 } // End of preconditioner_solve
807
808 // Ensure build of required objects (BUT only for CRDoubleMatrix objects)
811} // End of namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
General purpose block triangular preconditioner. By default this operates as an upper triangular prec...
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
unsigned long nnz() const
Return the number of nonzero entries (the local nnz)
Definition: matrices.h:1096
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:1002
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
General purpose block tridiagonal preconditioner. By default SuperLUPreconditioner (or SuperLUDistPre...
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to r.
The block preconditioner form of GMRES. This version extracts the blocks from the global systems and ...
SpaceTimeNavierStokesSubsidiaryPreconditioner * navier_stokes_subsidiary_preconditioner_pt() const
Handle to the Navier-Stokes subsidiary block preconditioner DRAIG: Make sure the desired const-ness i...
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....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
An interface to allow SuperLU to be used as an (exact) Preconditioner.
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 split(const DoubleVector &in_vector, Vector< DoubleVector * > &out_vector_pt)
Split a DoubleVector into the out DoubleVectors. Let vec_A be the in Vector, and let vec_B and vec_C ...
void concatenate(const Vector< DoubleVector * > &in_vector_pt, DoubleVector &out_vector)
Concatenate DoubleVectors. Takes a Vector of DoubleVectors. If the out vector is built,...
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...
void get_memory_usage_in_bytes(double *lu_factor_memory, double *total_memory)
Definition: superlu.c:85