distributed_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 #include "libmesh/libmesh_common.h"
21 
22 // C++ includes
23 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
24 #include <cmath> // for std::abs
25 #include <limits> // std::numeric_limits<T>::min()
26 
27 // Local Includes
29 #include "libmesh/dense_vector.h"
31 #include "libmesh/parallel.h"
32 #include "libmesh/tensor_tools.h"
33 
34 namespace libMesh
35 {
36 
37 
38 
39 //--------------------------------------------------------------------------
40 // DistributedVector methods
41 template <typename T>
43 {
44  // This function must be run on all processors at once
45  parallel_object_only();
46 
47  libmesh_assert (this->initialized());
48  libmesh_assert_equal_to (_values.size(), _local_size);
49  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
50 
51  T local_sum = 0.;
52 
53  for (numeric_index_type i=0; i<local_size(); i++)
54  local_sum += _values[i];
55 
56  this->comm().sum(local_sum);
57 
58  return local_sum;
59 }
60 
61 
62 
63 template <typename T>
65 {
66  // This function must be run on all processors at once
67  parallel_object_only();
68 
69  libmesh_assert (this->initialized());
70  libmesh_assert_equal_to (_values.size(), _local_size);
71  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
72 
73  double local_l1 = 0.;
74 
75  for (numeric_index_type i=0; i<local_size(); i++)
76  local_l1 += std::abs(_values[i]);
77 
78  this->comm().sum(local_l1);
79 
80  return local_l1;
81 }
82 
83 
84 
85 template <typename T>
87 {
88  // This function must be run on all processors at once
89  parallel_object_only();
90 
91  libmesh_assert (this->initialized());
92  libmesh_assert_equal_to (_values.size(), _local_size);
93  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
94 
95  double local_l2 = 0.;
96 
97  for (numeric_index_type i=0; i<local_size(); i++)
98  local_l2 += TensorTools::norm_sq(_values[i]);
99 
100  this->comm().sum(local_l2);
101 
102  return std::sqrt(local_l2);
103 }
104 
105 
106 
107 template <typename T>
109 {
110  // This function must be run on all processors at once
111  parallel_object_only();
112 
113  libmesh_assert (this->initialized());
114  libmesh_assert_equal_to (_values.size(), _local_size);
115  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
116 
117  Real local_linfty = 0.;
118 
119  for (numeric_index_type i=0; i<local_size(); i++)
120  local_linfty = std::max(local_linfty,
121  static_cast<Real>(std::abs(_values[i]))
122  ); // Note we static_cast so that both
123  // types are the same, as required
124  // by std::max
125 
126  this->comm().max(local_linfty);
127 
128  return local_linfty;
129 }
130 
131 
132 
133 template <typename T>
135 {
136  libmesh_assert (this->closed());
137  libmesh_assert (this->initialized());
138  libmesh_assert_equal_to (_values.size(), _local_size);
139  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
140 
141  add(1., v);
142 
143  return *this;
144 }
145 
146 
147 
148 template <typename T>
150 {
151  libmesh_assert (this->closed());
152  libmesh_assert (this->initialized());
153  libmesh_assert_equal_to (_values.size(), _local_size);
154  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
155 
156  add(-1., v);
157 
158  return *this;
159 }
160 
161 
162 
163 template <typename T>
165 {
166  libmesh_assert_equal_to(size(), v.size());
167 
168  DistributedVector<T> & v_vec = libmesh_cast_ref<DistributedVector<T>&>(v);
169 
170  std::size_t size = _values.size();
171 
172  for(std::size_t i=0; i<size; i++)
173  _values[i] = _values[i] / v_vec._values[i];
174 
175  return *this;
176 }
177 
178 
179 
180 
181 template <typename T>
182 void DistributedVector<T>::add_vector (const std::vector<T>& v,
183  const std::vector<numeric_index_type>& dof_indices)
184 {
185  libmesh_assert (!v.empty());
186  libmesh_assert_equal_to (v.size(), dof_indices.size());
187  libmesh_assert (this->initialized());
188  libmesh_assert_equal_to (_values.size(), _local_size);
189  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
190 
191  for (std::size_t i=0; i<v.size(); i++)
192  add (dof_indices[i], v[i]);
193 }
194 
195 
196 
197 template <typename T>
199  const std::vector<numeric_index_type>& dof_indices)
200 {
201  libmesh_assert_equal_to (V.size(), dof_indices.size());
202  libmesh_assert (this->initialized());
203  libmesh_assert_equal_to (_values.size(), _local_size);
204  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
205 
206  for (std::size_t i=0; i<V.size(); i++)
207  add (dof_indices[i], V(i));
208 }
209 
210 
211 
212 template <typename T>
214  const std::vector<numeric_index_type>& dof_indices)
215 {
216  libmesh_assert_equal_to (V.size(), dof_indices.size());
217  libmesh_assert (this->initialized());
218  libmesh_assert_equal_to (_values.size(), _local_size);
219  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
220 
221  for (unsigned int i=0; i<V.size(); i++)
222  add (dof_indices[i], V(i));
223 }
224 
225 
226 
227 
228 template <typename T>
230 {
231  for (numeric_index_type i=0; i<local_size(); i++)
232  {
233  // Don't divide by zero
234  libmesh_assert_not_equal_to (_values[i], T(0));
235 
236  _values[i] = 1. / _values[i];
237  }
238 }
239 
240 
241 
242 
243 template <typename T>
245 {
246  for (numeric_index_type i=0; i<local_size(); i++)
247  {
248  // Replace values by complex conjugate
249  _values[i] = libmesh_conj(_values[i]);
250  }
251 }
252 
253 
254 
255 
256 
257 template <typename T>
259 {
260  libmesh_assert (this->initialized());
261  libmesh_assert_equal_to (_values.size(), _local_size);
262  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
263 
264  for (numeric_index_type i=0; i<local_size(); i++)
265  _values[i] += v;
266 }
267 
268 
269 
270 template <typename T>
272 {
273  libmesh_assert (this->initialized());
274  libmesh_assert_equal_to (_values.size(), _local_size);
275  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
276 
277  add (1., v);
278 }
279 
280 
281 
282 template <typename T>
283 void DistributedVector<T>::add (const T a, const NumericVector<T>& v)
284 {
285  libmesh_assert (this->initialized());
286  libmesh_assert_equal_to (_values.size(), _local_size);
287  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
288 
289  add(a, v);
290 }
291 
292 
293 
294 template <typename T>
295 void DistributedVector<T>::insert (const std::vector<T>& v,
296  const std::vector<numeric_index_type>& dof_indices)
297 {
298  libmesh_assert (!v.empty());
299  libmesh_assert_equal_to (v.size(), dof_indices.size());
300  libmesh_assert (this->initialized());
301  libmesh_assert_equal_to (_values.size(), _local_size);
302  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
303 
304  for (std::size_t i=0; i<v.size(); i++)
305  this->set (dof_indices[i], v[i]);
306 }
307 
308 
309 
310 template <typename T>
312  const std::vector<numeric_index_type>& dof_indices)
313 {
314  libmesh_assert_equal_to (V.size(), dof_indices.size());
315  libmesh_assert (this->initialized());
316  libmesh_assert_equal_to (_values.size(), _local_size);
317  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
318 
319  for (std::size_t i=0; i<V.size(); i++)
320  this->set (dof_indices[i], V(i));
321 }
322 
323 
324 
325 template <typename T>
327  const std::vector<numeric_index_type>& dof_indices)
328 {
329  libmesh_assert_equal_to (V.size(), dof_indices.size());
330  libmesh_assert (this->initialized());
331  libmesh_assert_equal_to (_values.size(), _local_size);
332  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
333 
334  for (unsigned int i=0; i<V.size(); i++)
335  this->set (dof_indices[i], V(i));
336 }
337 
338 
339 
340 template <typename T>
342  const std::vector<numeric_index_type>& dof_indices)
343 {
344  libmesh_assert_equal_to (V.size(), dof_indices.size());
345  libmesh_assert (this->initialized());
346  libmesh_assert_equal_to (_values.size(), _local_size);
347  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
348 
349  for (unsigned int i=0; i<V.size(); i++)
350  this->set (dof_indices[i], V(i));
351 }
352 
353 
354 
355 template <typename T>
356 void DistributedVector<T>::scale (const T factor)
357 {
358  libmesh_assert (this->initialized());
359  libmesh_assert_equal_to (_values.size(), _local_size);
360  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
361 
362  for (std::size_t i=0; i<local_size(); i++)
363  _values[i] *= factor;
364 }
365 
366 template <typename T>
368 {
369  libmesh_assert (this->initialized());
370  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
371 
372  for (std::size_t i=0; i<local_size(); i++)
373  this->set(i,std::abs(_values[i]));
374 }
375 
376 
377 
378 
379 
380 template <typename T>
382 {
383  // This function must be run on all processors at once
384  parallel_object_only();
385 
386  // Make sure the NumericVector passed in is really a DistributedVector
387  const DistributedVector<T>* v = libmesh_cast_ptr<const DistributedVector<T>*>(&V);
388 
389  // Make sure that the two vectors are distributed in the same way.
390  libmesh_assert_equal_to ( this->first_local_index(), v->first_local_index() );
391  libmesh_assert_equal_to ( this->last_local_index(), v->last_local_index() );
392 
393  // The result of dotting together the local parts of the vector.
394  T local_dot = 0;
395 
396  for (std::size_t i=0; i<this->local_size(); i++)
397  local_dot += this->_values[i] * v->_values[i];
398 
399  // The local dot products are now summed via MPI
400  this->comm().sum(local_dot);
401 
402  return local_dot;
403 }
404 
405 
406 
407 template <typename T>
410 {
411  libmesh_assert (this->initialized());
412  libmesh_assert_equal_to (_values.size(), _local_size);
413  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
414 
415  for (std::size_t i=0; i<local_size(); i++)
416  _values[i] = s;
417 
418  return *this;
419 }
420 
421 
422 
423 template <typename T>
426 {
427  // Make sure the NumericVector passed in is really a DistributedVector
428  const DistributedVector<T>* v = libmesh_cast_ptr<const DistributedVector<T>*>(&v_in);
429 
430  *this = *v;
431 
432  return *this;
433 }
434 
435 
436 
437 template <typename T>
440 {
442  this->_is_closed = v._is_closed;
443 
444  _global_size = v._global_size;
445  _local_size = v._local_size;
446  _first_local_index = v._first_local_index;
447  _last_local_index = v._last_local_index;
448 
449  if (v.local_size() == this->local_size())
450  {
451  _values = v._values;
452  }
453  else
454  {
455  libmesh_error();
456  }
457 
458  return *this;
459 }
460 
461 
462 
463 template <typename T>
465 DistributedVector<T>::operator = (const std::vector<T>& v)
466 {
467  libmesh_assert (this->initialized());
468  libmesh_assert_equal_to (_values.size(), _local_size);
469  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
470 
471  if (v.size() == local_size())
472  _values = v;
473 
474  else if (v.size() == size())
475  for (std::size_t i=first_local_index(); i<last_local_index(); i++)
476  _values[i-first_local_index()] = v[i];
477 
478  else
479  {
480  libmesh_error();
481  }
482 
483 
484  return *this;
485 }
486 
487 
488 
489 template <typename T>
491 
492 {
493  libmesh_assert (this->initialized());
494  libmesh_assert_equal_to (_values.size(), _local_size);
495  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
496 
497  DistributedVector<T>* v_local = libmesh_cast_ptr<DistributedVector<T>*>(&v_local_in);
498 
499  v_local->_first_local_index = 0;
500 
501  v_local->_global_size =
502  v_local->_local_size =
503  v_local->_last_local_index = size();
504 
505  v_local->_is_initialized =
506  v_local->_is_closed = true;
507 
508  // Call localize on the vector's values. This will help
509  // prevent code duplication
510  localize (v_local->_values);
511 
512 #ifndef LIBMESH_HAVE_MPI
513 
514  libmesh_assert_equal_to (local_size(), size());
515 
516 #endif
517 }
518 
519 
520 
521 template <typename T>
523  const std::vector<numeric_index_type>&) const
524 {
525  libmesh_assert (this->initialized());
526  libmesh_assert_equal_to (_values.size(), _local_size);
527  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
528 
529  // TODO: We don't yet support the send list; this is inefficient:
530  localize (v_local_in);
531 }
532 
533 
534 
535 template <typename T>
537  const numeric_index_type last_local_idx,
538  const std::vector<numeric_index_type>& send_list)
539 {
540  // Only good for serial vectors
541  libmesh_assert_equal_to (this->size(), this->local_size());
542  libmesh_assert_greater (last_local_idx, first_local_idx);
543  libmesh_assert_less_equal (send_list.size(), this->size());
544  libmesh_assert_less (last_local_idx, this->size());
545 
546  const numeric_index_type my_size = this->size();
547  const numeric_index_type my_local_size = (last_local_idx - first_local_idx + 1);
548 
549  // Don't bother for serial cases
550  if ((first_local_idx == 0) &&
551  (my_local_size == my_size))
552  return;
553 
554 
555  // Build a parallel vector, initialize it with the local
556  // parts of (*this)
557  DistributedVector<T> parallel_vec(this->comm());
558 
559  parallel_vec.init (my_size, my_local_size, true, PARALLEL);
560 
561  // Copy part of *this into the parallel_vec
562  for (numeric_index_type i=first_local_idx; i<=last_local_idx; i++)
563  parallel_vec._values[i-first_local_idx] = _values[i];
564 
565  // localize like normal
566  parallel_vec.localize (*this, send_list);
567 }
568 
569 
570 
571 template <typename T>
572 void DistributedVector<T>::localize (std::vector<T>& v_local) const
573 {
574  // This function must be run on all processors at once
575  parallel_object_only();
576 
577  libmesh_assert (this->initialized());
578  libmesh_assert_equal_to (_values.size(), _local_size);
579  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
580 
581  v_local = this->_values;
582 
583  this->comm().allgather (v_local);
584 
585 #ifndef LIBMESH_HAVE_MPI
586  libmesh_assert_equal_to (local_size(), size());
587 #endif
588 }
589 
590 
591 
592 template <typename T>
593 void DistributedVector<T>::localize_to_one (std::vector<T>& v_local,
594  const processor_id_type pid) const
595 {
596  // This function must be run on all processors at once
597  parallel_object_only();
598 
599  libmesh_assert (this->initialized());
600  libmesh_assert_equal_to (_values.size(), _local_size);
601  libmesh_assert_equal_to ((_last_local_index - _first_local_index), _local_size);
602 
603  v_local = this->_values;
604 
605  this->comm().gather (pid, v_local);
606 
607 #ifndef LIBMESH_HAVE_MPI
608  libmesh_assert_equal_to (local_size(), size());
609 #endif
610 }
611 
612 
613 
614 template <typename T>
616  const NumericVector<T>&)
617 //void DistributedVector<T>::pointwise_mult (const NumericVector<T>& vec1,
618 // const NumericVector<T>& vec2)
619 {
620  libmesh_not_implemented();
621 }
622 
623 
624 
625 //--------------------------------------------------------------
626 // Explicit instantiations
627 template class DistributedVector<Number>;
628 
629 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo