matrix_vector_product.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: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 #include "matrix_vector_product.h"
31 
32 namespace oomph
33 {
34 
35 
36  //============================================================================
37  /// \short Setup the matrix vector product operator.
38  /// WARNING: This class is wrapper to Trilinos Epetra matrix vector
39  /// multiply methods, if Trilinos is not installed then this class will
40  /// function as expected, but there will be no computational speed gain.
41  /// By default the Epetra_CrsMatrix::multiply(...) are employed.
42  /// The optional argument col_dist_pt is the distribution of:
43  /// x if using multiply(...) or y if using multiply_transpose(...)
44  /// where this is A x = y. By default, this is assumed to the uniformly
45  /// distributed based on matrix_pt->ncol().
46  //============================================================================
48  const LinearAlgebraDistribution* col_dist_pt)
49  {
50  // clean memory
51  this->clean_up_memory();
52 
53  // (re)build distribution_pt
54  this->build_distribution(matrix_pt->distribution_pt());
55 
56  // store the number of columns
57  Ncol = matrix_pt->ncol();
58 
59  // determine whether we are using trilinos
60  Using_trilinos=false;
61 #ifdef OOMPH_HAS_TRILINOS
62 #ifdef OOMPH_HAS_MPI
64  {
65  Using_trilinos=true;
66  }
67 #else
68  Using_trilinos=true;
69 #endif
70 #endif
71 
72  // create the column distribution map, if a distribution is not provided,
73  // create a uniformly distributed based on matrix_pt->ncol().
74  // Otherwise, use the provided distribution.
75  if(col_dist_pt == 0)
76  {
78  (matrix_pt->distribution_pt()->communicator_pt(),
79  matrix_pt->ncol(),matrix_pt->distribution_pt()->distributed());
80  }
81  else
82  {
84  }
85 
86  // setup the operator
87  if (Using_trilinos)
88  {
89 #ifdef OOMPH_HAS_TRILINOS
90  double t_start = TimingHelpers::timer();
93  (matrix_pt,Column_distribution_pt);
94  double t_end = TimingHelpers::timer();
95  oomph_info << "Time to build epetra matrix [sec] : "
96  << t_end - t_start << std::endl;
97 #endif
98  }
99  else
100  {
101  double t_start = TimingHelpers::timer();
102  Oomph_matrix_pt = new CRDoubleMatrix(*matrix_pt);
103  double t_end = TimingHelpers::timer();
104  oomph_info << "Time to copy CRDoubleMatrix [sec] : "
105  << t_end - t_start << std::endl;
106  }
107  }
108 
109  //============================================================================
110  /// \short Apply the operator to the vector x and return the result in
111  /// the vector y
112  //============================================================================
114  DoubleVector& y) const
115  {
116 #ifdef PARANOID
117  // check that the distribution of x is setup
118  if (!x.built())
119  {
120  std::ostringstream error_message_stream;
121  error_message_stream
122  << "The distribution of the vector x must be setup";
123  throw OomphLibError(error_message_stream.str(),
124  OOMPH_CURRENT_FUNCTION,
125  OOMPH_EXCEPTION_LOCATION);
126  }
127 
128  // Check to see if the distribution of the matrix column is the same as the
129  // distribution of the vector to be operated on.
130  if (*this->Column_distribution_pt != *x.distribution_pt())
131  {
132  std::ostringstream error_message_stream;
133  error_message_stream
134  << "The distribution of the x Vector is not the same as"
135  << " the column distribution.";
136  throw OomphLibError(error_message_stream.str(),
137  OOMPH_CURRENT_FUNCTION,
138  OOMPH_EXCEPTION_LOCATION);
139  }
140 
141  // if y is setup then it should have the same distribution as
142  // the matrix used to set this up.
143  if (y.built())
144  {
145  if (!(*y.distribution_pt() == *this->distribution_pt()))
146  {
147  std::ostringstream error_message_stream;
148  error_message_stream
149  << "The y vector is setup and therefore must have the same "
150  << "distribution as the matrix used to set up the MatrixVectorProduct";
151  throw OomphLibError(error_message_stream.str(),
152  OOMPH_CURRENT_FUNCTION,
153  OOMPH_EXCEPTION_LOCATION);
154  }
155  }
156 #endif
157 
158  // if y is not setup then setup the distribution
159  if (!y.built())
160  {
161  // Resize and initialize the solution vector
162  y.build(this->distribution_pt(),0.0);
163  }
164 
165  // apply the operator
166  if (Using_trilinos)
167  {
168 #ifdef OOMPH_HAS_TRILINOS
170 #endif
171  }
172  else
173  {
175  }
176  }
177 
178  //============================================================================
179  /// \short Apply the transpose of the operator to the vector x and return
180  /// the result in the vector y
181  //============================================================================
183  DoubleVector& y) const
184  {
185 #ifdef PARANOID
186  // check that the distribution of x is setup
187  if (!x.built())
188  {
189  std::ostringstream error_message_stream;
190  error_message_stream
191  << "The distribution of the vector x must be setup";
192  throw OomphLibError(error_message_stream.str(),
193  OOMPH_CURRENT_FUNCTION,
194  OOMPH_EXCEPTION_LOCATION);
195  }
196  // Check to see if x.size() = ncol()
197  if (*this->distribution_pt() != *x.distribution_pt())
198  {
199  std::ostringstream error_message_stream;
200  error_message_stream
201  << "This class assumes that the y vector has a uniform "
202  << "distributed distribution.";
203  throw OomphLibError(error_message_stream.str(),
204  OOMPH_CURRENT_FUNCTION,
205  OOMPH_EXCEPTION_LOCATION);
206  }
207  // if y is setup then it should have the same distribution as x
208  if (y.built())
209  {
210  if (!(*y.distribution_pt() == *this->Column_distribution_pt))
211  {
212  std::ostringstream error_message_stream;
213  error_message_stream
214  << "The y vector is setup and therefore must have the same "
215  << "distribution as the vector x";
216  throw OomphLibError(error_message_stream.str(),
217  OOMPH_CURRENT_FUNCTION,
218  OOMPH_EXCEPTION_LOCATION);
219  }
220  }
221 #endif
222 
223  // if y is not setup then setup the distribution
224  if (!y.built())
225  {
226  // Resize and initialize the solution vector
227  y.build(this->Column_distribution_pt,0.0);
228  }
229 
230  // apply the transpose operator
231  if (Using_trilinos)
232  {
233 #ifdef OOMPH_HAS_TRILINOS
235 #endif
236  }
237  else
238  {
240  }
241  }
242 
243 #ifdef OOMPH_HAS_TRILINOS
244  //============================================================================
245  /// \short Apply the operator to the vector x and return
246  /// the result in the vector y (helper function)
247  //============================================================================
249  DoubleVector& y) const
250  {
251  // convert x to a Trilinos Epetra vector.
252  // x is const so it much be copied.
253  Epetra_Vector* epetra_x_pt =
255 
256  // create Trilinos vector for soln ('viewing' the contents of the oomph-lib
257  // matrix)
258  Epetra_Vector* epetra_soln_pt =
260 
261  // do the multiply
262 #ifdef PARANOID
263  int epetra_error_flag = 0;
264  epetra_error_flag =
265 #endif
266  Epetra_matrix_pt->Multiply(false,*epetra_x_pt,
267  *epetra_soln_pt);
268 
269  // throw error if there is an epetra error
270 #ifdef PARANOID
271  if (epetra_error_flag != 0)
272  {
273  std::ostringstream error_message;
274  error_message
275  << "Epetra Matrix Vector Multiply Error : epetra_error_flag = "
276  << epetra_error_flag;
277  throw OomphLibError(error_message.str(),
278  OOMPH_CURRENT_FUNCTION,
279  OOMPH_EXCEPTION_LOCATION);
280  }
281 #endif
282 
283  // return solution
285 
286  // clean up
287  delete epetra_x_pt;
288  delete epetra_soln_pt;
289  }
290 
291  //============================================================================
292  /// \short Apply the transpose of the operator to the vector x and return
293  /// the result in the vector y (helper function)
294  //============================================================================
296  (const DoubleVector& x,DoubleVector& y) const
297  {
298  // convert x to a Trilinos Epetra vector.
299  // x is const so it much be copied.
300  Epetra_Vector* epetra_x_pt =
302 
303  // create Trilinos vector for soln ('viewing' the contents of the oomph-lib
304  // matrix)
305  Epetra_Vector* epetra_soln_pt =
307 
308  // do the multiply
309 #ifdef PARANOID
310  int epetra_error_flag = 0;
311  epetra_error_flag =
312 #endif
313  Epetra_matrix_pt->Multiply(true,*epetra_x_pt,
314  *epetra_soln_pt);
315 
316  // throw error if there is an epetra error
317 #ifdef PARANOID
318  if (epetra_error_flag != 0)
319  {
320  std::ostringstream error_message;
321  error_message
322  << "Epetra Matrix Vector Multiply Error : epetra_error_flag = "
323  << epetra_error_flag;
324  throw OomphLibError(error_message.str(),
325  OOMPH_CURRENT_FUNCTION,
326  OOMPH_EXCEPTION_LOCATION);
327  }
328 #endif
329 
330  // copy to solution vector
332 
333  // clean up
334  delete epetra_x_pt;
335  delete epetra_soln_pt;
336  }
337 #endif
338 }
bool built() const
unsigned Ncol
number of columns of the matrix
OomphInfo oomph_info
LinearAlgebraDistribution * Column_distribution_pt
The distribution of: x if using multiply(...) or y if using multiply_transpose(...) where this is A x = y.
static bool mpi_has_been_initialised()
return true if MPI has been initialised
bool distributed() const
access function to the distributed - indicates whether the distribution is serial or distributed ...
Epetra_CrsMatrix * create_distributed_epetra_matrix(const CRDoubleMatrix *oomph_matrix_pt, const LinearAlgebraDistribution *dist_pt)
create an Epetra_CrsMatrix from an oomph-lib CRDoubleMatrix. If oomph_matrix_pt is NOT distributed (i...
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
void trilinos_multiply_helper(const DoubleVector &x, DoubleVector &y) const
Helper function for multiply(...)
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:998
void trilinos_multiply_transpose_helper(const DoubleVector &x, DoubleVector &y) const
Helper function for multiply_transpose(...)
void multiply(const DoubleVector &x, DoubleVector &y) const
Apply the operator to the vector x and return the result in the vector y.
Epetra_CrsMatrix * Epetra_matrix_pt
The Epetra version of the matrix.
void multiply_transpose(const DoubleVector &x, DoubleVector &soln) const
Multiply the transposed matrix by the vector x: soln=A^T x.
Definition: matrices.cc:1907
void setup(CRDoubleMatrix *matrix_pt, const LinearAlgebraDistribution *col_dist_pt=0)
Setup the matrix vector product operator. WARNING: This class is wrapper to Trilinos Epetra matrix ve...
double timer()
returns the time in seconds after some point in past
bool Using_trilinos
boolean indicating whether we are using trilinos to perform matvec
Epetra_Vector * create_distributed_epetra_vector(const DoubleVector &oomph_vec)
create an Epetra_Vector from an oomph-lib DoubleVector. If oomph_vec is NOT distributed (i...
CRDoubleMatrix * Oomph_matrix_pt
an oomph-lib matrix
void multiply(const DoubleVector &x, DoubleVector &soln) const
Multiply the matrix by the vector x: soln=Ax.
Definition: matrices.cc:1805
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
void multiply_transpose(const DoubleVector &x, DoubleVector &y) const
Apply the transpose of the operator to the vector x and return the result in the vector y...
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
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
void clean_up_memory()
clear the memory
void copy_to_oomphlib_vector(const Epetra_Vector *epetra_vec_pt, DoubleVector &oomph_vec)
Helper function to copy the contents of a Trilinos vector to an oomph-lib distributed vector...
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:872