numeric_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 <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
22 #include <cmath> // for std::abs
23 #include <limits>
24 
25 // Local Includes
26 #include "libmesh/numeric_vector.h"
28 #include "libmesh/laspack_vector.h"
30 #include "libmesh/petsc_vector.h"
32 #include "libmesh/shell_matrix.h"
33 #include "libmesh/tensor_tools.h"
34 
35 namespace libMesh
36 {
37 
38 
39 
40 //------------------------------------------------------------------
41 // NumericVector methods
42 
43 // Full specialization for Real datatypes
44 template <typename T>
45 AutoPtr<NumericVector<T> >
47 {
48  // Build the appropriate vector
49  switch (solver_package)
50  {
51 
52 
53 #ifdef LIBMESH_HAVE_LASPACK
54  case LASPACK_SOLVERS:
55  {
57  return ap;
58  }
59 #endif
60 
61 
62 #ifdef LIBMESH_HAVE_PETSC
63  case PETSC_SOLVERS:
64  {
66  return ap;
67  }
68 #endif
69 
70 
71 #ifdef LIBMESH_HAVE_TRILINOS
72  case TRILINOS_SOLVERS:
73  {
75  return ap;
76  }
77 #endif
78 
79 
80 #ifdef LIBMESH_HAVE_EIGEN
81  case EIGEN_SOLVERS:
82  {
84  return ap;
85  }
86 #endif
87 
88 
89  default:
91  return ap;
92 
93  }
94 
95  AutoPtr<NumericVector<T> > ap(NULL);
96  return ap;
97 }
98 
99 
100 #ifndef LIBMESH_DISABLE_COMMWORLD
101 template <typename T>
104 {
105  libmesh_deprecated();
106  return NumericVector<T>::build(CommWorld, solver_package);
107 }
108 #endif
109 
110 
111 template <typename T>
113  const Real threshold) const
114 {
115  libmesh_assert (this->initialized());
116  libmesh_assert (other_vector.initialized());
117  libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
118  libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
119 
120  int first_different_i = std::numeric_limits<int>::max();
121  numeric_index_type i = first_local_index();
122 
123  do
124  {
125  if ( std::abs( (*this)(i) - other_vector(i) ) > threshold )
126  first_different_i = i;
127  else
128  i++;
129  }
130  while (first_different_i==std::numeric_limits<int>::max()
131  && i<last_local_index());
132 
133  // Find the correct first differing index in parallel
134  this->comm().min(first_different_i);
135 
136  if (first_different_i == std::numeric_limits<int>::max())
137  return -1;
138 
139  return first_different_i;
140 }
141 
142 
143 template <typename T>
145  const Real threshold) const
146 {
147  libmesh_assert (this->initialized());
148  libmesh_assert (other_vector.initialized());
149  libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
150  libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
151 
152  int first_different_i = std::numeric_limits<int>::max();
153  numeric_index_type i = first_local_index();
154 
155  do
156  {
157  if ( std::abs( (*this)(i) - other_vector(i) ) > threshold *
158  std::max(std::abs((*this)(i)), std::abs(other_vector(i))))
159  first_different_i = i;
160  else
161  i++;
162  }
163  while (first_different_i==std::numeric_limits<int>::max()
164  && i<last_local_index());
165 
166  // Find the correct first differing index in parallel
167  this->comm().min(first_different_i);
168 
169  if (first_different_i == std::numeric_limits<int>::max())
170  return -1;
171 
172  return first_different_i;
173 }
174 
175 
176 template <typename T>
178  const Real threshold) const
179 {
180  libmesh_assert (this->initialized());
181  libmesh_assert (other_vector.initialized());
182  libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
183  libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
184 
185  int first_different_i = std::numeric_limits<int>::max();
186  numeric_index_type i = first_local_index();
187 
188  const Real my_norm = this->linfty_norm();
189  const Real other_norm = other_vector.linfty_norm();
190  const Real abs_threshold = std::max(my_norm, other_norm) * threshold;
191 
192  do
193  {
194  if ( std::abs( (*this)(i) - other_vector(i) ) > abs_threshold )
195  first_different_i = i;
196  else
197  i++;
198  }
199  while (first_different_i==std::numeric_limits<int>::max()
200  && i<last_local_index());
201 
202  // Find the correct first differing index in parallel
203  this->comm().min(first_different_i);
204 
205  if (first_different_i == std::numeric_limits<int>::max())
206  return -1;
207 
208  return first_different_i;
209 }
210 
211 /*
212 // Full specialization for float datatypes (DistributedVector wants this)
213 
214 template <>
215 int NumericVector<float>::compare (const NumericVector<float> &other_vector,
216  const Real threshold) const
217 {
218  libmesh_assert (this->initialized());
219  libmesh_assert (other_vector.initialized());
220  libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
221  libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
222 
223  int rvalue = -1;
224  numeric_index_type i = first_local_index();
225 
226  do
227  {
228  if ( std::abs( (*this)(i) - other_vector(i) ) > threshold )
229  rvalue = i;
230  else
231  i++;
232  }
233  while (rvalue==-1 && i<last_local_index());
234 
235  return rvalue;
236 }
237 
238 // Full specialization for double datatypes
239 template <>
240 int NumericVector<double>::compare (const NumericVector<double> &other_vector,
241  const Real threshold) const
242 {
243  libmesh_assert (this->initialized());
244  libmesh_assert (other_vector.initialized());
245  libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
246  libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
247 
248  int rvalue = -1;
249  numeric_index_type i = first_local_index();
250 
251  do
252  {
253  if ( std::abs( (*this)(i) - other_vector(i) ) > threshold )
254  rvalue = i;
255  else
256  i++;
257  }
258  while (rvalue==-1 && i<last_local_index());
259 
260  return rvalue;
261 }
262 
263 #ifdef LIBMESH_DEFAULT_TRIPLE_PRECISION
264 // Full specialization for long double datatypes
265 template <>
266 int NumericVector<long double>::compare (const NumericVector<long double> &other_vector,
267  const Real threshold) const
268 {
269  libmesh_assert (this->initialized());
270  libmesh_assert (other_vector.initialized());
271  libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
272  libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
273 
274  int rvalue = -1;
275  numeric_index_type i = first_local_index();
276 
277  do
278  {
279  if ( std::abs( (*this)(i) - other_vector(i) ) > threshold )
280  rvalue = i;
281  else
282  i++;
283  }
284  while (rvalue==-1 && i<last_local_index());
285 
286  return rvalue;
287 }
288 #endif
289 
290 
291 // Full specialization for Complex datatypes
292 template <>
293 int NumericVector<Complex>::compare (const NumericVector<Complex> &other_vector,
294  const Real threshold) const
295 {
296  libmesh_assert (this->initialized());
297  libmesh_assert (other_vector.initialized());
298  libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
299  libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
300 
301  int rvalue = -1;
302  numeric_index_type i = first_local_index();
303 
304  do
305  {
306  if (( std::abs( (*this)(i).real() - other_vector(i).real() ) > threshold ) ||
307  ( std::abs( (*this)(i).imag() - other_vector(i).imag() ) > threshold ))
308  rvalue = i;
309  else
310  i++;
311  }
312  while (rvalue==-1 && i<this->last_local_index());
313 
314  return rvalue;
315 }
316 */
317 
318 
319 template <class T>
320 Real NumericVector<T>::subset_l1_norm (const std::set<numeric_index_type> & indices) const
321 {
322  const NumericVector<T> & v = *this;
323 
324  std::set<numeric_index_type>::const_iterator it = indices.begin();
325  const std::set<numeric_index_type>::const_iterator it_end = indices.end();
326 
327  Real norm = 0;
328 
329  for(; it!=it_end; ++it)
330  norm += std::abs(v(*it));
331 
332  this->comm().sum(norm);
333 
334  return norm;
335 }
336 
337 template <class T>
338 Real NumericVector<T>::subset_l2_norm (const std::set<numeric_index_type> & indices) const
339 {
340  const NumericVector<T> & v = *this;
341 
342  std::set<numeric_index_type>::const_iterator it = indices.begin();
343  const std::set<numeric_index_type>::const_iterator it_end = indices.end();
344 
345  Real norm = 0;
346 
347  for(; it!=it_end; ++it)
348  norm += TensorTools::norm_sq(v(*it));
349 
350  this->comm().sum(norm);
351 
352  return std::sqrt(norm);
353 }
354 
355 template <class T>
356 Real NumericVector<T>::subset_linfty_norm (const std::set<numeric_index_type> & indices) const
357 {
358  const NumericVector<T> & v = *this;
359 
360  std::set<numeric_index_type>::const_iterator it = indices.begin();
361  const std::set<numeric_index_type>::const_iterator it_end = indices.end();
362 
363  Real norm = 0;
364 
365  for(; it!=it_end; ++it)
366  {
367  Real value = std::abs(v(*it));
368  if(value > norm)
369  norm = value;
370  }
371 
372  this->comm().max(norm);
373 
374  return norm;
375 }
376 
377 
378 
379 template <typename T>
381  const ShellMatrix<T>& a)
382 {
383  a.vector_mult_add(*this,v);
384 }
385 
386 
387 
388 //------------------------------------------------------------------
389 // Explicit instantiations
390 template class NumericVector<Number>;
391 
392 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo