Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 : // Local Includes
20 : #include "libmesh/numeric_vector.h"
21 : #include "libmesh/distributed_vector.h"
22 : #include "libmesh/laspack_vector.h"
23 : #include "libmesh/eigen_sparse_vector.h"
24 : #include "libmesh/petsc_vector.h"
25 : #include "libmesh/trilinos_epetra_vector.h"
26 : #include "libmesh/shell_matrix.h"
27 : #include "libmesh/tensor_tools.h"
28 : #include "libmesh/enum_solver_package.h"
29 : #include "libmesh/int_range.h"
30 :
31 : // gzstream for reading/writing compressed files as a stream
32 : #ifdef LIBMESH_HAVE_GZSTREAM
33 : # include "gzstream.h"
34 : #endif
35 :
36 : // C++ includes
37 : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
38 : #include <cmath> // for std::abs
39 : #include <fstream>
40 : #include <iostream>
41 : #include <sstream>
42 : #include <limits>
43 : #include <memory>
44 :
45 : #ifdef LIBMESH_HAVE_CXX11_REGEX
46 : #include <regex>
47 : #endif
48 :
49 :
50 : namespace libMesh
51 : {
52 :
53 :
54 :
55 : //------------------------------------------------------------------
56 : // NumericVector methods
57 :
58 : // Full specialization for Real datatypes
59 : template <typename T>
60 : std::unique_ptr<NumericVector<T>>
61 1869824 : NumericVector<T>::build(const Parallel::Communicator & comm,
62 : SolverPackage solver_package,
63 : ParallelType parallel_type)
64 : {
65 : // Build the appropriate vector
66 1869824 : switch (solver_package)
67 : {
68 :
69 : #ifdef LIBMESH_HAVE_LASPACK
70 0 : case LASPACK_SOLVERS:
71 0 : return std::make_unique<LaspackVector<T>>(comm, parallel_type);
72 : #endif
73 :
74 : #ifdef LIBMESH_HAVE_PETSC
75 1857081 : case PETSC_SOLVERS:
76 1802289 : return std::make_unique<PetscVector<T>>(comm, parallel_type);
77 : #endif
78 :
79 : #ifdef LIBMESH_TRILINOS_HAVE_EPETRA
80 : case TRILINOS_SOLVERS:
81 : return std::make_unique<EpetraVector<T>>(comm, parallel_type);
82 : #endif
83 :
84 : #ifdef LIBMESH_HAVE_EIGEN
85 12743 : case EIGEN_SOLVERS:
86 12743 : return std::make_unique<EigenSparseVector<T>>(comm, parallel_type);
87 : #endif
88 :
89 0 : default:
90 0 : return std::make_unique<DistributedVector<T>>(comm, parallel_type);
91 : }
92 : }
93 :
94 :
95 :
96 : template <typename T>
97 69 : void NumericVector<T>::set_type(ParallelType t)
98 : {
99 : // Check for no-op
100 69 : if (_type == t)
101 0 : return;
102 :
103 : // If the NumericVector is not yet initialized, then it is generally
104 : // safe to change the ParallelType, with minor restrictions.
105 69 : if (!this->initialized())
106 : {
107 : // If ghosted vectors are not enabled and the user requested a
108 : // GHOSTED vector, fall back on SERIAL.
109 : #ifndef LIBMESH_ENABLE_GHOSTED
110 : if (t == GHOSTED)
111 : {
112 : _type = SERIAL;
113 : return;
114 : }
115 : #endif
116 :
117 69 : _type = t;
118 69 : return;
119 : }
120 :
121 : // If we made it here, then the NumericVector was already
122 : // initialized and we don't currently allow the ParallelType to be
123 : // changed, although this could potentially be added later.
124 0 : libmesh_not_implemented();
125 : }
126 :
127 : template <typename T>
128 25516 : void NumericVector<T>::insert (const T * v,
129 : const std::vector<numeric_index_type> & dof_indices)
130 : {
131 256 : libmesh_assert (v);
132 :
133 188142 : for (auto i : index_range(dof_indices))
134 162992 : this->set (dof_indices[i], v[i]);
135 25516 : }
136 :
137 :
138 :
139 : template <typename T>
140 0 : void NumericVector<T>::insert (const NumericVector<T> & V,
141 : const std::vector<numeric_index_type> & dof_indices)
142 : {
143 0 : libmesh_assert_equal_to (V.size(), dof_indices.size());
144 0 : libmesh_assert (V.readable());
145 :
146 0 : for (auto i : index_range(dof_indices))
147 0 : this->set (dof_indices[i], V(i));
148 0 : }
149 :
150 :
151 :
152 : template <typename T>
153 0 : int NumericVector<T>::compare (const NumericVector<T> & other_vector,
154 : const Real threshold) const
155 : {
156 0 : libmesh_assert(this->compatible(other_vector));
157 :
158 0 : int first_different_i = std::numeric_limits<int>::max();
159 0 : numeric_index_type i = first_local_index();
160 :
161 0 : while (first_different_i==std::numeric_limits<int>::max()
162 0 : && i<last_local_index())
163 : {
164 0 : if (std::abs((*this)(i) - other_vector(i)) > threshold)
165 0 : first_different_i = i;
166 : else
167 0 : i++;
168 : }
169 :
170 : // Find the correct first differing index in parallel
171 0 : this->comm().min(first_different_i);
172 :
173 0 : if (first_different_i == std::numeric_limits<int>::max())
174 0 : return -1;
175 :
176 0 : return first_different_i;
177 : }
178 :
179 :
180 : template <typename T>
181 0 : int NumericVector<T>::local_relative_compare (const NumericVector<T> & other_vector,
182 : const Real threshold) const
183 : {
184 0 : libmesh_assert(this->compatible(other_vector));
185 :
186 0 : int first_different_i = std::numeric_limits<int>::max();
187 0 : numeric_index_type i = first_local_index();
188 :
189 0 : do
190 : {
191 0 : if (std::abs((*this)(i) - other_vector(i)) > threshold *
192 0 : std::max(std::abs((*this)(i)), std::abs(other_vector(i))))
193 0 : first_different_i = i;
194 : else
195 0 : i++;
196 : }
197 0 : while (first_different_i==std::numeric_limits<int>::max()
198 0 : && i<last_local_index());
199 :
200 : // Find the correct first differing index in parallel
201 0 : this->comm().min(first_different_i);
202 :
203 0 : if (first_different_i == std::numeric_limits<int>::max())
204 0 : return -1;
205 :
206 0 : return first_different_i;
207 : }
208 :
209 :
210 : template <typename T>
211 0 : int NumericVector<T>::global_relative_compare (const NumericVector<T> & other_vector,
212 : const Real threshold) const
213 : {
214 0 : libmesh_assert(this->compatible(other_vector));
215 :
216 0 : int first_different_i = std::numeric_limits<int>::max();
217 0 : numeric_index_type i = first_local_index();
218 :
219 0 : const Real my_norm = this->linfty_norm();
220 0 : const Real other_norm = other_vector.linfty_norm();
221 0 : const Real abs_threshold = std::max(my_norm, other_norm) * threshold;
222 :
223 0 : do
224 : {
225 0 : if (std::abs((*this)(i) - other_vector(i) ) > abs_threshold)
226 0 : first_different_i = i;
227 : else
228 0 : i++;
229 : }
230 0 : while (first_different_i==std::numeric_limits<int>::max()
231 0 : && i<last_local_index());
232 :
233 : // Find the correct first differing index in parallel
234 0 : this->comm().min(first_different_i);
235 :
236 0 : if (first_different_i == std::numeric_limits<int>::max())
237 0 : return -1;
238 :
239 0 : return first_different_i;
240 : }
241 :
242 :
243 :
244 : template <typename T>
245 292 : void NumericVector<T>::print_matlab(const std::string & filename) const
246 : {
247 12 : parallel_object_only();
248 :
249 12 : libmesh_assert (this->initialized());
250 :
251 292 : const numeric_index_type first_dof = this->first_local_index(),
252 292 : end_dof = this->last_local_index();
253 :
254 : // We'll print the vector from processor 0 to make sure
255 : // it's serialized properly
256 304 : if (this->processor_id() == 0)
257 : {
258 164 : std::unique_ptr<std::ofstream> file;
259 :
260 174 : if (filename != "")
261 328 : file = std::make_unique<std::ofstream>(filename.c_str());
262 :
263 174 : std::ostream & os = (filename == "") ? libMesh::out : *file;
264 :
265 : // Print a header similar to PETSC_VIEWER_ASCII_MATLAB
266 174 : os << "%Vec Object: () " << this->n_processors() << " MPI processes\n"
267 30 : << "% type: " << (this->n_processors() > 1 ? "mpi" : "seq") << "\n"
268 330 : << "Vec = [\n";
269 :
270 1914 : for (numeric_index_type i : make_range(end_dof))
271 3340 : os << (*this)(i) << '\n';
272 :
273 20 : std::vector<T> cbuf;
274 292 : for (auto p : IntRange<processor_id_type>(1, this->n_processors()))
275 : {
276 118 : this->comm().receive(p, cbuf);
277 :
278 1920 : for (auto c : cbuf)
279 1802 : os << c << '\n';
280 : }
281 :
282 164 : os << "];\n" << std::endl;
283 154 : }
284 : else
285 : {
286 120 : std::vector<T> cbuf(end_dof - first_dof);
287 :
288 1920 : for (numeric_index_type i : make_range(end_dof - first_dof))
289 1802 : cbuf[i] = (*this)(first_dof+i);
290 118 : this->comm().send(0,cbuf);
291 : }
292 292 : }
293 :
294 :
295 :
296 : template <typename T>
297 428 : void NumericVector<T>::read_matlab(const std::string & filename)
298 : {
299 32 : LOG_SCOPE("read_matlab()", "NumericVector");
300 :
301 : #ifndef LIBMESH_HAVE_CXX11_REGEX
302 : libmesh_not_implemented(); // What is your compiler?!? Email us!
303 : libmesh_ignore(filename);
304 : #else
305 16 : parallel_object_only();
306 :
307 428 : const bool gzipped_file = Utility::ends_with(filename, ".gz");
308 :
309 : // If we don't already have this size, we'll need to reinit, and
310 : // determine which entries each processor is in charge of.
311 32 : std::vector<numeric_index_type> first_entries, end_entries;
312 :
313 428 : numeric_index_type first_entry = 0,
314 428 : end_entry = 0,
315 428 : n = 0;
316 :
317 : // We'll use an istream here; it might be an ifstream if we're
318 : // opening a raw ASCII file or a gzstream if we're opening a
319 : // compressed one.
320 412 : std::unique_ptr<std::istream> file;
321 :
322 : // First read through the file, saving size and entry data
323 32 : std::vector<T> entries;
324 :
325 : {
326 : // We'll read the vector on processor 0 rather than try to juggle
327 : // parallel I/O.
328 444 : if (this->processor_id() == 0)
329 : {
330 : // We'll be using regular expressions to make ourselves slightly
331 : // more robust to formatting.
332 216 : const std::regex start_regex // assignment like "Vec_0x84000002_1 = ["
333 : ("\\s*\\w+\\s*=\\s*\\[");
334 216 : const std::regex end_regex // end of assignment
335 : ("^[^%]*\\]");
336 :
337 192 : if (gzipped_file)
338 : {
339 : #ifdef LIBMESH_HAVE_GZSTREAM
340 0 : auto inf = std::make_unique<igzstream>();
341 0 : libmesh_assert(inf);
342 0 : inf->open(filename.c_str(), std::ios::in);
343 0 : file = std::move(inf);
344 : #else
345 : libmesh_error_msg("ERROR: need gzstream to handle .gz files!!!");
346 : #endif
347 0 : }
348 : else
349 : {
350 204 : auto inf = std::make_unique<std::ifstream>();
351 12 : libmesh_assert(inf);
352 :
353 204 : std::string new_name = Utility::unzip_file(filename);
354 :
355 192 : inf->open(new_name.c_str(), std::ios::in);
356 180 : file = std::move(inf);
357 168 : }
358 :
359 204 : const std::string whitespace = " \t";
360 :
361 12 : bool have_started = false;
362 12 : bool have_ended = false;
363 :
364 12172 : for (std::string line; std::getline(*file, line);)
365 : {
366 212 : std::smatch sm;
367 :
368 : // First, try to match an entry. This is the most common
369 : // case so we won't rely on slow std::regex for it.
370 : // stringstream is at least an improvement over that.
371 6492 : std::istringstream l(line);
372 0 : T value;
373 212 : l >> value;
374 :
375 6292 : if (!l.fail())
376 : {
377 5524 : libmesh_error_msg_if
378 : (!have_started, "Confused by premature entries in vector file " << filename);
379 :
380 5524 : entries.push_back(value);
381 5524 : ++n;
382 : }
383 :
384 768 : else if (std::regex_search(line, start_regex))
385 12 : have_started = true;
386 :
387 576 : else if (std::regex_search(line, end_regex))
388 : {
389 12 : have_ended = true;
390 24 : break;
391 : }
392 : }
393 :
394 192 : libmesh_error_msg_if
395 : (!have_started, "Confused by missing assignment-beginning in vector file " << filename);
396 :
397 192 : libmesh_error_msg_if
398 : (!have_ended, "Confused by missing assignment-ending in vector file " << filename);
399 168 : }
400 :
401 428 : this->comm().broadcast(n);
402 :
403 444 : if (this->initialized() &&
404 428 : n == this->size())
405 : {
406 444 : first_entry = this->first_local_index(),
407 428 : end_entry = this->last_local_index();
408 : }
409 : else
410 : {
411 : // Determine which rows/columns each processor will be in charge of
412 0 : first_entry = this->processor_id() * n / this->n_processors(),
413 0 : end_entry = (this->processor_id()+1) * n / this->n_processors();
414 : }
415 :
416 428 : this->comm().gather(0, first_entry, first_entries);
417 428 : this->comm().gather(0, end_entry, end_entries);
418 :
419 : } // Done reading entry data and broadcasting vector size
420 :
421 : // If we're not already initialized compatibly with the file then
422 : // we'll initialize here
423 824 : bool need_init = !this->initialized() ||
424 444 : (this->size() != n) ||
425 428 : (this->local_size() != end_entry - first_entry);
426 :
427 428 : this->comm().max(need_init);
428 :
429 428 : if (need_init)
430 0 : this->init(n, end_entry - first_entry);
431 :
432 : // Set the vector values last. The iota call here is inefficient
433 : // but it's probably better than calling a single virtual function
434 : // per index.
435 444 : if (this->processor_id() == 0)
436 620 : for (auto p : make_range(this->n_processors()))
437 : {
438 428 : const numeric_index_type first_entry_p = first_entries[p];
439 428 : const numeric_index_type n_local = end_entries[p] - first_entry_p;
440 444 : std::vector<numeric_index_type> indices(n_local);
441 16 : std::iota(indices.begin(), indices.end(), first_entry_p);
442 444 : this->insert(entries.data() + first_entries[p],
443 : indices);
444 : }
445 :
446 428 : this->close();
447 : #endif
448 824 : }
449 :
450 :
451 : /*
452 : // Full specialization for float datatypes (DistributedVector wants this)
453 :
454 : template <>
455 : int NumericVector<float>::compare (const NumericVector<float> & other_vector,
456 : const Real threshold) const
457 : {
458 : libmesh_assert (this->initialized());
459 : libmesh_assert (other_vector.initialized());
460 : libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
461 : libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
462 :
463 : int rvalue = -1;
464 : numeric_index_type i = first_local_index();
465 :
466 : do
467 : {
468 : if (std::abs((*this)(i) - other_vector(i) ) > threshold)
469 : rvalue = i;
470 : else
471 : i++;
472 : }
473 : while (rvalue==-1 && i<last_local_index());
474 :
475 : return rvalue;
476 : }
477 :
478 : // Full specialization for double datatypes
479 : template <>
480 : int NumericVector<double>::compare (const NumericVector<double> & other_vector,
481 : const Real threshold) const
482 : {
483 : libmesh_assert (this->initialized());
484 : libmesh_assert (other_vector.initialized());
485 : libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
486 : libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
487 :
488 : int rvalue = -1;
489 : numeric_index_type i = first_local_index();
490 :
491 : do
492 : {
493 : if (std::abs((*this)(i) - other_vector(i) ) > threshold)
494 : rvalue = i;
495 : else
496 : i++;
497 : }
498 : while (rvalue==-1 && i<last_local_index());
499 :
500 : return rvalue;
501 : }
502 :
503 : #ifdef LIBMESH_DEFAULT_TRIPLE_PRECISION
504 : // Full specialization for long double datatypes
505 : template <>
506 : int NumericVector<long double>::compare (const NumericVector<long double> & other_vector,
507 : const Real threshold) const
508 : {
509 : libmesh_assert (this->initialized());
510 : libmesh_assert (other_vector.initialized());
511 : libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
512 : libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
513 :
514 : int rvalue = -1;
515 : numeric_index_type i = first_local_index();
516 :
517 : do
518 : {
519 : if (std::abs((*this)(i) - other_vector(i) ) > threshold)
520 : rvalue = i;
521 : else
522 : i++;
523 : }
524 : while (rvalue==-1 && i<last_local_index());
525 :
526 : return rvalue;
527 : }
528 : #endif
529 :
530 :
531 : // Full specialization for Complex datatypes
532 : template <>
533 : int NumericVector<Complex>::compare (const NumericVector<Complex> & other_vector,
534 : const Real threshold) const
535 : {
536 : libmesh_assert (this->initialized());
537 : libmesh_assert (other_vector.initialized());
538 : libmesh_assert_equal_to (this->first_local_index(), other_vector.first_local_index());
539 : libmesh_assert_equal_to (this->last_local_index(), other_vector.last_local_index());
540 :
541 : int rvalue = -1;
542 : numeric_index_type i = first_local_index();
543 :
544 : do
545 : {
546 : if ((std::abs((*this)(i).real() - other_vector(i).real()) > threshold) || (std::abs((*this)(i).imag() - other_vector(i).imag()) > threshold))
547 : rvalue = i;
548 : else
549 : i++;
550 : }
551 : while (rvalue==-1 && i<this->last_local_index());
552 :
553 : return rvalue;
554 : }
555 : */
556 :
557 :
558 : template <class T>
559 0 : Real NumericVector<T>::subset_l1_norm (const std::set<numeric_index_type> & indices) const
560 : {
561 0 : libmesh_assert (this->readable());
562 :
563 0 : const NumericVector<T> & v = *this;
564 :
565 0 : Real norm = 0;
566 :
567 0 : for (const auto & index : indices)
568 0 : norm += std::abs(v(index));
569 :
570 0 : this->comm().sum(norm);
571 :
572 0 : return norm;
573 : }
574 :
575 : template <class T>
576 0 : Real NumericVector<T>::subset_l2_norm (const std::set<numeric_index_type> & indices) const
577 : {
578 0 : libmesh_assert (this->readable());
579 :
580 0 : const NumericVector<T> & v = *this;
581 :
582 0 : Real norm = 0;
583 :
584 0 : for (const auto & index : indices)
585 0 : norm += TensorTools::norm_sq(v(index));
586 :
587 0 : this->comm().sum(norm);
588 :
589 0 : return std::sqrt(norm);
590 : }
591 :
592 : template <class T>
593 0 : Real NumericVector<T>::subset_linfty_norm (const std::set<numeric_index_type> & indices) const
594 : {
595 0 : libmesh_assert (this->readable());
596 :
597 0 : const NumericVector<T> & v = *this;
598 :
599 0 : Real norm = 0;
600 :
601 0 : for (const auto & index : indices)
602 : {
603 0 : Real value = std::abs(v(index));
604 0 : if (value > norm)
605 0 : norm = value;
606 : }
607 :
608 0 : this->comm().max(norm);
609 :
610 0 : return norm;
611 : }
612 :
613 :
614 :
615 : template <class T>
616 432 : Real NumericVector<T>::l2_norm_diff (const NumericVector<T> & v) const
617 : {
618 16 : libmesh_assert(this->compatible(v));
619 :
620 432 : Real norm = 0;
621 5996 : for (const auto i : make_range(this->first_local_index(), this->last_local_index()))
622 5564 : norm += TensorTools::norm_sq((*this)(i) - v(i));
623 :
624 432 : this->comm().sum(norm);
625 :
626 448 : return std::sqrt(norm);
627 : }
628 :
629 :
630 :
631 : template <class T>
632 1502 : Real NumericVector<T>::l1_norm_diff (const NumericVector<T> & v) const
633 : {
634 54 : libmesh_assert(this->compatible(v));
635 :
636 1502 : Real norm = 0;
637 27867 : for (const auto i : make_range(this->first_local_index(), this->last_local_index()))
638 26365 : norm += libMesh::l1_norm_diff((*this)(i), v(i));
639 :
640 1502 : this->comm().sum(norm);
641 :
642 1502 : return norm;
643 : }
644 :
645 :
646 :
647 : template <typename T>
648 1201601 : void NumericVector<T>::add_vector (const T * v,
649 : const std::vector<numeric_index_type> & dof_indices)
650 : {
651 0 : libmesh_assert(v);
652 :
653 12365358 : for (auto i : index_range(dof_indices))
654 11163757 : this->add (dof_indices[i], v[i]);
655 1201601 : }
656 :
657 :
658 :
659 : template <typename T>
660 0 : void NumericVector<T>::add_vector (const NumericVector<T> & v,
661 : const std::vector<numeric_index_type> & dof_indices)
662 : {
663 0 : libmesh_assert(v.readable());
664 :
665 0 : const std::size_t n = dof_indices.size();
666 0 : libmesh_assert_equal_to(v.size(), n);
667 0 : for (numeric_index_type i=0; i != n; i++)
668 0 : this->add (dof_indices[i], v(i));
669 0 : }
670 :
671 :
672 :
673 : template <typename T>
674 0 : void NumericVector<T>::add_vector (const NumericVector<T> & v,
675 : const ShellMatrix<T> & a)
676 : {
677 0 : libmesh_assert(this->compatible(v));
678 :
679 0 : a.vector_mult_add(*this,v);
680 0 : }
681 :
682 :
683 :
684 : template <typename T>
685 140 : bool NumericVector<T>::readable () const
686 : {
687 140 : return this->initialized() && this->closed();
688 : }
689 :
690 :
691 : template <typename T>
692 70 : bool NumericVector<T>::compatible (const NumericVector<T> & v) const
693 : {
694 140 : return this->readable() && v.readable() &&
695 70 : this->size() == v.size() &&
696 70 : this->local_size() == v.local_size() &&
697 210 : this->first_local_index() == v.first_local_index() &&
698 140 : this->last_local_index() == v.last_local_index();
699 : }
700 :
701 :
702 : //------------------------------------------------------------------
703 : // Explicit instantiations
704 : template class LIBMESH_EXPORT NumericVector<Number>;
705 :
706 : } // namespace libMesh
|
