libMesh
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
25 #include "libmesh/dense_subvector.h"
26 #include "libmesh/dense_vector.h"
27 #include "libmesh/eigen_sparse_vector.h"
28 #include "libmesh/eigen_sparse_matrix.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>
48 Real EigenSparseVector<T>::l1_norm () const
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>
59 Real EigenSparseVector<T>::l2_norm () const
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>
70 Real EigenSparseVector<T>::linfty_norm () const
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>
81 NumericVector<T> & EigenSparseVector<T>::operator += (const NumericVector<T> & v_in)
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>
96 NumericVector<T> & EigenSparseVector<T>::operator -= (const NumericVector<T> & v_in)
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>
110 NumericVector<T> & EigenSparseVector<T>::operator /= (NumericVector<T> & v_in)
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>
126 void EigenSparseVector<T>::reciprocal()
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>
142 void EigenSparseVector<T>::conjugate()
143 {
144  _vec = _vec.conjugate();
145 }
146 
147 
148 
149 template <typename T>
150 void EigenSparseVector<T>::add (const T v)
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>
163 void EigenSparseVector<T>::add (const NumericVector<T> & v_in)
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>
187 void EigenSparseVector<T>::add_vector (const NumericVector<T> & vec_in,
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>
203 void EigenSparseVector<T>::add_vector_transpose (const NumericVector<T> & vec_in,
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>
242 T EigenSparseVector<T>::dot (const NumericVector<T> & v_in) const
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_in);
248  libmesh_assert(v);
249 
250  return _vec.dot(v->_vec);
251 }
252 
253 
254 
255 template <typename T>
256 NumericVector<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>
270 NumericVector<T> &
271 EigenSparseVector<T>::operator = (const NumericVector<T> & v_in)
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>
287 EigenSparseVector<T> &
288 EigenSparseVector<T>::operator = (const EigenSparseVector<T> & v)
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>
306 NumericVector<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>
325 void EigenSparseVector<T>::localize (NumericVector<T> & v_local_in) const
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>
339 void EigenSparseVector<T>::localize (NumericVector<T> & v_local_in,
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>
368 void EigenSparseVector<T>::localize (const numeric_index_type libmesh_dbg_var(first_local_idx),
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>
408 void EigenSparseVector<T>::pointwise_mult (const NumericVector<T> & /*vec1*/,
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)
bool closed()
Checks that the library has been closed.
Definition: libmesh.C:281
double abs(double a)
void scale(MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
Scales the mesh.
uint8_t processor_id_type
Definition: id_types.h:99
The libMesh namespace provides an interface to certain functionality in the library.
void sum(T &r, const Communicator &comm=Communicator_World)
long double max(long double a, double b)
libmesh_assert(j)
dof_id_type numeric_index_type
Definition: id_types.h:92
PetscErrorCode Vec Mat libmesh_dbg_var(j)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Iterator & operator=(const Iterator &rhs)
Assignment operator.
bool initialized()
Checks that library initialization has been done.
Definition: libmesh.C:274
long double min(long double a, double b)