eigen_sparse_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 <algorithm> // for std::min
22 #include <limits>
23 
24 // Local Includes
26 #include "libmesh/dense_vector.h"
29 
30 
31 #ifdef LIBMESH_HAVE_EIGEN
32 
33 namespace libMesh
34 {
35 
36 template <typename T>
38 {
39  libmesh_assert (this->closed());
40  libmesh_assert (this->initialized());
41 
42  return _vec.sum();
43 }
44 
45 
46 
47 template <typename T>
49 {
50  libmesh_assert (this->closed());
51  libmesh_assert (this->initialized());
52 
53  return _vec.lpNorm<1>();
54 }
55 
56 
57 
58 template <typename T>
60 {
61  libmesh_assert (this->closed());
62  libmesh_assert (this->initialized());
63 
64  return _vec.lpNorm<2>();
65 }
66 
67 
68 
69 template <typename T>
71 {
72  libmesh_assert (this->closed());
73  libmesh_assert (this->initialized());
74 
75  return _vec.lpNorm<Eigen::Infinity>();
76 }
77 
78 
79 
80 template <typename T>
82 {
83  libmesh_assert (this->closed());
84 
85  const EigenSparseVector<T> & v = cast_ref<const EigenSparseVector<T> &>(v_in);
86 
87  _vec += v._vec;
88 
89  return *this;
90 }
91 
92 
93 
94 
95 template <typename T>
97 {
98  libmesh_assert (this->closed());
99 
100  const EigenSparseVector<T> & v = cast_ref<const EigenSparseVector<T> &>(v_in);
101 
102  _vec -= v._vec;
103 
104  return *this;
105 }
106 
107 
108 
109 template <typename T>
111 {
112  libmesh_assert (this->closed());
113  libmesh_assert_equal_to(size(), v_in.size());
114 
115  const EigenSparseVector<T> & v = cast_ref<const EigenSparseVector<T> &>(v_in);
116 
117  _vec = _vec.cwiseQuotient(v._vec);
118 
119  return *this;
120 }
121 
122 
123 
124 
125 template <typename T>
127 {
128 #ifndef NDEBUG
129  const numeric_index_type n = this->size();
130 
131  for (numeric_index_type i=0; i<n; i++)
132  // Don't divide by zero!
133  libmesh_assert_not_equal_to ((*this)(i), T(0));
134 #endif
135 
136  _vec = _vec.cwiseInverse();
137 }
138 
139 
140 
141 template <typename T>
143 {
144  _vec = _vec.conjugate();
145 }
146 
147 
148 
149 template <typename T>
151 {
152  _vec += EigenSV::Constant(this->size(), v);
153 
154 #ifndef NDEBUG
155  this->_is_closed = false;
156 #endif
157 }
158 
159 
160 
161 
162 template <typename T>
164 {
165  libmesh_assert (this->initialized());
166 
167  const EigenSparseVector<T> & v = cast_ref<const EigenSparseVector<T> &>(v_in);
168 
169  _vec += v._vec;
170 }
171 
172 
173 
174 template <typename T>
175 void EigenSparseVector<T>::add (const T a, const NumericVector<T> & v_in)
176 {
177  libmesh_assert (this->initialized());
178 
179  const EigenSparseVector<T> & v = cast_ref<const EigenSparseVector<T> &>(v_in);
180 
181  _vec += v._vec*a;
182 }
183 
184 
185 
186 template <typename T>
188  const SparseMatrix<T> & mat_in)
189 {
190  // Make sure the data passed in are really in Eigen types
191  const EigenSparseVector<T> * e_vec = cast_ptr<const EigenSparseVector<T> *>(&vec_in);
192  const EigenSparseMatrix<T> * mat = cast_ptr<const EigenSparseMatrix<T> *>(&mat_in);
193 
194  libmesh_assert(e_vec);
195  libmesh_assert(mat);
196 
197  _vec += mat->_mat*e_vec->_vec;
198 }
199 
200 
201 
202 template <typename T>
204  const SparseMatrix<T> & mat_in)
205 {
206  // Make sure the data passed in are really in Eigen types
207  const EigenSparseVector<T> * e_vec = cast_ptr<const EigenSparseVector<T> *>(&vec_in);
208  const EigenSparseMatrix<T> * mat = cast_ptr<const EigenSparseMatrix<T> *>(&mat_in);
209 
210  libmesh_assert(e_vec);
211  libmesh_assert(mat);
212 
213  _vec += mat->_mat.transpose()*e_vec->_vec;
214 }
215 
216 
217 
218 template <typename T>
219 void EigenSparseVector<T>::scale (const T factor)
220 {
221  libmesh_assert (this->initialized());
222 
223  _vec *= factor;
224 }
225 
226 
227 
228 template <typename T>
230 {
231  libmesh_assert (this->initialized());
232 
233  const numeric_index_type n = this->size();
234 
235  for (numeric_index_type i=0; i!=n; ++i)
236  this->set(i,std::abs((*this)(i)));
237 }
238 
239 
240 
241 template <typename T>
243 {
244  libmesh_assert (this->initialized());
245 
246  // Make sure the NumericVector passed in is really a EigenSparseVector
247  const EigenSparseVector<T> * v = cast_ptr<const EigenSparseVector<T> *>(&V);
248  libmesh_assert(v);
249 
250  return _vec.dot(v->_vec);
251 }
252 
253 
254 
255 template <typename T>
258 {
259  libmesh_assert (this->initialized());
260  libmesh_assert (this->closed());
261 
262  _vec.fill(s);
263 
264  return *this;
265 }
266 
267 
268 
269 template <typename T>
272 {
273  // Make sure the NumericVector passed in is really a EigenSparseVector
274  const EigenSparseVector<T> * v =
275  cast_ptr<const EigenSparseVector<T> *>(&v_in);
276 
277  libmesh_assert(v);
278 
279  *this = *v;
280 
281  return *this;
282 }
283 
284 
285 
286 template <typename T>
289 {
290  libmesh_assert (this->initialized());
291  libmesh_assert (v.closed());
292  libmesh_assert_equal_to (this->size(), v.size());
293 
294  _vec = v._vec;
295 
296 #ifndef NDEBUG
297  this->_is_closed = true;
298 #endif
299 
300  return *this;
301 }
302 
303 
304 
305 template <typename T>
307 EigenSparseVector<T>::operator = (const std::vector<T> & v)
308 {
313  if (this->size() == v.size())
314  for (numeric_index_type i=0; i<v.size(); i++)
315  this->set (i, v[i]);
316 
317  else
318  libmesh_error_msg("this->size() = " << this->size() << " must be equal to v.size() = " << v.size());
319 
320  return *this;
321 }
322 
323 
324 template <typename T>
326 {
327  // Make sure the NumericVector passed in is really a EigenSparseVector
328  EigenSparseVector<T> * v_local =
329  cast_ptr<EigenSparseVector<T> *>(&v_local_in);
330 
331  libmesh_assert(v_local);
332 
333  *v_local = *this;
334 }
335 
336 
337 
338 template <typename T>
340  const std::vector<numeric_index_type> & libmesh_dbg_var(send_list)) const
341 {
342  // Make sure the NumericVector passed in is really a EigenSparseVector
343  EigenSparseVector<T> * v_local =
344  cast_ptr<EigenSparseVector<T> *>(&v_local_in);
345 
346  libmesh_assert(v_local);
347  libmesh_assert_less_equal (send_list.size(), v_local->size());
348 
349  *v_local = *this;
350 }
351 
352 
353 
354 template <typename T>
355 void EigenSparseVector<T>::localize (std::vector<T> & v_local,
356  const std::vector<numeric_index_type> & indices) const
357 {
358  // EigenSparseVectors are serial, so we can just copy values
359  v_local.resize(indices.size());
360 
361  for (numeric_index_type i=0; i<v_local.size(); i++)
362  v_local[i] = (*this)(indices[i]);
363 }
364 
365 
366 
367 template <typename T>
369  const numeric_index_type libmesh_dbg_var(last_local_idx),
370  const std::vector<numeric_index_type> & libmesh_dbg_var(send_list))
371 {
372  libmesh_assert_equal_to (first_local_idx, 0);
373  libmesh_assert_equal_to (last_local_idx+1, this->size());
374 
375  libmesh_assert_less_equal (send_list.size(), this->size());
376 
377 #ifndef NDEBUG
378  this->_is_closed = true;
379 #endif
380 }
381 
382 
383 
384 template <typename T>
385 void EigenSparseVector<T>::localize (std::vector<T> & v_local) const
386 
387 {
388  v_local.resize(this->size());
389 
390  for (numeric_index_type i=0; i<v_local.size(); i++)
391  v_local[i] = (*this)(i);
392 }
393 
394 
395 
396 template <typename T>
397 void EigenSparseVector<T>::localize_to_one (std::vector<T> & v_local,
398  const processor_id_type libmesh_dbg_var(pid)) const
399 {
400  libmesh_assert_equal_to (pid, 0);
401 
402  this->localize (v_local);
403 }
404 
405 
406 
407 template <typename T>
409  const NumericVector<T> & /*vec2*/)
410 {
411  libmesh_not_implemented();
412 }
413 
414 
415 
416 template <typename T>
418 {
419  libmesh_assert (this->initialized());
420  if (!this->size())
422 
423 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
424  Real the_max = libmesh_real((*this)(0));
425 
426  const numeric_index_type n = this->size();
427 
428  for (numeric_index_type i=1; i<n; i++)
429  the_max = std::max (the_max, libmesh_real((*this)(i)));
430 
431  return the_max;
432 #else
433  return libmesh_real(_vec.maxCoeff());
434 #endif
435 }
436 
437 
438 
439 template <typename T>
441 {
442  libmesh_assert (this->initialized());
443  if (!this->size())
445 
446 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
447  Real the_min = libmesh_real((*this)(0));
448 
449  const numeric_index_type n = this->size();
450 
451  for (numeric_index_type i=1; i<n; i++)
452  the_min = std::min (the_min, libmesh_real((*this)(i)));
453 
454  return the_min;
455 #else
456  return libmesh_real(_vec.minCoeff());
457 #endif
458 }
459 
460 
461 //------------------------------------------------------------------
462 // Explicit instantiations
463 template class EigenSparseVector<Number>;
464 
465 } // namespace libMesh
466 
467 
468 #endif // #ifdef LIBMESH_HAVE_EIGEN
T libmesh_real(T a)
virtual NumericVector< T > & operator=(const T s) libmesh_override
virtual bool closed() const
bool closed()
Definition: libmesh.C:279
double abs(double a)
virtual void add_vector_transpose(const NumericVector< T > &, const SparseMatrix< T > &) libmesh_override
virtual numeric_index_type size() const =0
virtual void localize_to_one(std::vector< T > &v_local, const processor_id_type proc_id=0) const libmesh_override
virtual NumericVector< T > & operator-=(const NumericVector< T > &v) libmesh_override
virtual Real min() const libmesh_override
virtual void pointwise_mult(const NumericVector< T > &vec1, const NumericVector< T > &vec2) libmesh_override
uint8_t processor_id_type
Definition: id_types.h:99
virtual void reciprocal() libmesh_override
long double max(long double a, double b)
virtual NumericVector< T > & operator/=(NumericVector< T > &v_in) libmesh_override
libmesh_assert(j)
virtual void scale(const T factor) libmesh_override
virtual void localize(std::vector< T > &v_local) const libmesh_override
virtual NumericVector< T > & operator+=(const NumericVector< T > &v) libmesh_override
dof_id_type numeric_index_type
Definition: id_types.h:92
virtual void add(const numeric_index_type i, const T value) libmesh_override
PetscErrorCode Vec Mat libmesh_dbg_var(j)
virtual void conjugate() libmesh_override
virtual void abs() libmesh_override
virtual T sum() const libmesh_override
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual numeric_index_type size() const libmesh_override
virtual T dot(const NumericVector< T > &V) const libmesh_override
bool initialized()
Definition: libmesh.C:272
virtual Real max() const libmesh_override
virtual Real linfty_norm() const libmesh_override
long double min(long double a, double b)
virtual Real l2_norm() const libmesh_override
virtual void add_vector(const NumericVector< T > &, const SparseMatrix< T > &) libmesh_override
virtual Real l1_norm() const libmesh_override