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