libMesh
numeric_vector.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // C++ includes
21 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
22 #include <cmath> // for std::abs
23 #include <limits>
24 
25 // Local Includes
26 #include "libmesh/numeric_vector.h"
27 #include "libmesh/distributed_vector.h"
28 #include "libmesh/laspack_vector.h"
29 #include "libmesh/eigen_sparse_vector.h"
30 #include "libmesh/petsc_vector.h"
31 #include "libmesh/trilinos_epetra_vector.h"
32 #include "libmesh/shell_matrix.h"
33 #include "libmesh/tensor_tools.h"
34 
35 namespace libMesh
36 {
37 
38 
39 
40 //------------------------------------------------------------------
41 // NumericVector methods
42 
43 // Full specialization for Real datatypes
44 template <typename T>
45 UniquePtr<NumericVector<T>>
47 {
48  // Build the appropriate vector
49  switch (solver_package)
50  {
51 
52 #ifdef LIBMESH_HAVE_LASPACK
53  case LASPACK_SOLVERS:
55 #endif
56 
57 #ifdef LIBMESH_HAVE_PETSC
58  case PETSC_SOLVERS:
59  return UniquePtr<NumericVector<T>>(new PetscVector<T>(comm, AUTOMATIC));
60 #endif
61 
62 #ifdef LIBMESH_TRILINOS_HAVE_EPETRA
63  case TRILINOS_SOLVERS:
64  return UniquePtr<NumericVector<T>>(new EpetraVector<T>(comm, AUTOMATIC));
65 #endif
66 
67 #ifdef LIBMESH_HAVE_EIGEN
68  case EIGEN_SOLVERS:
70 #endif
71 
72  default:
73  return UniquePtr<NumericVector<T>>(new DistributedVector<T>(comm, AUTOMATIC));
74  }
75 
76  libmesh_error_msg("We'll never get here!");
78 }
79 
80 
81 #ifndef LIBMESH_DISABLE_COMMWORLD
82 #ifdef LIBMESH_ENABLE_DEPRECATED
83 template <typename T>
86 {
87  libmesh_deprecated();
88  return NumericVector<T>::build(CommWorld, solver_package);
89 }
90 #endif
91 #endif
92 
93 
94 
95 template <typename T>
96 void NumericVector<T>::insert (const T * v,
97  const std::vector<numeric_index_type> & dof_indices)
98 {
99  for (numeric_index_type i=0; i<dof_indices.size(); i++)
100  this->set (dof_indices[i], v[i]);
101 }
102 
103 
104 
105 template <typename T>
107  const std::vector<numeric_index_type> & dof_indices)
108 {
109  libmesh_assert_equal_to (V.size(), dof_indices.size());
110 
111  for (numeric_index_type i=0; i<dof_indices.size(); i++)
112  this->set (dof_indices[i], V(i));
113 }
114 
115 
116 
117 template <typename T>
118 int NumericVector<T>::compare (const NumericVector<T> & other_vector,
119  const Real threshold) const
120 {
121  libmesh_assert (this->initialized());
122  libmesh_assert (other_vector.initialized());
123  libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
124  libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
125 
126  int first_different_i = std::numeric_limits<int>::max();
127  numeric_index_type i = first_local_index();
128 
129  do
130  {
131  if (std::abs((*this)(i) - other_vector(i)) > threshold)
132  first_different_i = i;
133  else
134  i++;
135  }
136  while (first_different_i==std::numeric_limits<int>::max()
137  && i<last_local_index());
138 
139  // Find the correct first differing index in parallel
140  this->comm().min(first_different_i);
141 
142  if (first_different_i == std::numeric_limits<int>::max())
143  return -1;
144 
145  return first_different_i;
146 }
147 
148 
149 template <typename T>
151  const Real threshold) const
152 {
153  libmesh_assert (this->initialized());
154  libmesh_assert (other_vector.initialized());
155  libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
156  libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
157 
158  int first_different_i = std::numeric_limits<int>::max();
159  numeric_index_type i = first_local_index();
160 
161  do
162  {
163  if (std::abs((*this)(i) - other_vector(i)) > threshold *
164  std::max(std::abs((*this)(i)), std::abs(other_vector(i))))
165  first_different_i = i;
166  else
167  i++;
168  }
169  while (first_different_i==std::numeric_limits<int>::max()
170  && i<last_local_index());
171 
172  // Find the correct first differing index in parallel
173  this->comm().min(first_different_i);
174 
175  if (first_different_i == std::numeric_limits<int>::max())
176  return -1;
177 
178  return first_different_i;
179 }
180 
181 
182 template <typename T>
184  const Real threshold) const
185 {
186  libmesh_assert (this->initialized());
187  libmesh_assert (other_vector.initialized());
188  libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
189  libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
190 
191  int first_different_i = std::numeric_limits<int>::max();
192  numeric_index_type i = first_local_index();
193 
194  const Real my_norm = this->linfty_norm();
195  const Real other_norm = other_vector.linfty_norm();
196  const Real abs_threshold = std::max(my_norm, other_norm) * threshold;
197 
198  do
199  {
200  if (std::abs((*this)(i) - other_vector(i) ) > abs_threshold)
201  first_different_i = i;
202  else
203  i++;
204  }
205  while (first_different_i==std::numeric_limits<int>::max()
206  && i<last_local_index());
207 
208  // Find the correct first differing index in parallel
209  this->comm().min(first_different_i);
210 
211  if (first_different_i == std::numeric_limits<int>::max())
212  return -1;
213 
214  return first_different_i;
215 }
216 
217 /*
218 // Full specialization for float datatypes (DistributedVector wants this)
219 
220 template <>
221 int NumericVector<float>::compare (const NumericVector<float> & other_vector,
222 const Real threshold) const
223 {
224 libmesh_assert (this->initialized());
225 libmesh_assert (other_vector.initialized());
226 libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
227 libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
228 
229 int rvalue = -1;
230 numeric_index_type i = first_local_index();
231 
232 do
233 {
234 if (std::abs((*this)(i) - other_vector(i) ) > threshold)
235 rvalue = i;
236 else
237 i++;
238 }
239 while (rvalue==-1 && i<last_local_index());
240 
241 return rvalue;
242 }
243 
244 // Full specialization for double datatypes
245 template <>
246 int NumericVector<double>::compare (const NumericVector<double> & other_vector,
247 const Real threshold) const
248 {
249 libmesh_assert (this->initialized());
250 libmesh_assert (other_vector.initialized());
251 libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
252 libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
253 
254 int rvalue = -1;
255 numeric_index_type i = first_local_index();
256 
257 do
258 {
259 if (std::abs((*this)(i) - other_vector(i) ) > threshold)
260 rvalue = i;
261 else
262 i++;
263 }
264 while (rvalue==-1 && i<last_local_index());
265 
266 return rvalue;
267 }
268 
269 #ifdef LIBMESH_DEFAULT_TRIPLE_PRECISION
270 // Full specialization for long double datatypes
271 template <>
272 int NumericVector<long double>::compare (const NumericVector<long double> & other_vector,
273 const Real threshold) const
274 {
275 libmesh_assert (this->initialized());
276 libmesh_assert (other_vector.initialized());
277 libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
278 libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
279 
280 int rvalue = -1;
281 numeric_index_type i = first_local_index();
282 
283 do
284 {
285 if (std::abs((*this)(i) - other_vector(i) ) > threshold)
286 rvalue = i;
287 else
288 i++;
289 }
290 while (rvalue==-1 && i<last_local_index());
291 
292 return rvalue;
293 }
294 #endif
295 
296 
297 // Full specialization for Complex datatypes
298 template <>
299 int NumericVector<Complex>::compare (const NumericVector<Complex> & other_vector,
300 const Real threshold) const
301 {
302 libmesh_assert (this->initialized());
303 libmesh_assert (other_vector.initialized());
304 libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
305 libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
306 
307 int rvalue = -1;
308 numeric_index_type i = first_local_index();
309 
310 do
311 {
312 if ((std::abs((*this)(i).real() - other_vector(i).real()) > threshold) || (std::abs((*this)(i).imag() - other_vector(i).imag()) > threshold))
313 rvalue = i;
314 else
315 i++;
316 }
317 while (rvalue==-1 && i<this->last_local_index());
318 
319 return rvalue;
320 }
321 */
322 
323 
324 template <class T>
325 Real NumericVector<T>::subset_l1_norm (const std::set<numeric_index_type> & indices) const
326 {
327  const NumericVector<T> & v = *this;
328 
329  std::set<numeric_index_type>::const_iterator it = indices.begin();
330  const std::set<numeric_index_type>::const_iterator it_end = indices.end();
331 
332  Real norm = 0;
333 
334  for (; it!=it_end; ++it)
335  norm += std::abs(v(*it));
336 
337  this->comm().sum(norm);
338 
339  return norm;
340 }
341 
342 template <class T>
343 Real NumericVector<T>::subset_l2_norm (const std::set<numeric_index_type> & indices) const
344 {
345  const NumericVector<T> & v = *this;
346 
347  std::set<numeric_index_type>::const_iterator it = indices.begin();
348  const std::set<numeric_index_type>::const_iterator it_end = indices.end();
349 
350  Real norm = 0;
351 
352  for (; it!=it_end; ++it)
353  norm += TensorTools::norm_sq(v(*it));
354 
355  this->comm().sum(norm);
356 
357  return std::sqrt(norm);
358 }
359 
360 template <class T>
361 Real NumericVector<T>::subset_linfty_norm (const std::set<numeric_index_type> & indices) const
362 {
363  const NumericVector<T> & v = *this;
364 
365  std::set<numeric_index_type>::const_iterator it = indices.begin();
366  const std::set<numeric_index_type>::const_iterator it_end = indices.end();
367 
368  Real norm = 0;
369 
370  for (; it!=it_end; ++it)
371  {
372  Real value = std::abs(v(*it));
373  if (value > norm)
374  norm = value;
375  }
376 
377  this->comm().max(norm);
378 
379  return norm;
380 }
381 
382 
383 
384 template <typename T>
386  const std::vector<numeric_index_type> & dof_indices)
387 {
388  int n = dof_indices.size();
389  for (int i=0; i<n; i++)
390  this->add (dof_indices[i], v[i]);
391 }
392 
393 
394 
395 template <typename T>
397  const std::vector<numeric_index_type> & dof_indices)
398 {
399  int n = dof_indices.size();
400  libmesh_assert_equal_to(v.size(), static_cast<unsigned>(n));
401  for (int i=0; i<n; i++)
402  this->add (dof_indices[i], v(i));
403 }
404 
405 
406 
407 template <typename T>
409  const ShellMatrix<T> & a)
410 {
411  a.vector_mult_add(*this,v);
412 }
413 
414 
415 
416 //------------------------------------------------------------------
417 // Explicit instantiations
418 template class NumericVector<Number>;
419 
420 } // namespace libMesh
virtual void insert(const T *v, const std::vector< numeric_index_type > &dof_indices)
Inserts the entries of v in *this at the locations specified by v.
double abs(double a)
Encapsulates the MPI_Comm object.
Definition: parallel.h:657
virtual int compare(const NumericVector< T > &other_vector, const Real threshold=TOLERANCE) const
virtual numeric_index_type size() const =0
virtual numeric_index_type last_local_index() const =0
virtual Real subset_linfty_norm(const std::set< numeric_index_type > &indices) const
virtual int global_relative_compare(const NumericVector< T > &other_vector, const Real threshold=TOLERANCE) const
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
SolverPackage
Defines an enum for various linear solver packages.
Numeric vector.
Definition: dof_map.h:66
virtual bool initialized() const
The libMesh namespace provides an interface to certain functionality in the library.
long double max(long double a, double b)
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
dof_id_type numeric_index_type
Definition: id_types.h:92
virtual Real subset_l1_norm(const std::set< numeric_index_type > &indices) const
virtual int local_relative_compare(const NumericVector< T > &other_vector, const Real threshold=TOLERANCE) const
virtual Real linfty_norm() const =0
Parallel::FakeCommunicator CommWorld
The default libMesh communicator.
Definition: parallel.h:1433
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static const bool value
Definition: xdr_io.C:108
virtual numeric_index_type first_local_index() const =0
bool initialized()
Checks that library initialization has been done.
Definition: libmesh.C:274
virtual void vector_mult_add(NumericVector< T > &dest, const NumericVector< T > &arg) const =0
Multiplies the matrix with arg and adds the result to dest.
Generic shell matrix, i.e.
virtual Real subset_l2_norm(const std::set< numeric_index_type > &indices) const