eigen_sparse_vector.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 = libmesh_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 = libmesh_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 = libmesh_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 = libmesh_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 = libmesh_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 std::vector<T>& v,
188  const std::vector<numeric_index_type>& dof_indices)
189 {
190  libmesh_assert (!v.empty());
191  libmesh_assert_equal_to (v.size(), dof_indices.size());
192 
193  for (numeric_index_type i=0; i<v.size(); i++)
194  this->add (dof_indices[i], v[i]);
195 }
196 
197 
198 
199 template <typename T>
201  const std::vector<numeric_index_type>& dof_indices)
202 {
203  libmesh_assert_equal_to (V.size(), dof_indices.size());
204 
205  for (numeric_index_type i=0; i<V.size(); i++)
206  this->add (dof_indices[i], V(i));
207 }
208 
209 
210 
211 template <typename T>
213  const std::vector<numeric_index_type>& dof_indices)
214 {
215  libmesh_assert_equal_to (V.size(), dof_indices.size());
216 
217  for (unsigned int i=0; i<V.size(); i++)
218  this->add (dof_indices[i], V(i));
219 }
220 
221 
222 
223 template <typename T>
224 void EigenSparseVector<T>::insert (const std::vector<T>& v,
225  const std::vector<numeric_index_type>& dof_indices)
226 {
227  libmesh_assert (!v.empty());
228  libmesh_assert_equal_to (v.size(), dof_indices.size());
229 
230  for (numeric_index_type i=0; i<v.size(); i++)
231  this->set (dof_indices[i], v[i]);
232 }
233 
234 
235 
236 template <typename T>
238  const std::vector<numeric_index_type>& dof_indices)
239 {
240  libmesh_assert_equal_to (V.size(), dof_indices.size());
241 
242  for (numeric_index_type i=0; i<V.size(); i++)
243  this->set (dof_indices[i], V(i));
244 }
245 
246 
247 
248 template <typename T>
250  const std::vector<numeric_index_type>& dof_indices)
251 {
252  libmesh_assert_equal_to (V.size(), dof_indices.size());
253 
254  for (unsigned int i=0; i<V.size(); i++)
255  this->set (dof_indices[i], V(i));
256 }
257 
258 
259 
260 template <typename T>
262  const std::vector<numeric_index_type>& dof_indices)
263 {
264  libmesh_assert_equal_to (V.size(), dof_indices.size());
265 
266  for (unsigned int i=0; i<V.size(); i++)
267  this->set (dof_indices[i], V(i));
268 }
269 
270 
271 
272 template <typename T>
274  const SparseMatrix<T> &mat_in)
275 {
276  // Make sure the data passed in are really in Eigen types
277  const EigenSparseVector<T>* vec = libmesh_cast_ptr<const EigenSparseVector<T>*>(&vec_in);
278  const EigenSparseMatrix<T>* mat = libmesh_cast_ptr<const EigenSparseMatrix<T>*>(&mat_in);
279 
280  libmesh_assert(vec);
281  libmesh_assert(mat);
282 
283  _vec += mat->_mat*vec->_vec;
284 }
285 
286 
287 
288 template <typename T>
290  const SparseMatrix<T> &mat_in)
291 {
292  // Make sure the data passed in are really in Eigen types
293  const EigenSparseVector<T>* vec = libmesh_cast_ptr<const EigenSparseVector<T>*>(&vec_in);
294  const EigenSparseMatrix<T>* mat = libmesh_cast_ptr<const EigenSparseMatrix<T>*>(&mat_in);
295 
296  libmesh_assert(vec);
297  libmesh_assert(mat);
298 
299  _vec += mat->_mat.transpose()*vec->_vec;
300 }
301 
302 
303 
304 template <typename T>
305 void EigenSparseVector<T>::scale (const T factor)
306 {
307  libmesh_assert (this->initialized());
308 
309  _vec *= factor;
310 }
311 
312 
313 
314 template <typename T>
316 {
317  libmesh_assert (this->initialized());
318 
319  const numeric_index_type n = this->size();
320 
321  for (numeric_index_type i=0; i!=n; ++i)
322  this->set(i,std::abs((*this)(i)));
323 }
324 
325 
326 
327 template <typename T>
329 {
330  libmesh_assert (this->initialized());
331 
332  // Make sure the NumericVector passed in is really a EigenSparseVector
333  const EigenSparseVector<T>* v = libmesh_cast_ptr<const EigenSparseVector<T>*>(&V);
334  libmesh_assert(v);
335 
336  return _vec.dot(v->_vec);
337 }
338 
339 
340 
341 template <typename T>
344 {
345  libmesh_assert (this->initialized());
346  libmesh_assert (this->closed());
347 
348  _vec.fill(s);
349 
350  return *this;
351 }
352 
353 
354 
355 template <typename T>
358 {
359  // Make sure the NumericVector passed in is really a EigenSparseVector
360  const EigenSparseVector<T>* v =
361  libmesh_cast_ptr<const EigenSparseVector<T>*>(&v_in);
362 
363  libmesh_assert(v);
364 
365  *this = *v;
366 
367  return *this;
368 }
369 
370 
371 
372 template <typename T>
375 {
376  libmesh_assert (this->initialized());
377  libmesh_assert (v.closed());
378  libmesh_assert_equal_to (this->size(), v.size());
379 
380  _vec = v._vec;
381 
382 #ifndef NDEBUG
383  this->_is_closed = true;
384 #endif
385 
386  return *this;
387 }
388 
389 
390 
391 template <typename T>
393 EigenSparseVector<T>::operator = (const std::vector<T>& v)
394 {
399  if (this->size() == v.size())
400  for (numeric_index_type i=0; i<v.size(); i++)
401  this->set (i, v[i]);
402 
403  else
404  libmesh_error();
405 
406  return *this;
407 }
408 
409 
410 template <typename T>
412 {
413  // Make sure the NumericVector passed in is really a EigenSparseVector
414  EigenSparseVector<T>* v_local =
415  libmesh_cast_ptr<EigenSparseVector<T>*>(&v_local_in);
416 
417  libmesh_assert(v_local);
418 
419  *v_local = *this;
420 }
421 
422 
423 
424 template <typename T>
426  const std::vector<numeric_index_type>& libmesh_dbg_var(send_list)) const
427 {
428  // Make sure the NumericVector passed in is really a EigenSparseVector
429  EigenSparseVector<T>* v_local =
430  libmesh_cast_ptr<EigenSparseVector<T>*>(&v_local_in);
431 
432  libmesh_assert(v_local);
433  libmesh_assert_less_equal (send_list.size(), v_local->size());
434 
435  *v_local = *this;
436 }
437 
438 
439 
440 template <typename T>
441 void EigenSparseVector<T>::localize (const numeric_index_type libmesh_dbg_var(first_local_idx),
442  const numeric_index_type libmesh_dbg_var(last_local_idx),
443  const std::vector<numeric_index_type>& libmesh_dbg_var(send_list))
444 {
445  libmesh_assert_equal_to (first_local_idx, 0);
446  libmesh_assert_equal_to (last_local_idx+1, this->size());
447 
448  libmesh_assert_less_equal (send_list.size(), this->size());
449 
450 #ifndef NDEBUG
451  this->_is_closed = true;
452 #endif
453 }
454 
455 
456 
457 template <typename T>
458 void EigenSparseVector<T>::localize (std::vector<T>& v_local) const
459 
460 {
461  v_local.resize(this->size());
462 
463  for (numeric_index_type i=0; i<v_local.size(); i++)
464  v_local[i] = (*this)(i);
465 }
466 
467 
468 
469 template <typename T>
470 void EigenSparseVector<T>::localize_to_one (std::vector<T>& v_local,
471  const processor_id_type libmesh_dbg_var(pid)) const
472 {
473  libmesh_assert_equal_to (pid, 0);
474 
475  this->localize (v_local);
476 }
477 
478 
479 
480 template <typename T>
482  const NumericVector<T>& /*vec2*/)
483 {
484  libmesh_not_implemented();
485 }
486 
487 
488 
489 template <typename T>
491 {
492  libmesh_assert (this->initialized());
493  if (!this->size())
495 
496 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
497  Real the_max = libmesh_real((*this)(0));
498 
499  const numeric_index_type n = this->size();
500 
501  for (numeric_index_type i=1; i<n; i++)
502  the_max = std::max (the_max, libmesh_real((*this)(i)));
503 
504  return the_max;
505 #else
506  return libmesh_real(_vec.maxCoeff());
507 #endif
508 }
509 
510 
511 
512 template <typename T>
514 {
515  libmesh_assert (this->initialized());
516  if (!this->size())
518 
519 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
520  Real the_min = libmesh_real((*this)(0));
521 
522  const numeric_index_type n = this->size();
523 
524  for (numeric_index_type i=1; i<n; i++)
525  the_min = std::min (the_min, libmesh_real((*this)(i)));
526 
527  return the_min;
528 #else
529  return libmesh_real(_vec.minCoeff());
530 #endif
531 }
532 
533 
534 //------------------------------------------------------------------
535 // Explicit instantiations
536 template class EigenSparseVector<Number>;
537 
538 } // namespace libMesh
539 
540 
541 #endif // #ifdef LIBMESH_HAVE_EIGEN

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo