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// Version 1.0; svn revision $LastChangedRevision: 1132 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-02-23 05:43:26 +0000 (Tue, 23 Feb 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 
31 #include "oomph_definitions.h"
32 #include "preconditioner.h"
33 #include "block_preconditioner.h"
34 #include "matrices.h"
36 
37 namespace oomph
38 {
39 
40  //============================================================================
41  /// setup for the block diagonal preconditioner
42  //============================================================================
43  template<typename MATRIX>
45  {
46  // clean the memory
47  this->clean_up_memory();
48 
49  // Note: Subsidiary block preconditioners don't really use their
50  // mesh pointers (since the lookup schemes inferred from them are
51  // set up by their masters). Generally we insist that (for uniformity)
52  // a mesh should be set, but sometimes subsidiary preconditioners are
53  // set up on the fly and we can't sensibly set meshes, so we're forgiving...)
54  if (this->is_master_block_preconditioner())
55  {
56 #ifdef PARANOID
57  if( this->gp_nmesh() == 0)
58  {
59  std::ostringstream err_msg;
60  err_msg << "There are no meshes set.\n"
61  << "Did you remember to call add_mesh(...)?";
62  throw OomphLibError(err_msg.str(),
63  OOMPH_CURRENT_FUNCTION,
64  OOMPH_EXCEPTION_LOCATION);
65  }
66 #endif
67 
68  // Set all meshes if this is master block preconditioner
69  this->gp_preconditioner_set_all_meshes();
70  }
71 
72  // Set up the block look up schemes
73  this->gp_preconditioner_block_setup();
74 
75  // number of types of blocks.
76  unsigned nblock_types = this->nblock_types();
77 
78  // Create any subsidiary preconditioners needed
79  this->fill_in_subsidiary_preconditioners(nblock_types);
80 
81  // If using two level parallelisation then we need to use a
82  // PrecondtionerArray which requires very different setup. ??ds possibly
83  // it should have it's own class?
84  if(Use_two_level_parallelisation)
85  {
86  // Get the blocks. We have to use new because you can't have containers
87  // of matrices because the copy constructor is buggy (so we create a
88  // container of pointers instead). ??ds
89  Vector<CRDoubleMatrix*> block_diagonal_matrix_pt(nblock_types, 0);
90  for(unsigned i=0; i<nblock_types; i++)
91  {
92  block_diagonal_matrix_pt[i] = new CRDoubleMatrix;
93  this->get_block(i, get_other_diag_ds(i, nblock_types),
94  *block_diagonal_matrix_pt[i]);
95  }
96 
97  // Construct the preconditioner array
98  Preconditioner_array_pt = new PreconditionerArray;
99 
100  Preconditioner_array_pt->
101  setup_preconditioners(block_diagonal_matrix_pt,
102  this->Subsidiary_preconditioner_pt,
103  this->comm_pt());
104 
105  // and delete the blocks
106  for(unsigned i=0; i<nblock_types; i++)
107  {
108  delete block_diagonal_matrix_pt[i];
109  block_diagonal_matrix_pt[i] = 0;
110  }
111 
112  // Preconditioner array is weird, it calls delete on all the
113  // preconditioners you give it and requires new ones each time!
114  this->Subsidiary_preconditioner_pt.clear();
115  }
116 
117 
118  // Otherwise just set up each block's preconditioner in order
119  else
120  {
121  for(unsigned i=0; i<nblock_types; i++)
122  {
123  // Get the block
124  unsigned j = get_other_diag_ds(i, nblock_types);
125  CRDoubleMatrix block = this->get_block(i, j);
126 
127  // Set up preconditioner (i.e. approximately solve the block + store solution)
128  double superlusetup_start = TimingHelpers::timer();
129  this->Subsidiary_preconditioner_pt[i]->setup(&block);
130  double superlusetup_end = TimingHelpers::timer();
131  oomph_info << "Took " << superlusetup_end - superlusetup_start
132  << "s to setup."<< std::endl;
133  }
134  }
135 
136  }
137 
138  //=============================================================================
139  /// Preconditioner solve for the block diagonal preconditioner
140  //=============================================================================
141  template<typename MATRIX>
144  {
145  // Cache umber of block types
146  unsigned n_block = this->nblock_types();
147 
148  // Get the right hand side vector (b in Ax = b) in block form.
149  Vector<DoubleVector> block_r;
150  this->get_block_vectors(r, block_r);
151 
152  // Make sure z vector is built
153  if (!z.built()) { z.build(this->distribution_pt(),0.0); }
154 
155  // Vector of vectors for storage of solution in block form.
156  Vector<DoubleVector> block_z(n_block);
157 
158  if (Use_two_level_parallelisation)
159  {
160  Preconditioner_array_pt->solve_preconditioners(block_r, block_z);
161  }
162  else
163  {
164  // solve each diagonal block
165  for (unsigned i = 0; i < n_block; i++)
166  {
167  double t_start=0.0;
168  if (Doc_time_during_preconditioner_solve)
169  {
170  t_start=TimingHelpers::timer();
171  }
172  this->Subsidiary_preconditioner_pt[i]->preconditioner_solve(block_r[i],
173  block_z[i]);
174  if (Doc_time_during_preconditioner_solve)
175  {
176  oomph_info << "Time for application of " << i
177  << "-th block preconditioner: "
178  << TimingHelpers::timer()-t_start
179  << std::endl;
180  }
181  }
182  }
183 
184  // copy solution in block vectors block_z back to z
185  this->return_block_vectors(block_z,z);
186  }
187 
188 
189  //============================================================================
190  /// setup for the block triangular preconditioner
191  //============================================================================
192  template<typename MATRIX>
195  {
196  // clean the memory
197  this->clean_up_memory();
198 
199  // Subsidiary preconditioners don't really need the meshes
200  if (this->is_master_block_preconditioner())
201  {
202 #ifdef PARANOID
203  if( this->gp_nmesh() == 0)
204  {
205  std::ostringstream err_msg;
206  err_msg << "There are no meshes set.\n"
207  << "Did you remember to call add_mesh(...)?";
208  throw OomphLibError(err_msg.str(),
209  OOMPH_CURRENT_FUNCTION,
210  OOMPH_EXCEPTION_LOCATION);
211  }
212 #endif
213 
214  // Set all meshes if this is master block preconditioner
215  this->gp_preconditioner_set_all_meshes();
216  }
217 
218  // Set up the block look up schemes
219  this->gp_preconditioner_block_setup();
220 
221  // number of block types
222  unsigned nblock_types = this->nblock_types();
223 
224  // storage for the off diagonal matrix vector products
225  Off_diagonal_matrix_vector_products.resize(nblock_types,nblock_types,0);
226 
227  // Fill in any null subsidiary preconditioners
228  this->fill_in_subsidiary_preconditioners(nblock_types);
229 
230  // build the preconditioners and matrix vector products
231  for(unsigned i = 0; i < nblock_types; i++)
232  {
233  // Get the block and set up the preconditioner.
234  {
235  CRDoubleMatrix block_matrix = this->get_block(i,i);
236  this->Subsidiary_preconditioner_pt[i]
237  ->setup(&block_matrix);
238  }
239 
240  // next setup the off diagonal mat vec operators
241  unsigned l, u;
242  if(Upper_triangular)
243  {
244  l = i+1;
245  u = nblock_types;
246  }
247  else
248  {
249  l = 0;
250  u = i;
251  }
252 
253  for(unsigned j = l; j < u; j++)
254  {
255  // Get the block
256  CRDoubleMatrix block_matrix = this->get_block(i,j);
257 
258  // Copy the block into a "multiplier" class. If trilinos is being
259  // used this should also be faster than oomph-lib's multiphys.
260  Off_diagonal_matrix_vector_products(i,j) = new MatrixVectorProduct();
261 
262  this->setup_matrix_vector_product(
263  Off_diagonal_matrix_vector_products(i,j),
264  &block_matrix,j);
265  }
266  }
267  }
268 
269  //=============================================================================
270  /// Preconditioner solve for the block triangular preconditioner
271  //=============================================================================
272  template<typename MATRIX> void BlockTriangularPreconditioner<MATRIX>::
274  {
275  // Cache number of block types
276  unsigned n_block = this->nblock_types();
277 
278  //
279  int start, end, step;
280  if(Upper_triangular)
281  {
282  start = n_block-1;
283  end = -1;
284  step = -1;
285  }
286  else
287  {
288  start = 0;
289  end = n_block;
290  step = 1;
291  }
292 
293  // vector of vectors for each section of residual vector
294  Vector<DoubleVector> block_r;
295 
296  // rearrange the vector r into the vector of block vectors block_r
297  this->get_block_vectors(r,block_r);
298 
299  // vector of vectors for the solution block vectors
300  Vector<DoubleVector> block_z(n_block);
301 
302  //
303  for (int i = start; i != end; i+=step)
304  {
305  // ??ds ugly, fix this?
306  if(dynamic_cast<BlockPreconditioner<CRDoubleMatrix>*>
307  (this->Subsidiary_preconditioner_pt[i]) == 0)
308  {
309  // solve on the block
310  this->Subsidiary_preconditioner_pt[i]->
311  preconditioner_solve(block_r[i], block_z[i]);
312  }
313  else
314  {
315  // This is pretty inefficient but there is no other way to do this with
316  // block sub preconditioners as far as I can tell because they demand
317  // to have the full r and z vectors...
318 
319  DoubleVector block_z_with_size_of_full_z(r.distribution_pt()),
320  r_updated(r.distribution_pt());
321 
322  // Construct the new r vector (with the update given by back subs
323  // below).
324  this->return_block_vectors(block_r, r_updated);
325 
326  // Solve (blocking is handled inside the block preconditioner).
327  this->Subsidiary_preconditioner_pt[i]->
328  preconditioner_solve(r_updated, block_z_with_size_of_full_z);
329 
330  // Extract this block's z vector because we need it for the next step
331  // (and throw away block_z_with_size_of_full_z).
332  this->get_block_vector(i, block_z_with_size_of_full_z, block_z[i]);
333  }
334 
335  // substitute
336  for (int j = i + step; j !=end; j+=step)
337  {
338  DoubleVector temp;
339  Off_diagonal_matrix_vector_products(j,i)->multiply(block_z[i],temp);
340  block_r[j] -= temp;
341  }
342  }
343 
344  // copy solution in block vectors block_r back to z
345  this->return_block_vectors(block_z,z);
346  }
347 
348 
349 
350 
351  //=============================================================================
352  /// Setup for the exact block preconditioner
353  //=============================================================================
354  template<typename MATRIX>
356  {
357 
358  // If this is a master block preconditioner,
359  // we give the mesh pointers to the BlockPreconditioner base class.
360  // Subsidiary preconditioners don't need meshes...
361  if(this->is_master_block_preconditioner())
362  {
363 #ifdef PARANOID
364  if( this->gp_nmesh() == 0)
365  {
366  std::ostringstream err_msg;
367  err_msg << "There are no meshes set.\n"
368  << "Did you remember to call add_mesh(...)?";
369  throw OomphLibError(err_msg.str(),
370  OOMPH_CURRENT_FUNCTION,
371  OOMPH_EXCEPTION_LOCATION);
372  }
373 #endif
374 
375  // Set all meshes if this is master block preconditioner
376  this->gp_preconditioner_set_all_meshes();
377  }
378 
379  // Set up the block look up schemes
380  this->gp_preconditioner_block_setup();
381 
382  // get the number of DOF types
383  unsigned nblock_types = this->nblock_types();
384 
385  // Build the preconditioner matrix
386  VectorMatrix<BlockSelector> required_blocks(nblock_types,nblock_types);
387 
388  // boolean indicating if we want the block, stored for readability.
389  const bool want_block = true;
390  for (unsigned b_i = 0; b_i < nblock_types; b_i++)
391  {
392  for (unsigned b_j = 0; b_j < nblock_types; b_j++)
393  {
394  required_blocks[b_i][b_j].select_block(b_i,b_j,want_block);
395  }
396  }
397 
398  // Get the concatenated blocks
399  CRDoubleMatrix exact_block_matrix
400  = this->get_concatenated_block(required_blocks);
401 
402  // Setup the preconditioner.
403  this->fill_in_subsidiary_preconditioners(1); // Only need one preconditioner
404  preconditioner_pt()->setup(&exact_block_matrix);
405  }
406 
407  //=============================================================================
408  /// Preconditioner solve for the exact block preconditioner
409  //=============================================================================
410  template<typename MATRIX>
413  {
414  // get the block ordered components of the r vector for this preconditioner
415  DoubleVector block_order_r;
416  this->get_block_ordered_preconditioner_vector(r,block_order_r);
417 
418  // vector for solution
419  DoubleVector block_order_z;
420 
421  // apply the preconditioner
422  preconditioner_pt()->preconditioner_solve(block_order_r,block_order_z);
423 
424  // copy solution back to z vector
425  this->return_block_ordered_preconditioner_vector(block_order_z,z);
426  }
427 
433 
434 }
bool built() const
cstr elem_len * i
Definition: cfortran.h:607
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to r.
OomphInfo oomph_info
PreconditionerArray - NOTE - first implementation, a number of assumptions / simplifications were mad...
Block "anti-diagonal" preconditioner, i.e. same as block diagonal but along the other diagonal of the...
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to r.
void start(const unsigned &i)
(Re-)start i-th timer
double timer()
returns the time in seconds after some point in past
void preconditioner_solve(const DoubleVector &r, DoubleVector &z)
Apply preconditioner to r.
Block diagonal preconditioner. By default SuperLU is used to solve the subsidiary systems...
Matrix vector product helper class - primarily a wrapper to Trilinos's Epetra matrix vector product m...
A vector in the mathematical sense, initially developed for linear algebra type applications. If MPI then this vector can be distributed - its distribution is described by the LinearAlgebraDistribution object at Distribution_pt. Data is stored in a C-style pointer vector (double*)
Definition: double_vector.h:61
General purpose block triangular preconditioner By default this is Upper triangular. By default SuperLUPreconditioner (or SuperLUDistPreconditioner) is used to solve the subsidiary systems, but other preconditioners can be used by setting them using passing a pointer to a function of type SubsidiaryPreconditionerFctPt to the method subsidiary_preconditioner_function_pt().
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:872
virtual void setup()
Setup the preconditioner.