numeric_vector.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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"
28 #include "libmesh/laspack_vector.h"
30 #include "libmesh/petsc_vector.h"
32 #include "libmesh/shell_matrix.h"
33 #include "libmesh/tensor_tools.h"
34 #include "libmesh/auto_ptr.h" // libmesh_make_unique
36 
37 namespace libMesh
38 {
39 
40 
41 
42 //------------------------------------------------------------------
43 // NumericVector methods
44 
45 // Full specialization for Real datatypes
46 template <typename T>
47 std::unique_ptr<NumericVector<T>>
49 {
50  // Build the appropriate vector
51  switch (solver_package)
52  {
53 
54 #ifdef LIBMESH_HAVE_LASPACK
55  case LASPACK_SOLVERS:
56  return libmesh_make_unique<LaspackVector<T>>(comm, AUTOMATIC);
57 #endif
58 
59 #ifdef LIBMESH_HAVE_PETSC
60  case PETSC_SOLVERS:
61  return libmesh_make_unique<PetscVector<T>>(comm, AUTOMATIC);
62 #endif
63 
64 #ifdef LIBMESH_TRILINOS_HAVE_EPETRA
65  case TRILINOS_SOLVERS:
66  return libmesh_make_unique<EpetraVector<T>>(comm, AUTOMATIC);
67 #endif
68 
69 #ifdef LIBMESH_HAVE_EIGEN
70  case EIGEN_SOLVERS:
71  return libmesh_make_unique<EigenSparseVector<T>>(comm, AUTOMATIC);
72 #endif
73 
74  default:
75  return libmesh_make_unique<DistributedVector<T>>(comm, AUTOMATIC);
76  }
77 }
78 
79 
80 
81 template <typename T>
82 void NumericVector<T>::insert (const T * v,
83  const std::vector<numeric_index_type> & dof_indices)
84 {
85  for (numeric_index_type i=0; i<dof_indices.size(); i++)
86  this->set (dof_indices[i], v[i]);
87 }
88 
89 
90 
91 template <typename T>
93  const std::vector<numeric_index_type> & dof_indices)
94 {
95  libmesh_assert_equal_to (V.size(), dof_indices.size());
96 
97  for (numeric_index_type i=0; i<dof_indices.size(); i++)
98  this->set (dof_indices[i], V(i));
99 }
100 
101 
102 
103 template <typename T>
104 int NumericVector<T>::compare (const NumericVector<T> & other_vector,
105  const Real threshold) const
106 {
107  libmesh_assert (this->initialized());
108  libmesh_assert (other_vector.initialized());
109  libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
110  libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
111 
112  int first_different_i = std::numeric_limits<int>::max();
113  numeric_index_type i = first_local_index();
114 
115  do
116  {
117  if (std::abs((*this)(i) - other_vector(i)) > threshold)
118  first_different_i = i;
119  else
120  i++;
121  }
122  while (first_different_i==std::numeric_limits<int>::max()
123  && i<last_local_index());
124 
125  // Find the correct first differing index in parallel
126  this->comm().min(first_different_i);
127 
128  if (first_different_i == std::numeric_limits<int>::max())
129  return -1;
130 
131  return first_different_i;
132 }
133 
134 
135 template <typename T>
137  const Real threshold) const
138 {
139  libmesh_assert (this->initialized());
140  libmesh_assert (other_vector.initialized());
141  libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
142  libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
143 
144  int first_different_i = std::numeric_limits<int>::max();
145  numeric_index_type i = first_local_index();
146 
147  do
148  {
149  if (std::abs((*this)(i) - other_vector(i)) > threshold *
150  std::max(std::abs((*this)(i)), std::abs(other_vector(i))))
151  first_different_i = i;
152  else
153  i++;
154  }
155  while (first_different_i==std::numeric_limits<int>::max()
156  && i<last_local_index());
157 
158  // Find the correct first differing index in parallel
159  this->comm().min(first_different_i);
160 
161  if (first_different_i == std::numeric_limits<int>::max())
162  return -1;
163 
164  return first_different_i;
165 }
166 
167 
168 template <typename T>
170  const Real threshold) const
171 {
172  libmesh_assert (this->initialized());
173  libmesh_assert (other_vector.initialized());
174  libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
175  libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
176 
177  int first_different_i = std::numeric_limits<int>::max();
178  numeric_index_type i = first_local_index();
179 
180  const Real my_norm = this->linfty_norm();
181  const Real other_norm = other_vector.linfty_norm();
182  const Real abs_threshold = std::max(my_norm, other_norm) * threshold;
183 
184  do
185  {
186  if (std::abs((*this)(i) - other_vector(i) ) > abs_threshold)
187  first_different_i = i;
188  else
189  i++;
190  }
191  while (first_different_i==std::numeric_limits<int>::max()
192  && i<last_local_index());
193 
194  // Find the correct first differing index in parallel
195  this->comm().min(first_different_i);
196 
197  if (first_different_i == std::numeric_limits<int>::max())
198  return -1;
199 
200  return first_different_i;
201 }
202 
203 /*
204 // Full specialization for float datatypes (DistributedVector wants this)
205 
206 template <>
207 int NumericVector<float>::compare (const NumericVector<float> & other_vector,
208 const Real threshold) const
209 {
210 libmesh_assert (this->initialized());
211 libmesh_assert (other_vector.initialized());
212 libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
213 libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
214 
215 int rvalue = -1;
216 numeric_index_type i = first_local_index();
217 
218 do
219 {
220 if (std::abs((*this)(i) - other_vector(i) ) > threshold)
221 rvalue = i;
222 else
223 i++;
224 }
225 while (rvalue==-1 && i<last_local_index());
226 
227 return rvalue;
228 }
229 
230 // Full specialization for double datatypes
231 template <>
232 int NumericVector<double>::compare (const NumericVector<double> & other_vector,
233 const Real threshold) const
234 {
235 libmesh_assert (this->initialized());
236 libmesh_assert (other_vector.initialized());
237 libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
238 libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
239 
240 int rvalue = -1;
241 numeric_index_type i = first_local_index();
242 
243 do
244 {
245 if (std::abs((*this)(i) - other_vector(i) ) > threshold)
246 rvalue = i;
247 else
248 i++;
249 }
250 while (rvalue==-1 && i<last_local_index());
251 
252 return rvalue;
253 }
254 
255 #ifdef LIBMESH_DEFAULT_TRIPLE_PRECISION
256 // Full specialization for long double datatypes
257 template <>
258 int NumericVector<long double>::compare (const NumericVector<long double> & other_vector,
259 const Real threshold) const
260 {
261 libmesh_assert (this->initialized());
262 libmesh_assert (other_vector.initialized());
263 libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
264 libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
265 
266 int rvalue = -1;
267 numeric_index_type i = first_local_index();
268 
269 do
270 {
271 if (std::abs((*this)(i) - other_vector(i) ) > threshold)
272 rvalue = i;
273 else
274 i++;
275 }
276 while (rvalue==-1 && i<last_local_index());
277 
278 return rvalue;
279 }
280 #endif
281 
282 
283 // Full specialization for Complex datatypes
284 template <>
285 int NumericVector<Complex>::compare (const NumericVector<Complex> & other_vector,
286 const Real threshold) const
287 {
288 libmesh_assert (this->initialized());
289 libmesh_assert (other_vector.initialized());
290 libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
291 libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
292 
293 int rvalue = -1;
294 numeric_index_type i = first_local_index();
295 
296 do
297 {
298 if ((std::abs((*this)(i).real() - other_vector(i).real()) > threshold) || (std::abs((*this)(i).imag() - other_vector(i).imag()) > threshold))
299 rvalue = i;
300 else
301 i++;
302 }
303 while (rvalue==-1 && i<this->last_local_index());
304 
305 return rvalue;
306 }
307 */
308 
309 
310 template <class T>
311 Real NumericVector<T>::subset_l1_norm (const std::set<numeric_index_type> & indices) const
312 {
313  const NumericVector<T> & v = *this;
314 
315  Real norm = 0;
316 
317  for (const auto & index : indices)
318  norm += std::abs(v(index));
319 
320  this->comm().sum(norm);
321 
322  return norm;
323 }
324 
325 template <class T>
326 Real NumericVector<T>::subset_l2_norm (const std::set<numeric_index_type> & indices) const
327 {
328  const NumericVector<T> & v = *this;
329 
330  Real norm = 0;
331 
332  for (const auto & index : indices)
333  norm += TensorTools::norm_sq(v(index));
334 
335  this->comm().sum(norm);
336 
337  return std::sqrt(norm);
338 }
339 
340 template <class T>
341 Real NumericVector<T>::subset_linfty_norm (const std::set<numeric_index_type> & indices) const
342 {
343  const NumericVector<T> & v = *this;
344 
345  Real norm = 0;
346 
347  for (const auto & index : indices)
348  {
349  Real value = std::abs(v(index));
350  if (value > norm)
351  norm = value;
352  }
353 
354  this->comm().max(norm);
355 
356  return norm;
357 }
358 
359 
360 
361 template <typename T>
363  const std::vector<numeric_index_type> & dof_indices)
364 {
365  for (std::size_t i=0, n = dof_indices.size(); i != n; i++)
366  this->add (dof_indices[i], v[i]);
367 }
368 
369 
370 
371 template <typename T>
373  const std::vector<numeric_index_type> & dof_indices)
374 {
375  const std::size_t n = dof_indices.size();
376  libmesh_assert_equal_to(v.size(), n);
377  for (numeric_index_type i=0; i != n; i++)
378  this->add (dof_indices[i], v(i));
379 }
380 
381 
382 
383 template <typename T>
385  const ShellMatrix<T> & a)
386 {
387  a.vector_mult_add(*this,v);
388 }
389 
390 
391 
392 //------------------------------------------------------------------
393 // Explicit instantiations
394 template class NumericVector<Number>;
395 
396 } // namespace libMesh
virtual void insert(const T *v, const std::vector< numeric_index_type > &dof_indices)
double abs(double a)
virtual Real subset_l2_norm(const std::set< numeric_index_type > &indices) const
EIGEN_SOLVERS
Definition: libmesh.C:246
TRILINOS_SOLVERS
Definition: libmesh.C:244
virtual bool initialized() const
virtual numeric_index_type size() const =0
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: diff_context.h:40
virtual Real subset_linfty_norm(const std::set< numeric_index_type > &indices) const
long double max(long double a, double b)
virtual Real subset_l1_norm(const std::set< numeric_index_type > &indices) const
virtual void vector_mult_add(NumericVector< T > &dest, const NumericVector< T > &arg) const =0
dof_id_type numeric_index_type
Definition: id_types.h:92
LASPACK_SOLVERS
Definition: libmesh.C:248
virtual int global_relative_compare(const NumericVector< T > &other_vector, const Real threshold=TOLERANCE) const
static std::unique_ptr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
virtual numeric_index_type first_local_index() const =0
virtual int local_relative_compare(const NumericVector< T > &other_vector, const Real threshold=TOLERANCE) const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual int compare(const NumericVector< T > &other_vector, const Real threshold=TOLERANCE) const
static const bool value
Definition: xdr_io.C:109
bool initialized()
Definition: libmesh.C:258
virtual numeric_index_type last_local_index() const =0
virtual Real linfty_norm() const =0