petsc_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 
22 // Local Includes
23 #include "libmesh/petsc_vector.h"
24 #include "libmesh/petsc_matrix.h"
25 
26 #ifdef LIBMESH_HAVE_PETSC
27 
29 #include "libmesh/dense_vector.h"
30 #include "libmesh/parallel.h"
31 #include "libmesh/petsc_macro.h"
32 #include "libmesh/utility.h"
33 
34 namespace libMesh
35 {
36 
37 
38 
39 //-----------------------------------------------------------------------
40 // PetscVector members
41 
42 // void PetscVector<T>::init (const NumericVector<T>& v, const bool fast)
43 // {
44 // libmesh_error();
45 
46 // init (v.local_size(), v.size(), fast);
47 
48 // vec = libmesh_cast_ref<const PetscVector<T>&>(v).vec;
49 // }
50 
51 template <typename T>
53 {
54  this->_restore_array();
55  libmesh_assert(this->closed());
56 
58  PetscScalar value=0.;
59 
60  ierr = VecSum (_vec, &value);
61  LIBMESH_CHKERRABORT(ierr);
62 
63  return static_cast<T>(value);
64 }
65 
66 
67 template <typename T>
69 {
70  this->_restore_array();
71  libmesh_assert(this->closed());
72 
74  PetscReal value=0.;
75 
76  ierr = VecNorm (_vec, NORM_1, &value);
77  LIBMESH_CHKERRABORT(ierr);
78 
79  return static_cast<Real>(value);
80 }
81 
82 
83 
84 template <typename T>
86 {
87  this->_restore_array();
88  libmesh_assert(this->closed());
89 
91  PetscReal value=0.;
92 
93  ierr = VecNorm (_vec, NORM_2, &value);
94  LIBMESH_CHKERRABORT(ierr);
95 
96  return static_cast<Real>(value);
97 }
98 
99 
100 
101 
102 template <typename T>
104 {
105  this->_restore_array();
106  libmesh_assert(this->closed());
107 
109  PetscReal value=0.;
110 
111  ierr = VecNorm (_vec, NORM_INFINITY, &value);
112  LIBMESH_CHKERRABORT(ierr);
113 
114  return static_cast<Real>(value);
115 }
116 
117 
118 
119 
120 template <typename T>
123 {
124  this->_restore_array();
125  libmesh_assert(this->closed());
126 
127  this->add(1., v);
128 
129  return *this;
130 }
131 
132 
133 
134 template <typename T>
137 {
138  this->_restore_array();
139  libmesh_assert(this->closed());
140 
141  this->add(-1., v);
142 
143  return *this;
144 }
145 
146 
147 
148 template <typename T>
149 void PetscVector<T>::set (const numeric_index_type i, const T value)
150 {
151  this->_restore_array();
152  libmesh_assert_less (i, size());
153 
155  PetscInt i_val = static_cast<PetscInt>(i);
156  PetscScalar petsc_value = static_cast<PetscScalar>(value);
157 
158  ierr = VecSetValues (_vec, 1, &i_val, &petsc_value, INSERT_VALUES);
159  LIBMESH_CHKERRABORT(ierr);
160 
161  this->_is_closed = false;
162 }
163 
164 
165 
166 template <typename T>
168 {
169  PetscErrorCode ierr = 0;
170 
171  // VecReciprocal has been in PETSc since at least 2.3.3 days
172  ierr = VecReciprocal(_vec);
173  LIBMESH_CHKERRABORT(ierr);
174 }
175 
176 
177 
178 template <typename T>
180 {
181  PetscErrorCode ierr = 0;
182 
183  // We just call the PETSc VecConjugate
184  ierr = VecConjugate(_vec);
185  LIBMESH_CHKERRABORT(ierr);
186 }
187 
188 
189 
190 template <typename T>
191 void PetscVector<T>::add (const numeric_index_type i, const T value)
192 {
193  this->_restore_array();
194  libmesh_assert_less (i, size());
195 
197  PetscInt i_val = static_cast<PetscInt>(i);
198  PetscScalar petsc_value = static_cast<PetscScalar>(value);
199 
200  ierr = VecSetValues (_vec, 1, &i_val, &petsc_value, ADD_VALUES);
201  LIBMESH_CHKERRABORT(ierr);
202 
203  this->_is_closed = false;
204 }
205 
206 
207 
208 template <typename T>
209 void PetscVector<T>::add_vector (const std::vector<T>& v,
210  const std::vector<numeric_index_type>& dof_indices)
211 {
212  // If we aren't adding anything just return
213  if(v.empty() || dof_indices.empty())
214  return;
215 
216  this->_restore_array();
217  libmesh_assert_equal_to (v.size(), dof_indices.size());
218 
220  const PetscInt * i_val = reinterpret_cast<const PetscInt*>(&dof_indices[0]);
221  const PetscScalar * petsc_value = static_cast<const PetscScalar*>(&v[0]);
222 
223  ierr = VecSetValues (_vec, v.size(), i_val, petsc_value, ADD_VALUES);
224  LIBMESH_CHKERRABORT(ierr);
225 
226  this->_is_closed = false;
227 }
228 
229 
230 
231 template <typename T>
233  const std::vector<numeric_index_type>& dof_indices)
234 {
235  libmesh_assert_equal_to (V.size(), dof_indices.size());
236 
237  for (unsigned int i=0; i<V.size(); i++)
238  this->add (dof_indices[i], V(i));
239 }
240 
241 
242 
243 template <typename T>
245  const SparseMatrix<T>& A_in)
246 {
247  this->_restore_array();
248  // Make sure the data passed in are really of Petsc types
249  const PetscVector<T>* V = libmesh_cast_ptr<const PetscVector<T>*>(&V_in);
250  const PetscMatrix<T>* A = libmesh_cast_ptr<const PetscMatrix<T>*>(&A_in);
251 
253 
254  A->close();
255 
256  // The const_cast<> is not elegant, but it is required since PETSc
257  // is not const-correct.
258  ierr = MatMultAdd(const_cast<PetscMatrix<T>*>(A)->mat(), V->_vec, _vec, _vec);
259  LIBMESH_CHKERRABORT(ierr);
260 }
261 
262 
263 
264 template <typename T>
266  const std::vector<numeric_index_type>& dof_indices)
267 {
268  libmesh_assert_equal_to (V.size(), dof_indices.size());
269 
270  for (unsigned int i=0; i<V.size(); i++)
271  this->add (dof_indices[i], V(i));
272 }
273 
274 
275 template <typename T>
277  const SparseMatrix<T>& A_in)
278 {
279  this->_restore_array();
280  // Make sure the data passed in are really of Petsc types
281  const PetscVector<T>* V = libmesh_cast_ptr<const PetscVector<T>*>(&V_in);
282  const PetscMatrix<T>* A = libmesh_cast_ptr<const PetscMatrix<T>*>(&A_in);
283 
285 
286  A->close();
287 
288  // The const_cast<> is not elegant, but it is required since PETSc
289  // is not const-correct.
290  ierr = MatMultTransposeAdd(const_cast<PetscMatrix<T>*>(A)->mat(), V->_vec, _vec, _vec);
291  LIBMESH_CHKERRABORT(ierr);
292 }
293 
294 
295 #if PETSC_VERSION_LESS_THAN(3,1,0)
296 template <typename T>
298  const SparseMatrix<T>&)
299 {
300 
301  libMesh::out << "MatMultHermitianTranspose was introduced in PETSc 3.1.0,"
302  << "No one has made it backwards compatible with older "
303  << "versions of PETSc so far." << std::endl;
304  libmesh_error();
305 }
306 
307 #else
308 
309 template <typename T>
311  const SparseMatrix<T>& A_in)
312 {
313  this->_restore_array();
314  // Make sure the data passed in are really of Petsc types
315  const PetscVector<T>* V = libmesh_cast_ptr<const PetscVector<T>*>(&V_in);
316  const PetscMatrix<T>* A = libmesh_cast_ptr<const PetscMatrix<T>*>(&A_in);
317 
318  A->close();
319 
320  // Store a temporary copy since MatMultHermitianTransposeAdd doesn't seem to work
321  // TODO: Find out why MatMultHermitianTransposeAdd doesn't work, might be a PETSc bug?
322  AutoPtr< NumericVector<Number> > this_clone = this->clone();
323 
324  // The const_cast<> is not elegant, but it is required since PETSc
325  // is not const-correct.
326  PetscErrorCode ierr = MatMultHermitianTranspose(const_cast<PetscMatrix<T>*>(A)->mat(), V->_vec, _vec);
327  LIBMESH_CHKERRABORT(ierr);
328 
329  // Add the temporary copy to the matvec result
330  this->add(1., *this_clone);
331 }
332 #endif
333 
334 
335 template <typename T>
336 void PetscVector<T>::add (const T v_in)
337 {
338  this->_restore_array();
339 
340  PetscErrorCode ierr=0;
341  PetscScalar* values;
342  const PetscScalar v = static_cast<PetscScalar>(v_in);
343 
344  if(this->type() != GHOSTED)
345  {
346  const PetscInt n = static_cast<PetscInt>(this->local_size());
347  const PetscInt fli = static_cast<PetscInt>(this->first_local_index());
348 
349  for (PetscInt i=0; i<n; i++)
350  {
351  ierr = VecGetArray (_vec, &values);
352  LIBMESH_CHKERRABORT(ierr);
353 
354  PetscInt ig = fli + i;
355 
356  PetscScalar value = (values[i] + v);
357 
358  ierr = VecRestoreArray (_vec, &values);
359  LIBMESH_CHKERRABORT(ierr);
360 
361  ierr = VecSetValues (_vec, 1, &ig, &value, INSERT_VALUES);
362  LIBMESH_CHKERRABORT(ierr);
363  }
364  }
365  else
366  {
367  /* Vectors that include ghost values require a special
368  handling. */
369  Vec loc_vec;
370  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
371  LIBMESH_CHKERRABORT(ierr);
372 
373  PetscInt n=0;
374  ierr = VecGetSize(loc_vec, &n);
375  LIBMESH_CHKERRABORT(ierr);
376 
377  for (PetscInt i=0; i<n; i++)
378  {
379  ierr = VecGetArray (loc_vec, &values);
380  LIBMESH_CHKERRABORT(ierr);
381 
382  PetscScalar value = (values[i] + v);
383 
384  ierr = VecRestoreArray (loc_vec, &values);
385  LIBMESH_CHKERRABORT(ierr);
386 
387  ierr = VecSetValues (loc_vec, 1, &i, &value, INSERT_VALUES);
388  LIBMESH_CHKERRABORT(ierr);
389  }
390 
391  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
392  LIBMESH_CHKERRABORT(ierr);
393  }
394 
395  this->_is_closed = false;
396 }
397 
398 
399 
400 template <typename T>
402 {
403  this->add (1., v);
404 }
405 
406 
407 
408 template <typename T>
409 void PetscVector<T>::add (const T a_in, const NumericVector<T>& v_in)
410 {
411  this->_restore_array();
412 
413  PetscErrorCode ierr = 0;
414  PetscScalar a = static_cast<PetscScalar>(a_in);
415 
416  // Make sure the NumericVector passed in is really a PetscVector
417  const PetscVector<T>* v = libmesh_cast_ptr<const PetscVector<T>*>(&v_in);
418  v->_restore_array();
419 
420  libmesh_assert_equal_to (this->size(), v->size());
421 
422  if(this->type() != GHOSTED)
423  {
424 #if PETSC_VERSION_LESS_THAN(2,3,0)
425  // 2.2.x & earlier style
426  ierr = VecAXPY(&a, v->_vec, _vec);
427  LIBMESH_CHKERRABORT(ierr);
428 #else
429  // 2.3.x & later style
430  ierr = VecAXPY(_vec, a, v->_vec);
431  LIBMESH_CHKERRABORT(ierr);
432 #endif
433  }
434  else
435  {
436  Vec loc_vec;
437  Vec v_loc_vec;
438  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
439  LIBMESH_CHKERRABORT(ierr);
440  ierr = VecGhostGetLocalForm (v->_vec,&v_loc_vec);
441  LIBMESH_CHKERRABORT(ierr);
442 
443 #if PETSC_VERSION_LESS_THAN(2,3,0)
444  // 2.2.x & earlier style
445  ierr = VecAXPY(&a, v_loc_vec, loc_vec);
446  LIBMESH_CHKERRABORT(ierr);
447 #else
448  // 2.3.x & later style
449  ierr = VecAXPY(loc_vec, a, v_loc_vec);
450  LIBMESH_CHKERRABORT(ierr);
451 #endif
452 
453  ierr = VecGhostRestoreLocalForm (v->_vec,&v_loc_vec);
454  LIBMESH_CHKERRABORT(ierr);
455  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
456  LIBMESH_CHKERRABORT(ierr);
457  }
458 }
459 
460 
461 
462 template <typename T>
463 void PetscVector<T>::insert (const std::vector<T>& v,
464  const std::vector<numeric_index_type>& dof_indices)
465 {
466  libmesh_assert_equal_to (v.size(), dof_indices.size());
467 
468  for (unsigned int i=0; i<v.size(); i++)
469  this->set (dof_indices[i], v[i]);
470 }
471 
472 
473 
474 template <typename T>
476  const std::vector<numeric_index_type>& dof_indices)
477 {
478  libmesh_assert_equal_to (V.size(), dof_indices.size());
479 
480  for (unsigned int i=0; i<V.size(); i++)
481  this->set (dof_indices[i], V(i));
482 }
483 
484 
485 
486 template <typename T>
488  const std::vector<numeric_index_type>& dof_indices)
489 {
490  libmesh_assert_equal_to (V.size(), dof_indices.size());
491 
492  for (unsigned int i=0; i<V.size(); i++)
493  this->set (dof_indices[i], V(i));
494 }
495 
496 
497 
498 template <typename T>
500  const std::vector<numeric_index_type>& dof_indices)
501 {
502  libmesh_assert_equal_to (V.size(), dof_indices.size());
503 
504  for (unsigned int i=0; i<V.size(); i++)
505  this->set (dof_indices[i], V(i));
506 }
507 
508 
509 
510 template <typename T>
511 void PetscVector<T>::scale (const T factor_in)
512 {
513  this->_restore_array();
514 
515  PetscErrorCode ierr = 0;
516  PetscScalar factor = static_cast<PetscScalar>(factor_in);
517 
518  if(this->type() != GHOSTED)
519  {
520 #if PETSC_VERSION_LESS_THAN(2,3,0)
521  // 2.2.x & earlier style
522  ierr = VecScale(&factor, _vec);
523  LIBMESH_CHKERRABORT(ierr);
524 #else
525  // 2.3.x & later style
526  ierr = VecScale(_vec, factor);
527  LIBMESH_CHKERRABORT(ierr);
528 #endif
529  }
530  else
531  {
532  Vec loc_vec;
533  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
534  LIBMESH_CHKERRABORT(ierr);
535 
536 #if PETSC_VERSION_LESS_THAN(2,3,0)
537  // 2.2.x & earlier style
538  ierr = VecScale(&factor, loc_vec);
539  LIBMESH_CHKERRABORT(ierr);
540 #else
541  // 2.3.x & later style
542  ierr = VecScale(loc_vec, factor);
543  LIBMESH_CHKERRABORT(ierr);
544 #endif
545 
546  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
547  LIBMESH_CHKERRABORT(ierr);
548  }
549 }
550 
551 template <typename T>
553 {
554  PetscErrorCode ierr = 0;
555 
556  const PetscVector<T>* v_vec = libmesh_cast_ptr<const PetscVector<T>*>(&v);
557 
558  ierr = VecPointwiseDivide(_vec, _vec, v_vec->_vec);
559  LIBMESH_CHKERRABORT(ierr);
560 
561  return *this;
562 }
563 
564 template <typename T>
566 {
567  this->_restore_array();
568 
569  PetscErrorCode ierr = 0;
570 
571  if(this->type() != GHOSTED)
572  {
573  ierr = VecAbs(_vec);
574  LIBMESH_CHKERRABORT(ierr);
575  }
576  else
577  {
578  Vec loc_vec;
579  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
580  LIBMESH_CHKERRABORT(ierr);
581 
582  ierr = VecAbs(loc_vec);
583  LIBMESH_CHKERRABORT(ierr);
584 
585  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
586  LIBMESH_CHKERRABORT(ierr);
587  }
588 }
589 
590 template <typename T>
592 {
593  this->_restore_array();
594 
595  // Error flag
596  PetscErrorCode ierr = 0;
597 
598  // Return value
599  PetscScalar value=0.;
600 
601  // Make sure the NumericVector passed in is really a PetscVector
602  const PetscVector<T>* v = libmesh_cast_ptr<const PetscVector<T>*>(&V);
603 
604  // 2.3.x (at least) style. Untested for previous versions.
605  ierr = VecDot(this->_vec, v->_vec, &value);
606  LIBMESH_CHKERRABORT(ierr);
607 
608  return static_cast<T>(value);
609 }
610 
611 template <typename T>
613 {
614  this->_restore_array();
615 
616  // Error flag
617  PetscErrorCode ierr = 0;
618 
619  // Return value
620  PetscScalar value=0.;
621 
622  // Make sure the NumericVector passed in is really a PetscVector
623  const PetscVector<T>* v = libmesh_cast_ptr<const PetscVector<T>*>(&V);
624 
625  // 2.3.x (at least) style. Untested for previous versions.
626  ierr = VecTDot(this->_vec, v->_vec, &value);
627  LIBMESH_CHKERRABORT(ierr);
628 
629  return static_cast<T>(value);
630 }
631 
632 
633 template <typename T>
636 {
637  this->_restore_array();
638  libmesh_assert(this->closed());
639 
640  PetscErrorCode ierr = 0;
641  PetscScalar s = static_cast<PetscScalar>(s_in);
642 
643  if (this->size() != 0)
644  {
645  if(this->type() != GHOSTED)
646  {
647 #if PETSC_VERSION_LESS_THAN(2,3,0)
648  // 2.2.x & earlier style
649  ierr = VecSet(&s, _vec);
650  LIBMESH_CHKERRABORT(ierr);
651 #else
652  // 2.3.x & later style
653  ierr = VecSet(_vec, s);
654  LIBMESH_CHKERRABORT(ierr);
655 #endif
656  }
657  else
658  {
659  Vec loc_vec;
660  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
661  LIBMESH_CHKERRABORT(ierr);
662 
663 #if PETSC_VERSION_LESS_THAN(2,3,0)
664  // 2.2.x & earlier style
665  ierr = VecSet(&s, loc_vec);
666  LIBMESH_CHKERRABORT(ierr);
667 #else
668  // 2.3.x & later style
669  ierr = VecSet(loc_vec, s);
670  LIBMESH_CHKERRABORT(ierr);
671 #endif
672 
673  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
674  LIBMESH_CHKERRABORT(ierr);
675  }
676  }
677 
678  return *this;
679 }
680 
681 
682 
683 template <typename T>
686 {
687  // Make sure the NumericVector passed in is really a PetscVector
688  const PetscVector<T>* v = libmesh_cast_ptr<const PetscVector<T>*>(&v_in);
689 
690  *this = *v;
691 
692  return *this;
693 }
694 
695 
696 
697 template <typename T>
700 {
701  this->_restore_array();
702  v._restore_array();
703 
704  libmesh_assert_equal_to (this->size(), v.size());
705  libmesh_assert_equal_to (this->local_size(), v.local_size());
706  libmesh_assert (v.closed());
707 
708  PetscErrorCode ierr = 0;
709 
710  if (((this->type()==PARALLEL) && (v.type()==GHOSTED)) ||
711  ((this->type()==GHOSTED) && (v.type()==PARALLEL)) ||
712  ((this->type()==GHOSTED) && (v.type()==SERIAL)) ||
713  ((this->type()==SERIAL) && (v.type()==GHOSTED)))
714  {
715  /* Allow assignment of a ghosted to a parallel vector since this
716  causes no difficulty. See discussion in libmesh-devel of
717  June 24, 2010. */
718  ierr = VecCopy (v._vec, this->_vec);
719  LIBMESH_CHKERRABORT(ierr);
720  }
721  else
722  {
723  /* In all other cases, we assert that both vectors are of equal
724  type. */
725  libmesh_assert_equal_to (this->_type, v._type);
726  libmesh_assert (this->_global_to_local_map == v._global_to_local_map);
727 
728  if (v.size() != 0)
729  {
730  if(this->type() != GHOSTED)
731  {
732  ierr = VecCopy (v._vec, this->_vec);
733  LIBMESH_CHKERRABORT(ierr);
734  }
735  else
736  {
737  Vec loc_vec;
738  Vec v_loc_vec;
739  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
740  LIBMESH_CHKERRABORT(ierr);
741  ierr = VecGhostGetLocalForm (v._vec,&v_loc_vec);
742  LIBMESH_CHKERRABORT(ierr);
743 
744  ierr = VecCopy (v_loc_vec, loc_vec);
745  LIBMESH_CHKERRABORT(ierr);
746 
747  ierr = VecGhostRestoreLocalForm (v._vec,&v_loc_vec);
748  LIBMESH_CHKERRABORT(ierr);
749  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
750  LIBMESH_CHKERRABORT(ierr);
751  }
752  }
753  }
754 
755  close();
756 
757  return *this;
758 }
759 
760 
761 
762 template <typename T>
764 PetscVector<T>::operator = (const std::vector<T>& v)
765 {
766  this->_restore_array();
767 
768  const numeric_index_type nl = this->local_size();
769  const numeric_index_type ioff = this->first_local_index();
770  PetscErrorCode ierr=0;
771  PetscScalar* values;
772 
777  if (this->size() == v.size())
778  {
779  ierr = VecGetArray (_vec, &values);
780  LIBMESH_CHKERRABORT(ierr);
781 
782  for (numeric_index_type i=0; i<nl; i++)
783  values[i] = static_cast<PetscScalar>(v[i+ioff]);
784 
785  ierr = VecRestoreArray (_vec, &values);
786  LIBMESH_CHKERRABORT(ierr);
787  }
788 
793  else
794  {
795  libmesh_assert_equal_to (this->local_size(), v.size());
796 
797  ierr = VecGetArray (_vec, &values);
798  LIBMESH_CHKERRABORT(ierr);
799 
800  for (numeric_index_type i=0; i<nl; i++)
801  values[i] = static_cast<PetscScalar>(v[i]);
802 
803  ierr = VecRestoreArray (_vec, &values);
804  LIBMESH_CHKERRABORT(ierr);
805  }
806 
807  // Make sure ghost dofs are up to date
808  if (this->type() == GHOSTED)
809  this->close();
810 
811  return *this;
812 }
813 
814 
815 
816 template <typename T>
818 {
819  this->_restore_array();
820 
821  // Make sure the NumericVector passed in is really a PetscVector
822  PetscVector<T>* v_local = libmesh_cast_ptr<PetscVector<T>*>(&v_local_in);
823 
824  libmesh_assert(v_local);
825  libmesh_assert_equal_to (v_local->size(), this->size());
826 
827  PetscErrorCode ierr = 0;
828  const PetscInt n = this->size();
829 
830  IS is;
831  VecScatter scatter;
832 
833  // Create idx, idx[i] = i;
834  std::vector<PetscInt> idx(n); Utility::iota (idx.begin(), idx.end(), 0);
835 
836  // Create the index set & scatter object
837  ierr = ISCreateLibMesh(this->comm().get(), n, &idx[0], PETSC_USE_POINTER, &is);
838  LIBMESH_CHKERRABORT(ierr);
839 
840  ierr = VecScatterCreate(_vec, is,
841  v_local->_vec, is,
842  &scatter);
843  LIBMESH_CHKERRABORT(ierr);
844 
845  // Perform the scatter
846 #if PETSC_VERSION_LESS_THAN(2,3,3)
847 
848  ierr = VecScatterBegin(_vec, v_local->_vec, INSERT_VALUES,
849  SCATTER_FORWARD, scatter);
850  LIBMESH_CHKERRABORT(ierr);
851 
852  ierr = VecScatterEnd (_vec, v_local->_vec, INSERT_VALUES,
853  SCATTER_FORWARD, scatter);
854  LIBMESH_CHKERRABORT(ierr);
855 #else
856  // API argument order change in PETSc 2.3.3
857  ierr = VecScatterBegin(scatter, _vec, v_local->_vec,
858  INSERT_VALUES, SCATTER_FORWARD);
859  LIBMESH_CHKERRABORT(ierr);
860 
861  ierr = VecScatterEnd (scatter, _vec, v_local->_vec,
862  INSERT_VALUES, SCATTER_FORWARD);
863  LIBMESH_CHKERRABORT(ierr);
864 #endif
865 
866  // Clean up
867  ierr = LibMeshISDestroy (&is);
868  LIBMESH_CHKERRABORT(ierr);
869 
870  ierr = LibMeshVecScatterDestroy(&scatter);
871  LIBMESH_CHKERRABORT(ierr);
872 
873  // Make sure ghost dofs are up to date
874  if (v_local->type() == GHOSTED)
875  v_local->close();
876 }
877 
878 
879 
880 template <typename T>
882  const std::vector<numeric_index_type>& send_list) const
883 {
884  // FIXME: Workaround for a strange bug at large-scale.
885  // If we have ghosting, PETSc lets us just copy the solution, and
886  // doing so avoids a segfault?
887  if (v_local_in.type() == GHOSTED &&
888  this->type() == PARALLEL)
889  {
890  v_local_in = *this;
891  return;
892  }
893 
894  // Normal code path begins here
895 
896  this->_restore_array();
897 
898  // Make sure the NumericVector passed in is really a PetscVector
899  PetscVector<T>* v_local = libmesh_cast_ptr<PetscVector<T>*>(&v_local_in);
900 
901  libmesh_assert(v_local);
902  libmesh_assert_equal_to (v_local->size(), this->size());
903  libmesh_assert_less_equal (send_list.size(), v_local->size());
904 
905  PetscErrorCode ierr=0;
906  const numeric_index_type n_sl = send_list.size();
907 
908  IS is;
909  VecScatter scatter;
910 
911  std::vector<PetscInt> idx(n_sl + this->local_size());
912 
913  for (numeric_index_type i=0; i<n_sl; i++)
914  idx[i] = static_cast<PetscInt>(send_list[i]);
915  for (numeric_index_type i = 0; i != this->local_size(); ++i)
916  idx[n_sl+i] = i + this->first_local_index();
917 
918  // Create the index set & scatter object
919  if (idx.empty())
920  ierr = ISCreateLibMesh(this->comm().get(),
921  n_sl+this->local_size(), PETSC_NULL, PETSC_USE_POINTER, &is);
922  else
923  ierr = ISCreateLibMesh(this->comm().get(),
924  n_sl+this->local_size(), &idx[0], PETSC_USE_POINTER, &is);
925  LIBMESH_CHKERRABORT(ierr);
926 
927  ierr = VecScatterCreate(_vec, is,
928  v_local->_vec, is,
929  &scatter);
930  LIBMESH_CHKERRABORT(ierr);
931 
932 
933  // Perform the scatter
934 #if PETSC_VERSION_LESS_THAN(2,3,3)
935 
936  ierr = VecScatterBegin(_vec, v_local->_vec, INSERT_VALUES,
937  SCATTER_FORWARD, scatter);
938  LIBMESH_CHKERRABORT(ierr);
939 
940  ierr = VecScatterEnd (_vec, v_local->_vec, INSERT_VALUES,
941  SCATTER_FORWARD, scatter);
942  LIBMESH_CHKERRABORT(ierr);
943 
944 #else
945 
946  // API argument order change in PETSc 2.3.3
947  ierr = VecScatterBegin(scatter, _vec, v_local->_vec,
948  INSERT_VALUES, SCATTER_FORWARD);
949  LIBMESH_CHKERRABORT(ierr);
950 
951  ierr = VecScatterEnd (scatter, _vec, v_local->_vec,
952  INSERT_VALUES, SCATTER_FORWARD);
953  LIBMESH_CHKERRABORT(ierr);
954 
955 #endif
956 
957 
958  // Clean up
959  ierr = LibMeshISDestroy (&is);
960  LIBMESH_CHKERRABORT(ierr);
961 
962  ierr = LibMeshVecScatterDestroy(&scatter);
963  LIBMESH_CHKERRABORT(ierr);
964 
965  // Make sure ghost dofs are up to date
966  if (v_local->type() == GHOSTED)
967  v_local->close();
968 }
969 
970 
971 template <typename T>
972 void PetscVector<T>::localize (const numeric_index_type first_local_idx,
973  const numeric_index_type last_local_idx,
974  const std::vector<numeric_index_type>& send_list)
975 {
976  this->_restore_array();
977 
978  libmesh_assert_less_equal (send_list.size(), this->size());
979  libmesh_assert_less_equal (last_local_idx+1, this->size());
980 
981  const numeric_index_type my_size = this->size();
982  const numeric_index_type my_local_size = (last_local_idx - first_local_idx + 1);
983  PetscErrorCode ierr=0;
984 
985  // Don't bother for serial cases
986 // if ((first_local_idx == 0) &&
987 // (my_local_size == my_size))
988  // But we do need to stay in sync for degenerate cases
989  if (this->n_processors() == 1)
990  return;
991 
992 
993  // Build a parallel vector, initialize it with the local
994  // parts of (*this)
995  PetscVector<T> parallel_vec(this->comm(), PARALLEL);
996 
997  parallel_vec.init (my_size, my_local_size, true, PARALLEL);
998 
999 
1000  // Copy part of *this into the parallel_vec
1001  {
1002  IS is;
1003  VecScatter scatter;
1004 
1005  // Create idx, idx[i] = i+first_local_idx;
1006  std::vector<PetscInt> idx(my_local_size);
1007  Utility::iota (idx.begin(), idx.end(), first_local_idx);
1008 
1009  // Create the index set & scatter object
1010  ierr = ISCreateLibMesh(this->comm().get(), my_local_size,
1011  my_local_size ? &idx[0] : NULL, PETSC_USE_POINTER, &is);
1012  LIBMESH_CHKERRABORT(ierr);
1013 
1014  ierr = VecScatterCreate(_vec, is,
1015  parallel_vec._vec, is,
1016  &scatter);
1017  LIBMESH_CHKERRABORT(ierr);
1018 
1019  // Perform the scatter
1020 #if PETSC_VERSION_LESS_THAN(2,3,3)
1021 
1022  ierr = VecScatterBegin(_vec, parallel_vec._vec, INSERT_VALUES,
1023  SCATTER_FORWARD, scatter);
1024  LIBMESH_CHKERRABORT(ierr);
1025 
1026  ierr = VecScatterEnd (_vec, parallel_vec._vec, INSERT_VALUES,
1027  SCATTER_FORWARD, scatter);
1028  LIBMESH_CHKERRABORT(ierr);
1029 
1030 #else
1031 
1032  // API argument order change in PETSc 2.3.3
1033  ierr = VecScatterBegin(scatter, _vec, parallel_vec._vec,
1034  INSERT_VALUES, SCATTER_FORWARD);
1035  LIBMESH_CHKERRABORT(ierr);
1036 
1037  ierr = VecScatterEnd (scatter, _vec, parallel_vec._vec,
1038  INSERT_VALUES, SCATTER_FORWARD);
1039  LIBMESH_CHKERRABORT(ierr);
1040 
1041 #endif
1042 
1043  // Clean up
1044  ierr = LibMeshISDestroy (&is);
1045  LIBMESH_CHKERRABORT(ierr);
1046 
1047  ierr = LibMeshVecScatterDestroy(&scatter);
1048  LIBMESH_CHKERRABORT(ierr);
1049  }
1050 
1051  // localize like normal
1052  parallel_vec.close();
1053  parallel_vec.localize (*this, send_list);
1054  this->close();
1055 }
1056 
1057 
1058 
1059 template <typename T>
1060 void PetscVector<T>::localize (std::vector<T>& v_local) const
1061 {
1062  this->_restore_array();
1063 
1064  // This function must be run on all processors at once
1065  parallel_object_only();
1066 
1067  PetscErrorCode ierr=0;
1068  const PetscInt n = this->size();
1069  const PetscInt nl = this->local_size();
1070  PetscScalar *values;
1071 
1072  v_local.clear();
1073  v_local.resize(n, 0.);
1074 
1075  ierr = VecGetArray (_vec, &values);
1076  LIBMESH_CHKERRABORT(ierr);
1077 
1078  numeric_index_type ioff = first_local_index();
1079 
1080  for (PetscInt i=0; i<nl; i++)
1081  v_local[i+ioff] = static_cast<T>(values[i]);
1082 
1083  ierr = VecRestoreArray (_vec, &values);
1084  LIBMESH_CHKERRABORT(ierr);
1085 
1086  this->comm().sum(v_local);
1087 }
1088 
1089 
1090 
1091 // Full specialization for Real datatypes
1092 #ifdef LIBMESH_USE_REAL_NUMBERS
1093 
1094 template <>
1095 void PetscVector<Real>::localize_to_one (std::vector<Real>& v_local,
1096  const processor_id_type pid) const
1097 {
1098  this->_restore_array();
1099 
1100  PetscErrorCode ierr=0;
1101  const PetscInt n = size();
1102  const PetscInt nl = local_size();
1103  PetscScalar *values;
1104 
1105 
1106  // only one processor
1107  if (n_processors() == 1)
1108  {
1109  v_local.resize(n);
1110 
1111  ierr = VecGetArray (_vec, &values);
1112  LIBMESH_CHKERRABORT(ierr);
1113 
1114  for (PetscInt i=0; i<n; i++)
1115  v_local[i] = static_cast<Real>(values[i]);
1116 
1117  ierr = VecRestoreArray (_vec, &values);
1118  LIBMESH_CHKERRABORT(ierr);
1119  }
1120 
1121  // otherwise multiple processors
1122  else
1123  {
1124  if(pid == 0) // optimized version for localizing to 0
1125  {
1126  Vec vout;
1127  VecScatter ctx;
1128 
1129  ierr = VecScatterCreateToZero(_vec, &ctx, &vout);
1130  LIBMESH_CHKERRABORT(ierr);
1131 
1132  ierr = VecScatterBegin(ctx, _vec, vout, INSERT_VALUES, SCATTER_FORWARD);
1133  LIBMESH_CHKERRABORT(ierr);
1134  ierr = VecScatterEnd(ctx, _vec, vout, INSERT_VALUES, SCATTER_FORWARD);
1135  LIBMESH_CHKERRABORT(ierr);
1136 
1137  if(processor_id() == 0)
1138  {
1139  v_local.resize(n);
1140 
1141  ierr = VecGetArray (vout, &values);
1142  LIBMESH_CHKERRABORT(ierr);
1143 
1144  for (PetscInt i=0; i<n; i++)
1145  v_local[i] = static_cast<Real>(values[i]);
1146 
1147  ierr = VecRestoreArray (vout, &values);
1148  LIBMESH_CHKERRABORT(ierr);
1149  }
1150 
1151  ierr = LibMeshVecScatterDestroy(&ctx);
1152  LIBMESH_CHKERRABORT(ierr);
1153  ierr = LibMeshVecDestroy(&vout);
1154  LIBMESH_CHKERRABORT(ierr);
1155 
1156  }
1157  else
1158  {
1159  v_local.resize(n);
1160 
1161  numeric_index_type ioff = this->first_local_index();
1162  std::vector<Real> local_values (n, 0.);
1163 
1164  {
1165  ierr = VecGetArray (_vec, &values);
1166  LIBMESH_CHKERRABORT(ierr);
1167 
1168  for (PetscInt i=0; i<nl; i++)
1169  local_values[i+ioff] = static_cast<Real>(values[i]);
1170 
1171  ierr = VecRestoreArray (_vec, &values);
1172  LIBMESH_CHKERRABORT(ierr);
1173  }
1174 
1175 
1176  MPI_Reduce (&local_values[0], &v_local[0], n, MPI_REAL, MPI_SUM,
1177  pid, this->comm().get());
1178  }
1179  }
1180 }
1181 
1182 #endif
1183 
1184 
1185 // Full specialization for Complex datatypes
1186 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
1187 
1188 template <>
1189 void PetscVector<Complex>::localize_to_one (std::vector<Complex>& v_local,
1190  const processor_id_type pid) const
1191 {
1192  this->_restore_array();
1193 
1194  PetscErrorCode ierr=0;
1195  const PetscInt n = size();
1196  const PetscInt nl = local_size();
1197  PetscScalar *values;
1198 
1199 
1200  v_local.resize(n);
1201 
1202 
1203  for (PetscInt i=0; i<n; i++)
1204  v_local[i] = 0.;
1205 
1206  // only one processor
1207  if (n == nl)
1208  {
1209  ierr = VecGetArray (_vec, &values);
1210  LIBMESH_CHKERRABORT(ierr);
1211 
1212  for (PetscInt i=0; i<n; i++)
1213  v_local[i] = static_cast<Complex>(values[i]);
1214 
1215  ierr = VecRestoreArray (_vec, &values);
1216  LIBMESH_CHKERRABORT(ierr);
1217  }
1218 
1219  // otherwise multiple processors
1220  else
1221  {
1222  numeric_index_type ioff = this->first_local_index();
1223 
1224  /* in here the local values are stored, acting as send buffer for MPI
1225  * initialize to zero, since we collect using MPI_SUM
1226  */
1227  std::vector<Real> real_local_values(n, 0.);
1228  std::vector<Real> imag_local_values(n, 0.);
1229 
1230  {
1231  ierr = VecGetArray (_vec, &values);
1232  LIBMESH_CHKERRABORT(ierr);
1233 
1234  // provide my local share to the real and imag buffers
1235  for (PetscInt i=0; i<nl; i++)
1236  {
1237  real_local_values[i+ioff] = static_cast<Complex>(values[i]).real();
1238  imag_local_values[i+ioff] = static_cast<Complex>(values[i]).imag();
1239  }
1240 
1241  ierr = VecRestoreArray (_vec, &values);
1242  LIBMESH_CHKERRABORT(ierr);
1243  }
1244 
1245  /* have buffers of the real and imaginary part of v_local.
1246  * Once MPI_Reduce() collected all the real and imaginary
1247  * parts in these std::vector<Real>, the values can be
1248  * copied to v_local
1249  */
1250  std::vector<Real> real_v_local(n);
1251  std::vector<Real> imag_v_local(n);
1252 
1253  // collect entries from other proc's in real_v_local, imag_v_local
1254  MPI_Reduce (&real_local_values[0], &real_v_local[0], n,
1255  MPI_REAL, MPI_SUM,
1256  pid, this->comm().get());
1257 
1258  MPI_Reduce (&imag_local_values[0], &imag_v_local[0], n,
1259  MPI_REAL, MPI_SUM,
1260  pid, this->comm().get());
1261 
1262  // copy real_v_local and imag_v_local to v_local
1263  for (PetscInt i=0; i<n; i++)
1264  v_local[i] = Complex(real_v_local[i], imag_v_local[i]);
1265  }
1266 }
1267 
1268 #endif
1269 
1270 
1271 
1272 template <typename T>
1274  const NumericVector<T>& vec2)
1275 {
1276  this->_restore_array();
1277 
1278  PetscErrorCode ierr = 0;
1279 
1280  // Convert arguments to PetscVector*.
1281  const PetscVector<T>* vec1_petsc = libmesh_cast_ptr<const PetscVector<T>*>(&vec1);
1282  const PetscVector<T>* vec2_petsc = libmesh_cast_ptr<const PetscVector<T>*>(&vec2);
1283 
1284  // Call PETSc function.
1285 
1286 #if PETSC_VERSION_LESS_THAN(2,3,1)
1287 
1288  libMesh::out << "This method has been developed with PETSc 2.3.1. "
1289  << "No one has made it backwards compatible with older "
1290  << "versions of PETSc so far; however, it might work "
1291  << "without any change with some older version." << std::endl;
1292  libmesh_error();
1293 
1294 #else
1295 
1296  if(this->type() != GHOSTED)
1297  {
1298  ierr = VecPointwiseMult(this->vec(),
1299  const_cast<PetscVector<T>*>(vec1_petsc)->vec(),
1300  const_cast<PetscVector<T>*>(vec2_petsc)->vec());
1301  LIBMESH_CHKERRABORT(ierr);
1302  }
1303  else
1304  {
1305  Vec loc_vec;
1306  Vec v1_loc_vec;
1307  Vec v2_loc_vec;
1308  ierr = VecGhostGetLocalForm (_vec,&loc_vec);
1309  LIBMESH_CHKERRABORT(ierr);
1310  ierr = VecGhostGetLocalForm (const_cast<PetscVector<T>*>(vec1_petsc)->vec(),&v1_loc_vec);
1311  LIBMESH_CHKERRABORT(ierr);
1312  ierr = VecGhostGetLocalForm (const_cast<PetscVector<T>*>(vec2_petsc)->vec(),&v2_loc_vec);
1313  LIBMESH_CHKERRABORT(ierr);
1314 
1315  ierr = VecPointwiseMult(loc_vec,v1_loc_vec,v2_loc_vec);
1316  LIBMESH_CHKERRABORT(ierr);
1317 
1318  ierr = VecGhostRestoreLocalForm (const_cast<PetscVector<T>*>(vec1_petsc)->vec(),&v1_loc_vec);
1319  LIBMESH_CHKERRABORT(ierr);
1320  ierr = VecGhostRestoreLocalForm (const_cast<PetscVector<T>*>(vec2_petsc)->vec(),&v2_loc_vec);
1321  LIBMESH_CHKERRABORT(ierr);
1322  ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
1323  LIBMESH_CHKERRABORT(ierr);
1324  }
1325 
1326 #endif
1327 
1328 }
1329 
1330 
1331 
1332 template <typename T>
1333 void PetscVector<T>::print_matlab (const std::string name) const
1334 {
1335  this->_restore_array();
1336  libmesh_assert (this->closed());
1337 
1338  PetscErrorCode ierr=0;
1339  PetscViewer petsc_viewer;
1340 
1341 
1342  ierr = PetscViewerCreate (this->comm().get(),
1343  &petsc_viewer);
1344  LIBMESH_CHKERRABORT(ierr);
1345 
1350  if (name != "NULL")
1351  {
1352  ierr = PetscViewerASCIIOpen( this->comm().get(),
1353  name.c_str(),
1354  &petsc_viewer);
1355  LIBMESH_CHKERRABORT(ierr);
1356 
1357  ierr = PetscViewerSetFormat (petsc_viewer,
1358  PETSC_VIEWER_ASCII_MATLAB);
1359  LIBMESH_CHKERRABORT(ierr);
1360 
1361  ierr = VecView (_vec, petsc_viewer);
1362  LIBMESH_CHKERRABORT(ierr);
1363  }
1364 
1368  else
1369  {
1370  ierr = PetscViewerSetFormat (PETSC_VIEWER_STDOUT_WORLD,
1371  PETSC_VIEWER_ASCII_MATLAB);
1372  LIBMESH_CHKERRABORT(ierr);
1373 
1374  ierr = VecView (_vec, PETSC_VIEWER_STDOUT_WORLD);
1375  LIBMESH_CHKERRABORT(ierr);
1376  }
1377 
1378 
1382  ierr = LibMeshPetscViewerDestroy (&petsc_viewer);
1383  LIBMESH_CHKERRABORT(ierr);
1384 }
1385 
1386 
1387 
1388 
1389 
1390 template <typename T>
1392  const std::vector<numeric_index_type>& rows) const
1393 {
1394  this->_restore_array();
1395 
1396  // PETSc data structures
1397  IS parent_is, subvector_is;
1398  VecScatter scatter;
1399  PetscErrorCode ierr = 0;
1400 
1401  // Make sure the passed in subvector is really a PetscVector
1402  PetscVector<T>* petsc_subvector = libmesh_cast_ptr<PetscVector<T>*>(&subvector);
1403 
1404  // If the petsc_subvector is already initialized, we assume that the
1405  // user has already allocated the *correct* amount of space for it.
1406  // If not, we use the appropriate PETSc routines to initialize it.
1407  if (!petsc_subvector->initialized())
1408  {
1409  // Initialize the petsc_subvector to have enough space to hold
1410  // the entries which will be scattered into it. Note: such an
1411  // init() function (where we let PETSc decide the number of local
1412  // entries) is not currently offered by the PetscVector
1413  // class. Should we differentiate here between sequential and
1414  // parallel vector creation based on this->n_processors() ?
1415  ierr = VecCreateMPI(this->comm().get(),
1416  PETSC_DECIDE, // n_local
1417  rows.size(), // n_global
1418  &(petsc_subvector->_vec)); LIBMESH_CHKERRABORT(ierr);
1419 
1420  ierr = VecSetFromOptions (petsc_subvector->_vec); LIBMESH_CHKERRABORT(ierr);
1421 
1422  // Mark the subvector as initialized
1423  petsc_subvector->_is_initialized = true;
1424  }
1425  else
1426  {
1427  petsc_subvector->_restore_array();
1428  }
1429 
1430  // Use iota to fill an array with entries [0,1,2,3,4,...rows.size()]
1431  std::vector<PetscInt> idx(rows.size());
1432  Utility::iota (idx.begin(), idx.end(), 0);
1433 
1434  // Construct index sets
1435  ierr = ISCreateLibMesh(this->comm().get(),
1436  rows.size(),
1437  (PetscInt*) &rows[0],
1439  &parent_is); LIBMESH_CHKERRABORT(ierr);
1440 
1441  ierr = ISCreateLibMesh(this->comm().get(),
1442  rows.size(),
1443  (PetscInt*) &idx[0],
1445  &subvector_is); LIBMESH_CHKERRABORT(ierr);
1446 
1447  // Construct the scatter object
1448  ierr = VecScatterCreate(this->_vec,
1449  parent_is,
1450  petsc_subvector->_vec,
1451  subvector_is,
1452  &scatter); LIBMESH_CHKERRABORT(ierr);
1453 
1454  // Actually perform the scatter
1455 #if PETSC_VERSION_LESS_THAN(2,3,3)
1456  ierr = VecScatterBegin(this->_vec,
1457  petsc_subvector->_vec,
1458  INSERT_VALUES,
1459  SCATTER_FORWARD,
1460  scatter); LIBMESH_CHKERRABORT(ierr);
1461 
1462  ierr = VecScatterEnd(this->_vec,
1463  petsc_subvector->_vec,
1464  INSERT_VALUES,
1465  SCATTER_FORWARD,
1466  scatter); LIBMESH_CHKERRABORT(ierr);
1467 #else
1468  // API argument order change in PETSc 2.3.3
1469  ierr = VecScatterBegin(scatter,
1470  this->_vec,
1471  petsc_subvector->_vec,
1472  INSERT_VALUES,
1473  SCATTER_FORWARD); LIBMESH_CHKERRABORT(ierr);
1474 
1475  ierr = VecScatterEnd(scatter,
1476  this->_vec,
1477  petsc_subvector->_vec,
1478  INSERT_VALUES,
1479  SCATTER_FORWARD); LIBMESH_CHKERRABORT(ierr);
1480 #endif
1481 
1482  // Clean up
1483  ierr = LibMeshISDestroy(&parent_is); LIBMESH_CHKERRABORT(ierr);
1484  ierr = LibMeshISDestroy(&subvector_is); LIBMESH_CHKERRABORT(ierr);
1485  ierr = LibMeshVecScatterDestroy(&scatter); LIBMESH_CHKERRABORT(ierr);
1486 
1487 }
1488 
1489 
1490 
1491 
1492 //------------------------------------------------------------------
1493 // Explicit instantiations
1494 template class PetscVector<Number>;
1495 
1496 } // namespace libMesh
1497 
1498 
1499 
1500 #endif // #ifdef LIBMESH_HAVE_PETSC

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

Hosted By:
SourceForge.net Logo