double_vector.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26#include "double_vector.h"
27#include "matrices.h"
28
29
30namespace oomph
31{
32 //============================================================================
33 /// Just copys the argument DoubleVector
34 //============================================================================
35 void DoubleVector::build(const DoubleVector& old_vector)
36 {
37 if (!(*this == old_vector))
38 {
39 // the vector owns the internal data
40 Internal_values = true;
41
42 // reset the distribution and resize the data
43 this->build(old_vector.distribution_pt(), 0.0);
44
45 // copy the data
46 if (this->distribution_built())
47 {
48 unsigned nrow_local = this->nrow_local();
49 const double* old_vector_values = old_vector.values_pt();
50 std::copy(old_vector_values, old_vector_values + nrow_local, Values_pt);
51 }
52 }
53 }
54
55 //============================================================================
56 /// Assembles a DoubleVector with distribution dist, if v is specified
57 /// each row is set to v
58 //============================================================================
60 const double& v)
61 {
62 // clean the memory
63 this->clear();
64
65 // the vector owns the internal data
66 Internal_values = true;
67
68 // Set the distribution
69 this->build_distribution(dist_pt);
70
71 // update the values
72 if (dist_pt->built())
73 {
74 unsigned nrow_local = this->nrow_local();
75 Values_pt = new double[nrow_local];
76
77 std::fill_n(Values_pt, nrow_local, v);
78 Built = true;
79 }
80 else
81 {
82 Built = false;
83 }
84 }
85
86 //============================================================================
87 /// Assembles a DoubleVector with a distribution dist and coefficients
88 /// taken from the vector v.
89 /// Note. The vector v MUST be of length nrow()
90 //============================================================================
92 const Vector<double>& v)
93 {
94 // clean the memory
95 this->clear();
96
97 // the vector owns the internal data
98 Internal_values = true;
99
100 // Set the distribution
101 this->build_distribution(dist_pt);
102
103 // update the values
104 if (dist_pt->built())
105 {
106 // re-allocate memory which was deleted by clear()
107 unsigned nrow_local = this->nrow_local();
108 Values_pt = new double[nrow_local];
109
110 // use the initialise method to populate the vector
111 this->initialise(v);
112 Built = true;
113 }
114 else
115 {
116 Built = false;
117 }
118 }
119
120 //============================================================================
121 /// initialise the whole vector with value v
122 //============================================================================
123 void DoubleVector::initialise(const double& v)
124 {
125 if (Built)
126 {
127 // cache nrow local
128 unsigned nrow_local = this->nrow_local();
129
130 std::fill_n(Values_pt, nrow_local, v);
131 }
132 }
133
134 //============================================================================
135 /// initialise the vector with coefficient from the vector v.
136 /// Note: The vector v must be of length
137 //============================================================================
139 {
140#ifdef PARANOID
141 if (v.size() != this->nrow())
142 {
143 std::ostringstream error_message;
144 error_message << "The vector passed to initialise(...) must be of length "
145 << "nrow()";
146 throw OomphLibError(
147 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
148 }
149#endif
150 unsigned begin_first_row = this->first_row();
151 unsigned end = begin_first_row + this->nrow_local();
152
153 std::copy(v.begin() + begin_first_row, v.begin() + end, Values_pt);
154 }
155
156 //============================================================================
157 /// The contents of the vector are redistributed to match the new
158 /// distribution. In a non-MPI build this method works, but does nothing.
159 /// \b NOTE 1: The current distribution and the new distribution must have
160 /// the same number of global rows.
161 /// \b NOTE 2: The current distribution and the new distribution must have
162 /// the same Communicator.
163 //============================================================================
165 const LinearAlgebraDistribution* const& dist_pt)
166 {
167#ifdef OOMPH_HAS_MPI
168#ifdef PARANOID
169 if (!Internal_values)
170 {
171 // if this vector does not own the double* values then it cannot be
172 // distributed.
173 // note: this is not stictly necessary - would just need to be careful
174 // with delete[] below.
175 std::ostringstream error_message;
176 error_message << "This vector does not own its data (i.e. it has been "
177 << "passed in via set_external_values() and therefore "
178 << "cannot be redistributed";
179 throw OomphLibError(
180 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
181 }
182 // paranoid check that the nrows for both distributions is the
183 // same
184 if (dist_pt->nrow() != this->nrow())
185 {
186 std::ostringstream error_message;
187 error_message << "The number of global rows in the new distribution ("
188 << dist_pt->nrow() << ") is not equal to the number"
189 << " of global rows in the current distribution ("
190 << this->nrow() << ").\n";
191 throw OomphLibError(
192 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
193 }
194 // paranoid check that the current distribution and the new distribution
195 // have the same Communicator
196 OomphCommunicator temp_comm(*dist_pt->communicator_pt());
197 if (!(temp_comm == *this->distribution_pt()->communicator_pt()))
198 {
199 std::ostringstream error_message;
200 error_message << "The new distribution and the current distribution must "
201 << "have the same communicator.";
202 throw OomphLibError(
203 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
204 }
205#endif
206
207 // check the distributions are not the same
208 if (!((*this->distribution_pt()) == *dist_pt))
209 {
210 // get the rank and the number of processors
211 int my_rank = this->distribution_pt()->communicator_pt()->my_rank();
212 int nproc = this->distribution_pt()->communicator_pt()->nproc();
213
214 // if both vectors are distributed
215 if (this->distributed() && dist_pt->distributed())
216 {
217 // new nrow_local and first_row data
218 Vector<unsigned> new_first_row_data(nproc);
219 Vector<unsigned> new_nrow_local_data(nproc);
220 Vector<unsigned> current_first_row_data(nproc);
221 Vector<unsigned> current_nrow_local_data(nproc);
222 for (int i = 0; i < nproc; i++)
223 {
224 new_first_row_data[i] = dist_pt->first_row(i);
225 new_nrow_local_data[i] = dist_pt->nrow_local(i);
226 current_first_row_data[i] = this->first_row(i);
227 current_nrow_local_data[i] = this->nrow_local(i);
228 }
229
230 // compute which local rows are expected to be received from each
231 // processor / sent to each processor
232 Vector<unsigned> new_first_row_for_proc(nproc);
233 Vector<unsigned> new_nrow_local_for_proc(nproc);
234 Vector<unsigned> new_first_row_from_proc(nproc);
235 Vector<unsigned> new_nrow_local_from_proc(nproc);
236
237 // for every processor compute first_row and nrow_local that will
238 // will sent and received by this processor
239 for (int p = 0; p < nproc; p++)
240 {
241 // start with data to be sent
242 if ((new_first_row_data[p] < (current_first_row_data[my_rank] +
243 current_nrow_local_data[my_rank])) &&
244 (current_first_row_data[my_rank] <
245 (new_first_row_data[p] + new_nrow_local_data[p])))
246 {
247 new_first_row_for_proc[p] =
248 std::max(current_first_row_data[my_rank], new_first_row_data[p]);
249 new_nrow_local_for_proc[p] =
250 std::min((current_first_row_data[my_rank] +
251 current_nrow_local_data[my_rank]),
252 (new_first_row_data[p] + new_nrow_local_data[p])) -
253 new_first_row_for_proc[p];
254 }
255
256 // and data to be received
257 if ((new_first_row_data[my_rank] <
258 (current_first_row_data[p] + current_nrow_local_data[p])) &&
259 (current_first_row_data[p] <
260 (new_first_row_data[my_rank] + new_nrow_local_data[my_rank])))
261 {
262 new_first_row_from_proc[p] =
263 std::max(current_first_row_data[p], new_first_row_data[my_rank]);
264 new_nrow_local_from_proc[p] =
265 std::min(
266 (current_first_row_data[p] + current_nrow_local_data[p]),
267 (new_first_row_data[my_rank] + new_nrow_local_data[my_rank])) -
268 new_first_row_from_proc[p];
269 }
270 }
271
272 // temporary storage for the new data
273 double* temp_data = new double[new_nrow_local_data[my_rank]];
274
275 // "send to self" or copy Data that does not need to be sent else where
276 // to temp_data
277 if (new_nrow_local_for_proc[my_rank] != 0)
278 {
279 unsigned j =
280 new_first_row_for_proc[my_rank] - current_first_row_data[my_rank];
281 unsigned k =
282 new_first_row_for_proc[my_rank] - new_first_row_data[my_rank];
283 for (unsigned i = 0; i < new_nrow_local_for_proc[my_rank]; i++)
284 {
285 temp_data[k + i] = Values_pt[j + i];
286 }
287 }
288
289 // send and receive circularly
290 for (int p = 1; p < nproc; p++)
291 {
292 // next processor to send to
293 unsigned dest_p = (my_rank + p) % nproc;
294
295 // next processor to receive from
296 unsigned source_p = (nproc + my_rank - p) % nproc;
297
298 // send and receive the value
299 MPI_Status status;
300 MPI_Sendrecv(Values_pt + new_first_row_for_proc[dest_p] -
301 current_first_row_data[my_rank],
302 new_nrow_local_for_proc[dest_p],
303 MPI_DOUBLE,
304 dest_p,
305 1,
306 temp_data + new_first_row_from_proc[source_p] -
307 new_first_row_data[my_rank],
308 new_nrow_local_from_proc[source_p],
309 MPI_DOUBLE,
310 source_p,
311 1,
312 this->distribution_pt()->communicator_pt()->mpi_comm(),
313 &status);
314 }
315
316 // copy from temp data to Values_pt
317 delete[] Values_pt;
318 unsigned nrow_local = dist_pt->nrow_local();
319 Values_pt = new double[nrow_local];
320 for (unsigned i = 0; i < nrow_local; i++)
321 {
322 Values_pt[i] = temp_data[i];
323 }
324 delete[] temp_data;
325 }
326
327 // if this vector is distributed but the new distributed is global
328 else if (this->distributed() && !dist_pt->distributed())
329 {
330 // copy existing Values_pt to temp_data
331 unsigned nrow_local = this->nrow_local();
332 double* temp_data = new double[nrow_local];
333 for (unsigned i = 0; i < nrow_local; i++)
334 {
335 temp_data[i] = Values_pt[i];
336 }
337
338 // clear and resize Values_pt
339 delete[] Values_pt;
340 Values_pt = new double[this->nrow()];
341
342 // create a int vector of first rows
343 int* dist_first_row = new int[nproc];
344 int* dist_nrow_local = new int[nproc];
345 for (int p = 0; p < nproc; p++)
346 {
347 dist_first_row[p] = this->first_row(p);
348 dist_nrow_local[p] = this->nrow_local(p);
349 }
350
351 // gather the local vectors from all processors on all processors
352 int my_nrow_local(this->nrow_local());
353 MPI_Allgatherv(temp_data,
354 my_nrow_local,
355 MPI_DOUBLE,
356 Values_pt,
357 dist_nrow_local,
358 dist_first_row,
359 MPI_DOUBLE,
360 this->distribution_pt()->communicator_pt()->mpi_comm());
361
362 // update the distribution
363 this->build_distribution(dist_pt);
364
365 // delete the temp_data
366 delete[] temp_data;
367
368 // clean up
369 delete[] dist_first_row;
370 delete[] dist_nrow_local;
371 }
372
373 // if this vector is not distrubted but the target vector is
374 else if (!this->distributed() && dist_pt->distributed())
375 {
376 // cache the new nrow_local
377 unsigned nrow_local = dist_pt->nrow_local();
378
379 // and first_row
380 unsigned first_row = dist_pt->first_row();
381
382 // temp storage for the new data
383 double* temp_data = new double[nrow_local];
384
385 // copy the data
386 for (unsigned i = 0; i < nrow_local; i++)
387 {
388 temp_data[i] = Values_pt[first_row + i];
389 }
390
391 // copy to Values_pt
392 delete[] Values_pt;
393 Values_pt = temp_data;
394
395 // update the distribution
396 this->build_distribution(dist_pt);
397 }
398
399 // copy the Distribution
400 this->build_distribution(dist_pt);
401 }
402#endif
403 }
404
405 //============================================================================
406 /// [] access function to the (local) values of this vector
407 //============================================================================
409 {
410#ifdef RANGE_CHECKING
411 if (i >= int(this->nrow_local()))
412 {
413 std::ostringstream error_message;
414 error_message << "Range Error: " << i << " is not in the range (0,"
415 << this->nrow_local() - 1 << ")";
416 throw OomphLibError(
417 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
418 }
419#endif
420 return Values_pt[i];
421 }
422
423 //============================================================================
424 /// == operator
425 //============================================================================
427 {
428 // if v is not setup return false
429 if (v.built() && !this->built())
430 {
431 return false;
432 }
433 else if (!v.built() && this->built())
434 {
435 return false;
436 }
437 else if (!v.built() && !this->built())
438 {
439 return true;
440 }
441 else
442 {
443 const double* v_values_pt = v.values_pt();
444 unsigned nrow_local = this->nrow_local();
445 for (unsigned i = 0; i < nrow_local; i++)
446 {
447 if (Values_pt[i] != v_values_pt[i])
448 {
449 return false;
450 }
451 }
452 return true;
453 }
454 }
455
456 //============================================================================
457 /// += operator
458 //============================================================================
460 {
461#ifdef PARANOID
462 // PARANOID check that this vector is setup
463 if (!this->built())
464 {
465 std::ostringstream error_message;
466 error_message << "This vector must be setup.";
467 throw OomphLibError(
468 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
469 }
470 // PARANOID check that the vector v is setup
471 if (!v.built())
472 {
473 std::ostringstream error_message;
474 error_message << "The vector v must be setup.";
475 throw OomphLibError(
476 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
477 }
478 // PARANOID check that the vectors have the same distribution
479 if (!(*v.distribution_pt() == *this->distribution_pt()))
480 {
481 std::ostringstream error_message;
482 error_message << "The vector v and this vector must have the same "
483 << "distribution.";
484 throw OomphLibError(
485 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
486 }
487#endif //
488
489 // cache nrow_local
490 double* v_values_pt = v.values_pt();
491 unsigned nrow_local = this->nrow_local();
492
493 // Decided to keep this as a loop rather than use std::transform, because
494 // this is a very simple loop and should compile to the same code.
495 for (unsigned i = 0; i < nrow_local; i++)
496 {
497 Values_pt[i] += v_values_pt[i];
498 }
499 }
500
501 //============================================================================
502 /// -= operator
503 //============================================================================
505 {
506#ifdef PARANOID
507 // PARANOID check that this vector is setup
508 if (!this->distribution_built())
509 {
510 std::ostringstream error_message;
511 error_message << "This vector must be setup.";
512 throw OomphLibError(
513 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
514 }
515 // PARANOID check that the vector v is setup
516 if (!v.built())
517 {
518 std::ostringstream error_message;
519 error_message << "The vector v must be setup.";
520 throw OomphLibError(
521 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
522 }
523 // PARANOID check that the vectors have the same distribution
524 if (!(*v.distribution_pt() == *this->distribution_pt()))
525 {
526 std::ostringstream error_message;
527 error_message << "The vector v and this vector must have the same "
528 << "distribution.";
529 throw OomphLibError(
530 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
531 }
532#endif
533
534 // cache nrow_local
535 double* v_values_pt = v.values_pt();
536 unsigned nrow_local = this->nrow_local();
537
538 // Decided to keep this as a loop rather than use std::transform, because
539 // this is a very simple loop and should compile to the same code.
540 for (unsigned i = 0; i < nrow_local; i++)
541 {
542 Values_pt[i] -= v_values_pt[i];
543 }
544 }
545
546
547 //============================================================================
548 /// Multiply by double
549 //============================================================================
550 void DoubleVector::operator*=(const double& d)
551 {
552#ifdef PARANOID
553 if (!this->distribution_built())
554 {
555 std::ostringstream error_msg;
556 error_msg << "DoubleVector must be set up.";
557 throw OomphLibError(
558 error_msg.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
559 }
560#endif
561
562 // Decided to keep this as a loop rather than use std::transform, because
563 // this is a very simple loop and should compile to the same code.
564 for (unsigned i = 0, ni = this->nrow_local(); i < ni; i++)
565 {
566 Values_pt[i] *= d;
567 }
568 }
569
570 //============================================================================
571 /// Divide by double
572 //============================================================================
573 void DoubleVector::operator/=(const double& d)
574 {
575 // PARANOID checks are done inside operator *=
576
577 // Decided to keep this as a loop rather than use std::transform, because
578 // this is a very simple loop and should compile to the same code.
579 double divisor = (1.0 / d);
580 this->operator*=(divisor);
581 }
582
583 //============================================================================
584 /// [] access function to the (local) values of this vector
585 //============================================================================
586 const double& DoubleVector::operator[](int i) const
587 {
588#ifdef RANGE_CHECKING
589 if (i >= int(this->nrow_local()))
590 {
591 std::ostringstream error_message;
592 error_message << "Range Error: " << i << " is not in the range (0,"
593 << this->nrow_local() - 1 << ")";
594 throw OomphLibError(
595 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
596 }
597#endif
598 return Values_pt[i];
599 }
600
601 //============================================================================
602 /// returns the maximum coefficient
603 //============================================================================
604 double DoubleVector::max() const
605 {
606 // the number of local rows
607 unsigned nrow = this->nrow_local();
608
609 // get the local maximum
610 double max = 0.0;
611 for (unsigned i = 0; i < nrow; i++)
612 {
613 if (std::fabs(Values_pt[i]) > std::fabs(max))
614 {
615 max = std::fabs(Values_pt[i]);
616 }
617 }
618
619 // now return the maximum
620#ifdef OOMPH_HAS_MPI
621 // if this vector is not distributed then the local maximum is the global
622 // maximum
623 if (!this->distributed())
624 {
625 return max;
626 }
627 // else if the vector is distributed but only on a single processor
628 // then the local maximum is the global maximum
629 else if (this->distribution_pt()->communicator_pt()->nproc() == 1)
630 {
631 return max;
632 }
633 // otherwise use MPI_Allreduce to find the global maximum
634 else
635 {
636 double local_max = max;
637 MPI_Allreduce(&local_max,
638 &max,
639 1,
640 MPI_DOUBLE,
641 MPI_MAX,
642 this->distribution_pt()->communicator_pt()->mpi_comm());
643 return max;
644 }
645#else
646 return max;
647#endif
648 }
649
650 //============================================================================
651 /// output the contents of the vector
652 //============================================================================
653 void DoubleVector::output(std::ostream& outfile,
654 const int& output_precision) const
655 {
656 // temp pointer to values
657 double* temp;
658
659 // number of global row
660 unsigned nrow = this->nrow();
661
662#ifdef OOMPH_HAS_MPI
663
664 // number of local rows
665 int nrow_local = this->nrow_local();
666
667 // gather from all processors
668 if (this->distributed() &&
669 this->distribution_pt()->communicator_pt()->nproc() > 1)
670 {
671 // number of processors
672 int nproc = this->distribution_pt()->communicator_pt()->nproc();
673
674 // number of gobal row
675 unsigned nrow = this->nrow();
676
677 // get the vector of first_row s and nrow_local s
678 int* dist_first_row = new int[nproc];
679 int* dist_nrow_local = new int[nproc];
680 for (int p = 0; p < nproc; p++)
681 {
682 dist_first_row[p] = this->first_row(p);
683 dist_nrow_local[p] = this->nrow_local(p);
684 }
685
686 // gather
687 temp = new double[nrow];
688 MPI_Allgatherv(Values_pt,
690 MPI_DOUBLE,
691 temp,
692 dist_nrow_local,
693 dist_first_row,
694 MPI_DOUBLE,
695 this->distribution_pt()->communicator_pt()->mpi_comm());
696
697 // clean up
698 delete[] dist_first_row;
699 delete[] dist_nrow_local;
700 }
701 else
702 {
703 temp = Values_pt;
704 }
705#else
706 temp = Values_pt;
707#endif
708
709 // output
710 // Store the precision so we can revert it.
711 std::streamsize old_precision = 0;
712 if (output_precision > 0)
713 {
714 old_precision = outfile.precision();
715 outfile << std::setprecision(output_precision);
716 }
717
718 for (unsigned i = 0; i < nrow; i++)
719 {
720 outfile << i << " " << temp[i] << std::endl;
721 }
722
723 // Revert the precision.
724 if (output_precision > 0)
725 {
726 outfile << std::setprecision(old_precision);
727 }
728
729 // clean up if requires
730#ifdef OOMPH_HAS_MPI
731 if (this->distributed() &&
732 this->distribution_pt()->communicator_pt()->nproc() > 1)
733 {
734 delete[] temp;
735 }
736#endif
737 }
738
739 //============================================================================
740 /// output the local contents of the vector
741 //============================================================================
742 void DoubleVector::output_local_values(std::ostream& outfile,
743 const int& output_precision) const
744 {
745 // Number of local rows.
746 unsigned nrow_local = this->nrow_local();
747
748 // output
749 // Store the precision so we can revert it.
750 std::streamsize old_precision = 0;
751 if (output_precision > 0)
752 {
753 old_precision = outfile.precision();
754 outfile << std::setprecision(output_precision);
755 }
756
757 for (unsigned i = 0; i < nrow_local; i++)
758 {
759 outfile << i << " " << Values_pt[i] << std::endl;
760 }
761
762 // Revert the precision.
763 if (output_precision > 0)
764 {
765 outfile << std::setprecision(old_precision);
766 }
767 }
768
769 //============================================================================
770 /// output the local contents of the vector with the first row offset.
771 //============================================================================
773 std::ostream& outfile, const int& output_precision) const
774 {
775 // Number of local rows.
776 unsigned nrow_local = this->nrow_local();
777
778 // First row on this processor.
779 unsigned first_row = this->first_row();
780
781 // output
782 // Store the precision so we can revert it.
783 std::streamsize old_precision = 0;
784 if (output_precision > 0)
785 {
786 old_precision = outfile.precision();
787 outfile << std::setprecision(output_precision);
788 }
789
790 for (unsigned i = 0; i < nrow_local; i++)
791 {
792 outfile << (i + first_row) << " " << Values_pt[i] << std::endl;
793 }
794
795 // Revert the precision.
796 if (output_precision > 0)
797 {
798 outfile << std::setprecision(old_precision);
799 }
800 }
801
802 //============================================================================
803 /// compute the dot product of this vector with the vector vec
804 //============================================================================
805 double DoubleVector::dot(const DoubleVector& vec) const
806 {
807#ifdef PARANOID
808 // paranoid check that the vector is setup
809 if (!this->built())
810 {
811 std::ostringstream error_message;
812 error_message << "This vector must be setup.";
813 throw OomphLibError(
814 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
815 }
816 if (!vec.built())
817 {
818 std::ostringstream error_message;
819 error_message << "The input vector be setup.";
820 throw OomphLibError(
821 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
822 }
823 if (*this->distribution_pt() != *vec.distribution_pt())
824 {
825 std::ostringstream error_message;
826 error_message << "The distribution of this vector and the vector vec "
827 << "must be the same."
828 << "\n\n this: " << *this->distribution_pt()
829 << "\n vec: " << *vec.distribution_pt();
830 throw OomphLibError(
831 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
832 }
833#endif
834
835 // compute the local norm
836 unsigned nrow_local = this->nrow_local();
837 double n = 0.0;
838 const double* vec_values_pt = vec.values_pt();
839 for (unsigned i = 0; i < nrow_local; i++)
840 {
841 n += Values_pt[i] * vec_values_pt[i];
842 }
843
844 // if this vector is distributed and on multiple processors then gather
845#ifdef OOMPH_HAS_MPI
846 double n2 = n;
847 if (this->distributed() &&
848 this->distribution_pt()->communicator_pt()->nproc() > 1)
849 {
850 MPI_Allreduce(&n,
851 &n2,
852 1,
853 MPI_DOUBLE,
854 MPI_SUM,
855 this->distribution_pt()->communicator_pt()->mpi_comm());
856 }
857 n = n2;
858#endif
859
860 // and return;
861 return n;
862 }
863
864 //============================================================================
865 /// compute the 2 norm of this vector
866 //============================================================================
867 double DoubleVector::norm() const
868 {
869#ifdef PARANOID
870 // paranoid check that the vector is setup
871 if (!this->built())
872 {
873 std::ostringstream error_message;
874 error_message << "This vector must be setup.";
875 throw OomphLibError(
876 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
877 }
878#endif
879
880 // compute the local norm
881 unsigned nrow_local = this->nrow_local();
882 double n = 0;
883 for (unsigned i = 0; i < nrow_local; i++)
884 {
885 n += Values_pt[i] * Values_pt[i];
886 }
887
888 // if this vector is distributed and on multiple processors then gather
889#ifdef OOMPH_HAS_MPI
890 double n2 = n;
891 if (this->distributed() &&
892 this->distribution_pt()->communicator_pt()->nproc() > 1)
893 {
894 MPI_Allreduce(&n,
895 &n2,
896 1,
897 MPI_DOUBLE,
898 MPI_SUM,
899 this->distribution_pt()->communicator_pt()->mpi_comm());
900 }
901 n = n2;
902#endif
903
904 // sqrt the norm
905 n = sqrt(n);
906
907 // and return
908 return n;
909 }
910
911 //============================================================================
912 /// compute the A-norm using the matrix at matrix_pt
913 //============================================================================
914 double DoubleVector::norm(const CRDoubleMatrix* matrix_pt) const
915 {
916#ifdef PARANOID
917 // paranoid check that the vector is setup
918 if (!this->built())
919 {
920 std::ostringstream error_message;
921 error_message << "This vector must be setup.";
922 throw OomphLibError(
923 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
924 }
925 if (!matrix_pt->built())
926 {
927 std::ostringstream error_message;
928 error_message << "The input matrix be built.";
929 throw OomphLibError(
930 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
931 }
932 if (*this->distribution_pt() != *matrix_pt->distribution_pt())
933 {
934 std::ostringstream error_message;
935 error_message << "The distribution of this vector and the matrix at "
936 << "matrix_pt must be the same";
937 throw OomphLibError(
938 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
939 }
940#endif
941
942 // compute the matrix norm
943 DoubleVector x(this->distribution_pt(), 0.0);
944 matrix_pt->multiply(*this, x);
945 return sqrt(this->dot(x));
946 }
947
948 /// output operator
949 std::ostream& operator<<(std::ostream& out, const DoubleVector& v)
950 {
951 // Do the first value outside the loop to get the ", "s right.
952 out << "[" << v[0];
953
954 for (unsigned i = 1, ni = v.nrow_local(); i < ni; i++)
955 {
956 out << ", " << v[i];
957 }
958 out << "]";
959
960 return out;
961 }
962
963 //=================================================================
964 /// Namespace for helper functions for DoubleVectors
965 //=================================================================
966 namespace DoubleVectorHelpers
967 {
968 //===========================================================================
969 /// Concatenate DoubleVectors.
970 /// Takes a Vector of DoubleVectors. If the out vector is built, we will not
971 /// build a new distribution. Otherwise we build a uniform distribution.
972 ///
973 /// The rows of the out vector is seen "as it is" in the in vectors.
974 /// For example, if we have DoubleVectors with distributions A and B,
975 /// distributed across two processors (p0 and p1),
976 ///
977 /// A: [a0] (on p0) B: [b0] (on p0)
978 /// [a1] (on p1) [b1] (on P1),
979 ///
980 /// then the out_vector is
981 ///
982 /// [a0 (on p0)
983 /// a1] (on p0)
984 /// [b0] (on p1)
985 /// b1] (on p1),
986 ///
987 /// Communication is required between processors. The sum of the global
988 /// number of rows in the in vectors must equal to the global number of rows
989 /// in the out vector. This condition must be met if one is to supply an out
990 /// vector with a distribution, otherwise we can let the function generate
991 /// the out vector distribution itself.
992 //===========================================================================
993 void concatenate(const Vector<DoubleVector*>& in_vector_pt,
994 DoubleVector& out_vector)
995 {
996 // How many in vectors to concatenate?
997 unsigned nvectors = in_vector_pt.size();
998
999 // PARANIOD checks which involves the in vectors only
1000#ifdef PARANOID
1001 // Check that there is at least one vector.
1002 if (nvectors == 0)
1003 {
1004 std::ostringstream error_message;
1005 error_message << "There is no vector to concatenate...\n"
1006 << "Perhaps you forgot to fill in_vector_pt?\n";
1007 throw OomphLibError(error_message.str(),
1008 OOMPH_CURRENT_FUNCTION,
1009 OOMPH_EXCEPTION_LOCATION);
1010 }
1011
1012 // Does this vector need concatenating?
1013 if (nvectors == 1)
1014 {
1015 std::ostringstream warning_message;
1016 warning_message << "There is only one vector to concatenate...\n"
1017 << "This does not require concatenating...\n";
1018 OomphLibWarning(warning_message.str(),
1019 OOMPH_CURRENT_FUNCTION,
1020 OOMPH_EXCEPTION_LOCATION);
1021 }
1022
1023 // Check that all the DoubleVectors in in_vector_pt are built
1024 for (unsigned vec_i = 0; vec_i < nvectors; vec_i++)
1025 {
1026 if (!in_vector_pt[vec_i]->built())
1027 {
1028 std::ostringstream error_message;
1029 error_message << "The vector in position " << vec_i
1030 << " is not built.\n"
1031 << "I cannot concatenate an unbuilt vector.\n";
1032 throw OomphLibError(error_message.str(),
1033 OOMPH_CURRENT_FUNCTION,
1034 OOMPH_EXCEPTION_LOCATION);
1035 }
1036 }
1037#endif
1038
1039 // The communicator pointer for the first in vector.
1040 const OomphCommunicator* const comm_pt =
1041 in_vector_pt[0]->distribution_pt()->communicator_pt();
1042
1043 // Check if the first in vector is distributed.
1044 bool distributed = in_vector_pt[0]->distributed();
1045
1046 // If the out vector is not built, build it with a uniform distribution.
1047 if (!out_vector.built())
1048 {
1049 // Nrow for the out vector is the sum of the nrow of the in vectors.
1050 unsigned tmp_nrow = 0;
1051 for (unsigned vec_i = 0; vec_i < nvectors; vec_i++)
1052 {
1053 tmp_nrow += in_vector_pt[vec_i]->nrow();
1054 }
1055
1056 // Build the out vector with uniform distribution.
1057 out_vector.build(
1058 LinearAlgebraDistribution(comm_pt, tmp_nrow, distributed), 0.0);
1059 }
1060 else
1061 {
1062#ifdef PARANOID
1063 // Check that the sum of nrow of in vectors match the nrow in the out
1064 // vectors.
1065 unsigned in_nrow = 0;
1066 for (unsigned vec_i = 0; vec_i < nvectors; vec_i++)
1067 {
1068 in_nrow += in_vector_pt[vec_i]->nrow();
1069 }
1070
1071 if (in_nrow != out_vector.nrow())
1072 {
1073 std::ostringstream error_message;
1074 error_message << "The sum of nrow of the in vectors does not match\n"
1075 << "the nrow of the out vector.\n";
1076 throw OomphLibError(error_message.str(),
1077 OOMPH_CURRENT_FUNCTION,
1078 OOMPH_EXCEPTION_LOCATION);
1079 }
1080#endif
1081 }
1082
1083#ifdef PARANOID
1084 // Check that all communicators of the vectors to concatenate are the same
1085 // by comparing all communicators against the out vector.
1086 const OomphCommunicator out_comm =
1087 *(out_vector.distribution_pt()->communicator_pt());
1088
1089 for (unsigned vec_i = 0; vec_i < nvectors; vec_i++)
1090 {
1091 // Get the Communicator for the current vector.
1092 const OomphCommunicator in_comm =
1093 *(in_vector_pt[vec_i]->distribution_pt()->communicator_pt());
1094
1095 if (out_comm != in_comm)
1096 {
1097 std::ostringstream error_message;
1098 error_message << "The vector in position " << vec_i << " has a\n"
1099 << "different communicator from the out vector.\n";
1100 throw OomphLibError(error_message.str(),
1101 OOMPH_CURRENT_FUNCTION,
1102 OOMPH_EXCEPTION_LOCATION);
1103 }
1104 }
1105
1106 // Check that the distributed boolean is the same for all vectors.
1107 if (out_comm.nproc() != 1)
1108 {
1109 const bool out_distributed = out_vector.distributed();
1110 for (unsigned vec_i = 0; vec_i < nvectors; vec_i++)
1111 {
1112 if (out_distributed != in_vector_pt[vec_i]->distributed())
1113 {
1114 std::ostringstream error_message;
1115 error_message << "The vector in position " << vec_i << " has a\n"
1116 << "different distributed boolean from "
1117 << "the out vector.\n";
1118 throw OomphLibError(error_message.str(),
1119 OOMPH_CURRENT_FUNCTION,
1120 OOMPH_EXCEPTION_LOCATION);
1121 }
1122 }
1123 }
1124#endif
1125
1126
1127 // Now we do the concatenation.
1128 if ((comm_pt->nproc() == 1) || !distributed)
1129 {
1130 // Serial version of the code.
1131 // This is trivial, we simply loop through the in vectors and
1132 // fill in the out vector.
1133
1134 // Out vector index.
1135 unsigned out_i = 0;
1136
1137 // Out vector values.
1138 double* out_value_pt = out_vector.values_pt();
1139
1140 // Loop through the in vectors.
1141 for (unsigned vec_i = 0; vec_i < nvectors; vec_i++)
1142 {
1143 // Nrow of current in vector.
1144 unsigned in_nrow = in_vector_pt[vec_i]->nrow();
1145
1146 // In vector values.
1147 double* in_value_pt = in_vector_pt[vec_i]->values_pt();
1148
1149 // Loop through the entries of this in vector.
1150 for (unsigned i = 0; i < in_nrow; i++)
1151 {
1152 out_value_pt[out_i++] = in_value_pt[i];
1153 }
1154 }
1155 }
1156 // Otherwise we are dealing with a distributed vector.
1157 else
1158 {
1159#ifdef OOMPH_HAS_MPI
1160 // Get the number of processors
1161 unsigned nproc = comm_pt->nproc();
1162
1163 // My rank
1164 unsigned my_rank = comm_pt->my_rank();
1165
1166 // Storage for the data (per processor) to send
1167 Vector<Vector<double>> values_to_send(nproc);
1168
1169 // The sum of the nrow for the in vectors (so far). This is used as an
1170 // offset to calculate the global equation number in the out vector
1171 unsigned long sum_of_vec_nrow = 0;
1172
1173 // Loop over the in vectors and work out:
1174 // out_p: the rank the of receiving processor
1175 // out_local_eqn: the local equation number of the receiving processor
1176 //
1177 // Then put the value and out_local_eqn at out_p in values_to_send
1178
1179 LinearAlgebraDistribution* out_distribution_pt =
1180 out_vector.distribution_pt();
1181 for (unsigned in_vec_i = 0; in_vec_i < nvectors; in_vec_i++)
1182 {
1183 // Loop through the local equations
1184 unsigned in_vec_nrow_local = in_vector_pt[in_vec_i]->nrow_local();
1185 unsigned in_vec_first_row = in_vector_pt[in_vec_i]->first_row();
1186
1187 for (unsigned in_row_i = 0; in_row_i < in_vec_nrow_local; in_row_i++)
1188 {
1189 // Calculate the global equation number for this in_row_i
1190 unsigned out_global_eqn =
1191 in_row_i + in_vec_first_row + sum_of_vec_nrow;
1192
1193 // Get the processor that this global row belongs to.
1194 // The rank_of_global_row(...) function loops through all the
1195 // processors and does two unsigned comparisons. Since we have to do
1196 // this for every row, it may be better to store a list mapping for
1197 // very large number of processors.
1198 unsigned out_p =
1199 out_distribution_pt->rank_of_global_row(out_global_eqn);
1200 // unsigned out_p = out_distribution_pt
1201 // ->rank_of_global_row_map(out_global_eqn);
1202
1203 // Knowing out_p enables us to work out the out_first_row and
1204 // out_local_eqn.
1205 unsigned out_first_row = out_distribution_pt->first_row(out_p);
1206 unsigned out_local_eqn = out_global_eqn - out_first_row;
1207
1208 // Now push back the out_local_eqn and the value
1209 values_to_send[out_p].push_back(out_local_eqn);
1210 values_to_send[out_p].push_back(
1211 (*in_vector_pt[in_vec_i])[in_row_i]);
1212 }
1213
1214 // Update the offset.
1215 sum_of_vec_nrow += in_vector_pt[in_vec_i]->nrow();
1216 }
1217
1218 // Prepare to send the data!
1219
1220 // Storage for the number of data to be sent to each processor.
1221 Vector<int> send_n(nproc, 0);
1222
1223 // Storage for all the values to be send to each processor.
1224 Vector<double> send_values_data;
1225
1226 // Storage location within send_values_data
1227 Vector<int> send_displacement(nproc, 0);
1228
1229 // Get the total amount of data which needs to be sent, so we can
1230 // reserve space for it.
1231 unsigned total_ndata = 0;
1232 for (unsigned rank = 0; rank < nproc; rank++)
1233 {
1234 if (rank != my_rank)
1235 {
1236 total_ndata += values_to_send[rank].size();
1237 }
1238 }
1239
1240 // Now we don't have to re-allocate data/memory when push_back is
1241 // called. Nb. Using push_back without reserving memory may cause
1242 // multiple re-allocation behind the scenes, this is expensive.
1243 send_values_data.reserve(total_ndata);
1244
1245 // Loop over all the processors to "flat pack" the data for sending.
1246 for (unsigned rank = 0; rank < nproc; rank++)
1247 {
1248 // Set the offset for the current processor
1249 send_displacement[rank] = send_values_data.size();
1250
1251 // Don't bother to do anything if
1252 // the processor in the loop is the current processor.
1253 if (rank != my_rank)
1254 {
1255 // Put the values into the send data vector.
1256 unsigned n_data = values_to_send[rank].size();
1257 for (unsigned j = 0; j < n_data; j++)
1258 {
1259 send_values_data.push_back(values_to_send[rank][j]);
1260 } // Loop over the data
1261 } // if rank != my_rank
1262
1263 // Find the number of data to be added to the vector.
1264 send_n[rank] = send_values_data.size() - send_displacement[rank];
1265 } // Loop over processors
1266
1267 // Storage for the number of data to be received from each processor.
1268 Vector<int> receive_n(nproc, 0);
1269 MPI_Alltoall(&send_n[0],
1270 1,
1271 MPI_INT,
1272 &receive_n[0],
1273 1,
1274 MPI_INT,
1275 comm_pt->mpi_comm());
1276
1277 // Prepare the data to be received
1278 // by working out the displacement from the received data.
1279 Vector<int> receive_displacement(nproc, 0);
1280 int receive_data_count = 0;
1281 for (unsigned rank = 0; rank < nproc; rank++)
1282 {
1283 receive_displacement[rank] = receive_data_count;
1284 receive_data_count += receive_n[rank];
1285 }
1286
1287 // Now resize the receive buffer for all data from all processors.
1288 // Make sure that it has size of at least one.
1289 if (receive_data_count == 0)
1290 {
1291 receive_data_count++;
1292 }
1293 Vector<double> receive_values_data(receive_data_count);
1294
1295 // Make sure that the send buffer has size at least one
1296 // so that we don't get a segmentation fault.
1297 if (send_values_data.size() == 0)
1298 {
1299 send_values_data.resize(1);
1300 }
1301
1302 // Now send the data between all processors
1303 MPI_Alltoallv(&send_values_data[0],
1304 &send_n[0],
1305 &send_displacement[0],
1306 MPI_DOUBLE,
1307 &receive_values_data[0],
1308 &receive_n[0],
1309 &receive_displacement[0],
1310 MPI_DOUBLE,
1311 comm_pt->mpi_comm());
1312
1313 // Data from all other processors are stored in:
1314 // receive_values_data
1315 // Data already on this processor is stored in:
1316 // values_to_send[my_rank]
1317
1318 // Loop through the data on this processor.
1319 unsigned location_i = 0;
1320 unsigned my_values_to_send_size = values_to_send[my_rank].size();
1321 while (location_i < my_values_to_send_size)
1322 {
1323 out_vector[unsigned(values_to_send[my_rank][location_i])] =
1324 values_to_send[my_rank][location_i + 1];
1325
1326 location_i += 2;
1327 }
1328
1329 // Before we loop through the data on other processors, we need to check
1330 // if any data has been received.
1331 bool data_has_been_received = false;
1332 unsigned send_rank = 0;
1333 while (send_rank < nproc)
1334 {
1335 if (receive_n[send_rank] > 0)
1336 {
1337 data_has_been_received = true;
1338 break;
1339 }
1340 send_rank++;
1341 }
1342
1343 location_i = 0;
1344 if (data_has_been_received)
1345 {
1346 unsigned receive_values_data_size = receive_values_data.size();
1347 while (location_i < receive_values_data_size)
1348 {
1349 out_vector[unsigned(receive_values_data[location_i])] =
1350 receive_values_data[location_i + 1];
1351 location_i += 2;
1352 }
1353 }
1354#else
1355 {
1356 std::ostringstream error_message;
1357 error_message << "I don't know what to do with distributed vectors\n"
1358 << "without MPI... :(";
1359 throw OomphLibError(error_message.str(),
1360 OOMPH_CURRENT_FUNCTION,
1361 OOMPH_EXCEPTION_LOCATION);
1362 }
1363#endif
1364 }
1365 } // function concatenate
1366
1367 //===========================================================================
1368 /// Wrapper around the other concatenate(...) function.
1369 /// Be careful with Vector of vectors. If the DoubleVectors are resized,
1370 /// there could be reallocation of memory. If we wanted to use the function
1371 /// which takes a Vector of pointers to DoubleVectors, we would either have
1372 /// to invoke new and remember to delete, or create a temporary Vector to
1373 /// store pointers to the DoubleVector objects.
1374 /// This wrapper is meant to make life easier for the user by avoiding calls
1375 /// to new/delete AND without creating a temporary vector of pointers to
1376 /// DoubleVectors.
1377 /// If we had C++ 11, this would be so much nicer since we can use smart
1378 /// pointers which will delete themselves, so we do not have to remember
1379 /// to delete!
1380 //===========================================================================
1381 void concatenate(Vector<DoubleVector>& in_vector, DoubleVector& out_vector)
1382 {
1383 const unsigned n_in_vector = in_vector.size();
1384
1385 Vector<DoubleVector*> in_vector_pt(n_in_vector, 0);
1386
1387 for (unsigned i = 0; i < n_in_vector; i++)
1388 {
1389 in_vector_pt[i] = &in_vector[i];
1390 }
1391
1392 DoubleVectorHelpers::concatenate(in_vector_pt, out_vector);
1393 } // function concatenate
1394
1395 //===========================================================================
1396 /// Split a DoubleVector into the out DoubleVectors.
1397 /// Let vec_A be the in Vector, and let vec_B and vec_C be the out vectors.
1398 /// Then the splitting of vec_A is depicted below:
1399 /// vec_A: [a0 (on p0)
1400 /// a1] (on p0)
1401 /// [a2 (on p1)
1402 /// a3] (on p1)
1403 ///
1404 /// vec_B: [a0] (on p0) vec_C: [a2] (on p0)
1405 /// [a1] (on p1) [a3] (on p1)
1406 ///
1407 /// Communication is required between processors.
1408 /// The out_vector_pt must contain pointers to DoubleVector which has
1409 /// already been built with the correct distribution; the sum of the number
1410 /// of global row of the out vectors must be the same the number of global
1411 /// rows of the in vector.
1412 //===========================================================================
1413 void split(const DoubleVector& in_vector,
1414 Vector<DoubleVector*>& out_vector_pt)
1415 {
1416 // How many out vectors do we have?
1417 unsigned nvec = out_vector_pt.size();
1418#ifdef PARANOID
1419
1420 // Check that the in vector is built.
1421 if (!in_vector.built())
1422 {
1423 std::ostringstream error_message;
1424 error_message << "The in_vector is not built.\n"
1425 << "Please build it!.\n";
1426 throw OomphLibError(error_message.str(),
1427 OOMPH_CURRENT_FUNCTION,
1428 OOMPH_EXCEPTION_LOCATION);
1429 }
1430
1431 // Check that all the out vectors are built.
1432 for (unsigned vec_i = 0; vec_i < nvec; vec_i++)
1433 {
1434 if (!out_vector_pt[vec_i]->built())
1435 {
1436 std::ostringstream error_message;
1437 error_message << "The vector at position " << vec_i
1438 << " is not built.\n"
1439 << "Please build it!.\n";
1440 throw OomphLibError(error_message.str(),
1441 OOMPH_CURRENT_FUNCTION,
1442 OOMPH_EXCEPTION_LOCATION);
1443 }
1444 }
1445
1446 // Check that the sum of the nrow from out vectors is the same as the
1447 // nrow from in_vector.
1448 unsigned out_nrow_sum = 0;
1449 for (unsigned vec_i = 0; vec_i < nvec; vec_i++)
1450 {
1451 out_nrow_sum += out_vector_pt[vec_i]->nrow();
1452 }
1453
1454 if (in_vector.nrow() != out_nrow_sum)
1455 {
1456 std::ostringstream error_message;
1457 error_message << "The global number of rows in the in_vector\n"
1458 << "is not equal to the sum of the global nrows\n"
1459 << "of the in vectors.\n";
1460 throw OomphLibError(error_message.str(),
1461 OOMPH_CURRENT_FUNCTION,
1462 OOMPH_EXCEPTION_LOCATION);
1463 }
1464
1465 // Check that all communicators are the same. We use a communicator to
1466 // get the number of processors and my_rank. So we would like them to be
1467 // the same for in_vector and all out vectors.
1468 const OomphCommunicator in_vector_comm =
1469 *(in_vector.distribution_pt()->communicator_pt());
1470 for (unsigned vec_i = 0; vec_i < nvec; vec_i++)
1471 {
1472 const OomphCommunicator dist_i_comm =
1473 *(out_vector_pt[vec_i]->distribution_pt()->communicator_pt());
1474
1475 if (in_vector_comm != dist_i_comm)
1476 {
1477 std::ostringstream error_message;
1478 error_message << "The communicator for the distribution in the \n"
1479 << "position " << vec_i
1480 << " is not the same as the in_vector\n";
1481 throw OomphLibError(error_message.str(),
1482 OOMPH_CURRENT_FUNCTION,
1483 OOMPH_EXCEPTION_LOCATION);
1484 }
1485 }
1486
1487 // Check that the distributed boolean is the same for all vectors.
1488 bool para_distributed = in_vector.distributed();
1489
1490 for (unsigned vec_i = 0; vec_i < nvec; vec_i++)
1491 {
1492 if (para_distributed != out_vector_pt[vec_i]->distributed())
1493 {
1494 std::ostringstream error_message;
1495 error_message
1496 << "The vector in position " << vec_i << " does not \n"
1497 << " have the same distributed boolean as the in_vector\n";
1498 throw OomphLibError(error_message.str(),
1499 OOMPH_CURRENT_FUNCTION,
1500 OOMPH_EXCEPTION_LOCATION);
1501 }
1502 }
1503
1504#endif
1505
1506 // The communicator.
1507 const OomphCommunicator* const comm_pt =
1508 in_vector.distribution_pt()->communicator_pt();
1509
1510 // Is this distributed?
1511 bool distributed = in_vector.distributed();
1512
1513 // The serial code.
1514 if ((comm_pt->nproc() == 1) || !distributed)
1515 {
1516 // Serial version of the code: loop through all the out vectors and
1517 // insert the elements of in_vector.
1518
1519 // index for in vector, and in vector values.
1520 unsigned in_vec_i = 0;
1521 double* in_value_pt = in_vector.values_pt();
1522
1523 // Fill in the out vectors.
1524 for (unsigned out_vec_i = 0; out_vec_i < nvec; out_vec_i++)
1525 {
1526 // out vector nrow and values.
1527 unsigned out_nrow = out_vector_pt[out_vec_i]->nrow();
1528 double* out_value_pt = out_vector_pt[out_vec_i]->values_pt();
1529
1530 // Fill in the current out vector.
1531 for (unsigned out_val_i = 0; out_val_i < out_nrow; out_val_i++)
1532 {
1533 out_value_pt[out_val_i] = in_value_pt[in_vec_i++];
1534 }
1535 }
1536 }
1537 // Otherwise we are dealing with a distributed vector.
1538 else
1539 {
1540#ifdef OOMPH_HAS_MPI
1541 // For each entry in the in_vector, we need to work out:
1542 // 1) Which out vector this entry belongs to,
1543 // 2) which processor to send the data to and
1544 // 3) the local equation number in the out vector.
1545 //
1546 // We know the in_local_eqn, we can work out the in_global_eqn.
1547 //
1548 // From in_global_eqn we can work out the out vector and
1549 // the out_global_eqn.
1550 //
1551 // The out_global_eqn allows us to determine which processor to send to.
1552 // With the out_p (processor to send data to) and out vector, we get the
1553 // out_first_row which then allows us to work out the out_local_eqn.
1554
1555
1556 // Get the number of processors
1557 unsigned nproc = comm_pt->nproc();
1558
1559 // My rank
1560 unsigned my_rank = comm_pt->my_rank();
1561
1562 // Storage for the data (per processor) to send.
1563 Vector<Vector<double>> values_to_send(nproc);
1564
1565 // Sum of the nrow of the out vectors so far. This is used to work out
1566 // which out_vector a in_global_eqn belongs to.
1567 Vector<unsigned> sum_of_out_nrow(nvec + 1);
1568 for (unsigned vec_i = 0; vec_i < nvec; vec_i++)
1569 {
1570 sum_of_out_nrow[vec_i + 1] =
1571 sum_of_out_nrow[vec_i] + out_vector_pt[vec_i]->nrow();
1572 }
1573
1574 // Loop through the in_vector local values.
1575 unsigned in_nrow_local = in_vector.nrow_local();
1576 for (unsigned in_local_eqn = 0; in_local_eqn < in_nrow_local;
1577 in_local_eqn++)
1578 {
1579 // The global equation number of this row.
1580 unsigned in_global_eqn = in_local_eqn + in_vector.first_row();
1581
1582 // Which out_vector does this in_global_eqn belong to?
1583 unsigned out_vector_i = 0;
1584 while (in_global_eqn < sum_of_out_nrow[out_vector_i] ||
1585 in_global_eqn >= sum_of_out_nrow[out_vector_i + 1])
1586 {
1587 out_vector_i++;
1588 }
1589
1590 // The out_global_eqn
1591 // (this is the global equation in the current out vector)
1592 unsigned out_global_eqn =
1593 in_global_eqn - sum_of_out_nrow[out_vector_i];
1594
1595 // The processor to send this row to.
1596 unsigned out_p =
1597 out_vector_pt[out_vector_i]->distribution_pt()->rank_of_global_row(
1598 out_global_eqn);
1599
1600 // The local_eqn in the out_vector_i
1601 unsigned out_local_eqn =
1602 out_global_eqn -
1603 out_vector_pt[out_vector_i]->distribution_pt()->first_row(out_p);
1604
1605
1606 // Fill in the data to send
1607
1608 // Which out vector to put this data in.
1609 values_to_send[out_p].push_back(out_vector_i);
1610
1611 // The local equation of the data.
1612 values_to_send[out_p].push_back(out_local_eqn);
1613
1614 // The actual data.
1615 values_to_send[out_p].push_back(in_vector[in_local_eqn]);
1616 }
1617
1618 // Prepare to send the data!
1619
1620 // Storage for the number of data to be sent to each processor.
1621 Vector<int> send_n(nproc, 0);
1622
1623 // Storage for all the values to be send to each processor.
1624 Vector<double> send_values_data;
1625
1626 // Storage location within send_values_data
1627 Vector<int> send_displacement(nproc, 0);
1628
1629 // Get the total amount of data which needs to be sent, so we can
1630 // reserve space for it.
1631 unsigned total_ndata = 0;
1632 for (unsigned rank = 0; rank < nproc; rank++)
1633 {
1634 if (rank != my_rank)
1635 {
1636 total_ndata += values_to_send[rank].size();
1637 }
1638 }
1639
1640 // Now we don't have to re-allocate data/memory when push_back is
1641 // called. Nb. Using push_back without reserving memory may cause
1642 // multiple re-allocation behind the scenes, this is expensive.
1643 send_values_data.reserve(total_ndata);
1644
1645 // Loop over all the processors to "flat pack" the data for sending.
1646 for (unsigned rank = 0; rank < nproc; rank++)
1647 {
1648 // Set the offset for the current processor
1649 send_displacement[rank] = send_values_data.size();
1650
1651 // Don't bother to do anything if
1652 // the processor in the loop is the current processor.
1653 if (rank != my_rank)
1654 {
1655 // Put the values into the send data vector.
1656 unsigned n_data = values_to_send[rank].size();
1657 for (unsigned j = 0; j < n_data; j++)
1658 {
1659 send_values_data.push_back(values_to_send[rank][j]);
1660 } // Loop over the data
1661 } // if rank != my_rank
1662
1663 // Find the number of data to be added to the vector.
1664 send_n[rank] = send_values_data.size() - send_displacement[rank];
1665 } // Loop over processors
1666
1667 // Storage for the number of data to be received from each processor.
1668 Vector<int> receive_n(nproc, 0);
1669 MPI_Alltoall(&send_n[0],
1670 1,
1671 MPI_INT,
1672 &receive_n[0],
1673 1,
1674 MPI_INT,
1675 comm_pt->mpi_comm());
1676
1677 // Prepare the data to be received
1678 // by working out the displacement from the received data.
1679 Vector<int> receive_displacement(nproc, 0);
1680 int receive_data_count = 0;
1681 for (unsigned rank = 0; rank < nproc; rank++)
1682 {
1683 receive_displacement[rank] = receive_data_count;
1684 receive_data_count += receive_n[rank];
1685 }
1686
1687 // Now resize the receive buffer for all data from all processors.
1688 // Make sure that it has size of at least one.
1689 if (receive_data_count == 0)
1690 {
1691 receive_data_count++;
1692 }
1693 Vector<double> receive_values_data(receive_data_count);
1694
1695 // Make sure that the send buffer has size at least one
1696 // so that we don't get a segmentation fault.
1697 if (send_values_data.size() == 0)
1698 {
1699 send_values_data.resize(1);
1700 }
1701
1702 // Now send the data between all processors
1703 MPI_Alltoallv(&send_values_data[0],
1704 &send_n[0],
1705 &send_displacement[0],
1706 MPI_DOUBLE,
1707 &receive_values_data[0],
1708 &receive_n[0],
1709 &receive_displacement[0],
1710 MPI_DOUBLE,
1711 comm_pt->mpi_comm());
1712
1713 // Data from all other processors are stored in:
1714 // receive_values_data
1715 // Data already on this processor is stored in:
1716 // values_to_send[my_rank]
1717 //
1718
1719 // Index for values_to_send Vector.
1720 unsigned location_i = 0;
1721 // Loop through the data on this processor
1722 unsigned my_values_to_send_size = values_to_send[my_rank].size();
1723 while (location_i < my_values_to_send_size)
1724 {
1725 // The vector to put the values in.
1726 unsigned out_vector_i =
1727 unsigned(values_to_send[my_rank][location_i++]);
1728
1729 // Where to put the value.
1730 unsigned out_local_eqn =
1731 unsigned(values_to_send[my_rank][location_i++]);
1732
1733 // The actual value!
1734 double out_value = values_to_send[my_rank][location_i++];
1735
1736 // Insert the value in the out vector.
1737 (*out_vector_pt[out_vector_i])[out_local_eqn] = out_value;
1738 }
1739
1740 // Before we loop through the data on other processors, we need to check
1741 // if any data has been received. This is because the
1742 // receive_values_data has been resized to at least one, even if no data
1743 // is sent.
1744 bool data_has_been_received = false;
1745 unsigned send_rank = 0;
1746 while (send_rank < nproc)
1747 {
1748 if (receive_n[send_rank] > 0)
1749 {
1750 data_has_been_received = true;
1751 break;
1752 }
1753 send_rank++;
1754 }
1755
1756 // Reset the index, it is now being used to index the
1757 // receive_values_data vector.
1758 location_i = 0;
1759 if (data_has_been_received)
1760 {
1761 // Extract the data and put it into the out vector.
1762 unsigned receive_values_data_size = receive_values_data.size();
1763 while (location_i < receive_values_data_size)
1764 {
1765 // Which out vector to put the value in?
1766 unsigned out_vector_i = unsigned(receive_values_data[location_i++]);
1767
1768 // Where in the out vector to put the value?
1769 unsigned out_local_eqn =
1770 unsigned(receive_values_data[location_i++]);
1771
1772 // The value to put in.
1773 double out_value = receive_values_data[location_i++];
1774
1775 // Insert the value in the out vector.
1776 (*out_vector_pt[out_vector_i])[out_local_eqn] = out_value;
1777 }
1778 }
1779#else
1780 {
1781 std::ostringstream error_message;
1782 error_message << "You have a distributed vector but with no mpi...\n"
1783 << "I don't know what to do :( \n";
1784 throw OomphLibError(
1785 error_message.str(), "RYARAYERR", OOMPH_EXCEPTION_LOCATION);
1786 }
1787#endif
1788 }
1789 } // function split(...)
1790
1791 //===========================================================================
1792 /// Wrapper around the other split(...) function.
1793 /// Be careful with Vector of vectors. If the DoubleVectors are resized,
1794 /// there could be reallocation of memory. If we wanted to use the function
1795 /// which takes a Vector of pointers to DoubleVectors, we would either have
1796 /// to invoke new and remember to delete, or create a temporary Vector to
1797 /// store pointers to the DoubleVector objects.
1798 /// This wrapper is meant to make life easier for the user by avoiding calls
1799 /// to new/delete AND without creating a temporary vector of pointers to
1800 /// DoubleVectors.
1801 /// If we had C++ 11, this would be so much nicer since we can use smart
1802 /// pointers which will delete themselves, so we do not have to remember
1803 /// to delete!
1804 //===========================================================================
1805 void split(const DoubleVector& in_vector, Vector<DoubleVector>& out_vector)
1806 {
1807 const unsigned n_out_vector = out_vector.size();
1808 Vector<DoubleVector*> out_vector_pt(n_out_vector, 0);
1809
1810 for (unsigned i = 0; i < n_out_vector; i++)
1811 {
1812 out_vector_pt[i] = &out_vector[i];
1813 }
1814
1815 DoubleVectorHelpers::split(in_vector, out_vector_pt);
1816 } // function split(...)
1817
1818 //===========================================================================
1819 /// Concatenate DoubleVectors.
1820 /// Takes a Vector of DoubleVectors. If the out vector is built, we will not
1821 /// build a new distribution. Otherwise a new distribution will be built
1822 /// using LinearAlgebraDistribution::concatenate(...).
1823 ///
1824 /// The out vector has its rows permuted according to the individual
1825 /// distributions of the in vectors. For example, if we have DoubleVectors
1826 /// with distributions A and B, distributed across two processors
1827 /// (p0 and p1),
1828 ///
1829 /// A: [a0] (on p0) B: [b0] (on p0)
1830 /// [a1] (on p1) [b1] (on P1),
1831 ///
1832 /// then the out_vector is
1833 ///
1834 /// [a0 (on p0)
1835 /// b0] (on p0)
1836 /// [a1 (on p1)
1837 /// b1] (on p1),
1838 ///
1839 /// as opposed to
1840 ///
1841 /// [a0 (on p0)
1842 /// a1] (on p0)
1843 /// [b0 (on p1)
1844 /// b1] (on p1).
1845 ///
1846 /// Note (1): The out vector may not be uniformly distributed even
1847 /// if the in vectors have uniform distributions. The nrow_local of the
1848 /// out vector will be the sum of the nrow_local of the in vectors.
1849 /// Try this out with two distributions of global rows 3 and 5, uniformly
1850 /// distributed across two processors. Compare this against a distribution
1851 /// of global row 8 distributed across two processors.
1852 ///
1853 /// There are no MPI send and receive, the data stays on the processor
1854 /// as defined by the distributions from the in vectors.
1855 //===========================================================================
1857 const Vector<DoubleVector*>& in_vector_pt, DoubleVector& out_vector)
1858 {
1859 // How many in vectors do we want to concatenate?
1860 unsigned nvectors = in_vector_pt.size();
1861
1862 // PARANOID checks which involves the in vectors only.
1863#ifdef PARANOID
1864 // Check that there is at least one vector.
1865 if (nvectors == 0)
1866 {
1867 std::ostringstream error_message;
1868 error_message << "There is no vector to concatenate...\n"
1869 << "Perhaps you forgot to fill in_vector_pt?\n";
1870 throw OomphLibError(error_message.str(),
1871 OOMPH_CURRENT_FUNCTION,
1872 OOMPH_EXCEPTION_LOCATION);
1873 }
1874
1875 // Does this vector need concatenating?
1876 if (nvectors == 1)
1877 {
1878 std::ostringstream warning_message;
1879 warning_message << "There is only one vector to concatenate...\n"
1880 << "This does not require concatenating...\n";
1881 OomphLibWarning(warning_message.str(),
1882 OOMPH_CURRENT_FUNCTION,
1883 OOMPH_EXCEPTION_LOCATION);
1884 }
1885
1886 // Check that all the DoubleVectors in in_vector_pt are built
1887 for (unsigned vec_i = 0; vec_i < nvectors; vec_i++)
1888 {
1889 if (!in_vector_pt[vec_i]->built())
1890 {
1891 std::ostringstream error_message;
1892 error_message << "The vector in position " << vec_i
1893 << " is not built.\n"
1894 << "I cannot concatenate an unbuilt vector.\n";
1895 throw OomphLibError(error_message.str(),
1896 OOMPH_CURRENT_FUNCTION,
1897 OOMPH_EXCEPTION_LOCATION);
1898 }
1899 }
1900#endif
1901
1902 // If the out vector is not built, build it with the correct distribution.
1903 if (!out_vector.built())
1904 {
1905 Vector<LinearAlgebraDistribution*> in_distribution_pt(nvectors, 0);
1906 for (unsigned vec_i = 0; vec_i < nvectors; vec_i++)
1907 {
1908 in_distribution_pt[vec_i] = in_vector_pt[vec_i]->distribution_pt();
1909 }
1910
1911 LinearAlgebraDistribution tmp_distribution;
1913 tmp_distribution);
1914 out_vector.build(tmp_distribution, 0.0);
1915 }
1916
1917 // PARANOID checks which involves all in vectors and out vectors.
1918#ifdef PARANOID
1919
1920 // Check that all communicators of the vectors to concatenate are the same
1921 // by comparing all communicators against the out vector.
1922 const OomphCommunicator out_comm =
1923 *(out_vector.distribution_pt()->communicator_pt());
1924
1925 for (unsigned vec_i = 0; vec_i < nvectors; vec_i++)
1926 {
1927 // Get the Communicator for the current vector.
1928 const OomphCommunicator in_comm =
1929 *(in_vector_pt[vec_i]->distribution_pt()->communicator_pt());
1930
1931 if (out_comm != in_comm)
1932 {
1933 std::ostringstream error_message;
1934 error_message << "The vector in position " << vec_i << " has a\n"
1935 << "different communicator from the out vector.\n";
1936 throw OomphLibError(error_message.str(),
1937 OOMPH_CURRENT_FUNCTION,
1938 OOMPH_EXCEPTION_LOCATION);
1939 }
1940 }
1941
1942 // Check that the distributed boolean is the same for all vectors.
1943 if (out_comm.nproc() > 1)
1944 {
1945 const bool out_distributed = out_vector.distributed();
1946 for (unsigned vec_i = 0; vec_i < nvectors; vec_i++)
1947 {
1948 if (out_distributed != in_vector_pt[vec_i]->distributed())
1949 {
1950 std::ostringstream error_message;
1951 error_message << "The vector in position " << vec_i << " has a\n"
1952 << "different distributed boolean from the "
1953 << "out vector.\n";
1954 throw OomphLibError(error_message.str(),
1955 OOMPH_CURRENT_FUNCTION,
1956 OOMPH_EXCEPTION_LOCATION);
1957 }
1958 }
1959 }
1960
1961 // Check that the distribution from the out vector is indeed the
1962 // same as the one created by
1963 // LinearAlgebraDistributionHelpers::concatenate(...). This test is
1964 // redundant if the out_vector is not built to begin with.
1965
1966 // Create tmp_distribution, a concatenation of all distributions from
1967 // the in vectors.
1968 Vector<LinearAlgebraDistribution*> in_distribution_pt(nvectors, 0);
1969 for (unsigned vec_i = 0; vec_i < nvectors; vec_i++)
1970 {
1971 in_distribution_pt[vec_i] = in_vector_pt[vec_i]->distribution_pt();
1972 }
1973
1974 LinearAlgebraDistribution tmp_distribution;
1976 tmp_distribution);
1977 // The the distribution from the out vector.
1978 LinearAlgebraDistribution out_distribution =
1979 *(out_vector.distribution_pt());
1980
1981 // Compare them!
1982 if (tmp_distribution != out_distribution)
1983 {
1984 std::ostringstream error_message;
1985 error_message << "The distribution of the out vector is not correct.\n"
1986 << "Please call the function with a cleared out vector,\n"
1987 << "or compare the distribution of the out vector with\n"
1988 << "the distribution created by\n"
1989 << "LinearAlgebraDistributionHelpers::concatenate(...)\n";
1990 throw OomphLibError(error_message.str(),
1991 OOMPH_CURRENT_FUNCTION,
1992 OOMPH_EXCEPTION_LOCATION);
1993 }
1994
1995 // Do not need these distributions.
1996 tmp_distribution.clear();
1997 out_distribution.clear();
1998#endif
1999
2000
2001 unsigned out_value_offset = 0;
2002
2003 double* out_value_pt = out_vector.values_pt();
2004
2005 // Loop through the vectors.
2006 for (unsigned vec_i = 0; vec_i < nvectors; vec_i++)
2007 {
2008 // Get the nrow_local and
2009 // pointer to the values for the current in vector.
2010 unsigned in_vector_nrow_local = in_vector_pt[vec_i]->nrow_local();
2011 double* in_vector_value_pt = in_vector_pt[vec_i]->values_pt();
2012
2013 // Loop through the local values and inset them into the out_vector.
2014 for (unsigned val_i = 0; val_i < in_vector_nrow_local; val_i++)
2015 {
2016 out_value_pt[out_value_offset + val_i] = in_vector_value_pt[val_i];
2017 }
2018
2019 // Update the offset.
2020 out_value_offset += in_vector_nrow_local;
2021 }
2022 } // function concatenate_without_communication
2023
2024 //===========================================================================
2025 /// Wrapper around the other concatenate_without_communication(...)
2026 /// function.
2027 /// Be careful with Vector of vectors. If the DoubleVectors are resized,
2028 /// there could be reallocation of memory. If we wanted to use the function
2029 /// which takes a Vector of pointers to DoubleVectors, we would either have
2030 /// to invoke new and remember to delete, or create a temporary Vector to
2031 /// store pointers to the DoubleVector objects.
2032 /// This wrapper is meant to make life easier for the user by avoiding calls
2033 /// to new/delete AND without creating a temporary vector of pointers to
2034 /// DoubleVectors.
2035 /// If we had C++ 11, this would be so much nicer since we can use smart
2036 /// pointers which will delete themselves, so we do not have to remember
2037 /// to delete!
2038 //===========================================================================
2040 DoubleVector& out_vector)
2041 {
2042 const unsigned n_in_vector = in_vector.size();
2043
2044 Vector<DoubleVector*> in_vector_pt(n_in_vector, 0);
2045
2046 for (unsigned i = 0; i < n_in_vector; i++)
2047 {
2048 in_vector_pt[i] = &in_vector[i];
2049 }
2050
2052 out_vector);
2053 } // function concatenate_without_communication
2054
2055 //===========================================================================
2056 /// Split a DoubleVector into the out DoubleVectors.
2057 /// Data stays on its current processor, no data is sent between processors.
2058 /// This results in our vectors which are a permutation of the in vector.
2059 ///
2060 /// Let vec_A be the in Vector, and let vec_B and vec_C be the out vectors.
2061 /// Then the splitting of vec_A is depicted below:
2062 /// vec_A: [a0 (on p0)
2063 /// a1] (on p0)
2064 /// [a2 (on p1)
2065 /// a3] (on p1)
2066 ///
2067 /// vec_B: [a0] (on p0) vec_C: [a1] (on p0)
2068 /// [a2] (on p1) [a3] (on p1).
2069 ///
2070 /// This means that the distribution of the in vector MUST be a
2071 /// concatenation of the out vector distributions, refer to
2072 /// LinearAlgebraDistributionHelpers::concatenate(...) to concatenate
2073 /// distributions.
2074 //===========================================================================
2076 Vector<DoubleVector*>& out_vector_pt)
2077 {
2078 // How many out vectors do we need?
2079 unsigned nvec = out_vector_pt.size();
2080
2081#ifdef PARANOID
2082 // Check that in_vector is built
2083 if (!in_vector.built())
2084 {
2085 std::ostringstream error_message;
2086 error_message << "The in_vector is not built.\n"
2087 << "Please build it!.\n";
2088 throw OomphLibError(error_message.str(),
2089 OOMPH_CURRENT_FUNCTION,
2090 OOMPH_EXCEPTION_LOCATION);
2091 }
2092
2093 // Check that all out vectors are built.
2094 for (unsigned vec_i = 0; vec_i < nvec; vec_i++)
2095 {
2096 if (!out_vector_pt[vec_i]->built())
2097 {
2098 std::ostringstream error_message;
2099 error_message << "The vector at position " << vec_i
2100 << " is not built.\n"
2101 << "Please build it!.\n";
2102 throw OomphLibError(error_message.str(),
2103 OOMPH_CURRENT_FUNCTION,
2104 OOMPH_EXCEPTION_LOCATION);
2105 }
2106 }
2107
2108 // Check that the concatenation of distributions from the out vectors is
2109 // the same as the distribution from in_vector.
2110
2111 // Create the distribution from out_distribution.
2112 Vector<LinearAlgebraDistribution*> tmp_out_distribution_pt(nvec, 0);
2113 for (unsigned vec_i = 0; vec_i < nvec; vec_i++)
2114 {
2115 tmp_out_distribution_pt[vec_i] =
2116 out_vector_pt[vec_i]->distribution_pt();
2117 }
2118
2119 LinearAlgebraDistribution tmp_distribution;
2120 LinearAlgebraDistributionHelpers::concatenate(tmp_out_distribution_pt,
2121 tmp_distribution);
2122 // Compare the distributions
2123 if (tmp_distribution != *(in_vector.distribution_pt()))
2124 {
2125 std::ostringstream error_message;
2126 error_message << "The distribution from the in vector is incorrect.\n"
2127 << "It must be a concatenation of all the distributions\n"
2128 << "from the out vectors.\n";
2129 throw OomphLibError(error_message.str(),
2130 OOMPH_CURRENT_FUNCTION,
2131 OOMPH_EXCEPTION_LOCATION);
2132 }
2133
2134 // Clear the distribution.
2135 tmp_distribution.clear();
2136
2137 // Check that all communicators are the same. We use a communicator to
2138 // get the number of processors and my_rank. So we would like them to be
2139 // the same for the in vector and all the out vectors.
2140 const OomphCommunicator in_vector_comm =
2141 *(in_vector.distribution_pt()->communicator_pt());
2142 for (unsigned vec_i = 0; vec_i < nvec; vec_i++)
2143 {
2144 const OomphCommunicator vec_i_comm =
2145 *(out_vector_pt[vec_i]->distribution_pt()->communicator_pt());
2146
2147 if (in_vector_comm != vec_i_comm)
2148 {
2149 std::ostringstream error_message;
2150 error_message << "The communicator for the vector in position\n"
2151 << vec_i << " is not the same as the in_vector\n"
2152 << "communicator.";
2153 throw OomphLibError(error_message.str(),
2154 OOMPH_CURRENT_FUNCTION,
2155 OOMPH_EXCEPTION_LOCATION);
2156 }
2157 }
2158
2159 // Check that if the in vector is distributed, then all the out vectors
2160 // are also distributed.
2161 if (in_vector_comm.nproc() > 1)
2162 {
2163 bool in_distributed = in_vector.distributed();
2164 for (unsigned vec_i = 0; vec_i < nvec; vec_i++)
2165 {
2166 if (in_distributed != out_vector_pt[vec_i]->distributed())
2167 {
2168 std::ostringstream error_message;
2169 error_message << "The vector in position " << vec_i
2170 << " does not have\n"
2171 << "the same distributed boolean as the in vector";
2172 throw OomphLibError(error_message.str(),
2173 OOMPH_CURRENT_FUNCTION,
2174 OOMPH_EXCEPTION_LOCATION);
2175 }
2176 }
2177 }
2178#endif
2179
2180 // Loop through the sub vectors and insert the values from the
2181 // in vector.
2182 double* in_value_pt = in_vector.values_pt();
2183 unsigned in_value_offset = 0;
2184 for (unsigned vec_i = 0; vec_i < nvec; vec_i++)
2185 {
2186 // The nrow_local and values for the current out vector.
2187 unsigned out_nrow_local = out_vector_pt[vec_i]->nrow_local();
2188 double* out_value_pt = out_vector_pt[vec_i]->values_pt();
2189
2190 // Loop through the local values of out vector.
2191 for (unsigned val_i = 0; val_i < out_nrow_local; val_i++)
2192 {
2193 out_value_pt[val_i] = in_value_pt[in_value_offset + val_i];
2194 }
2195
2196 // Update the offset.
2197 in_value_offset += out_nrow_local;
2198 }
2199 } // function split_distribution_vector
2200
2201 //===========================================================================
2202 /// Wrapper around the other split_without_communication(...)
2203 /// function.
2204 /// Be careful with Vector of vectors. If the DoubleVectors are resized,
2205 /// there could be reallocation of memory. If we wanted to use the function
2206 /// which takes a Vector of pointers to DoubleVectors, we would either have
2207 /// to invoke new and remember to delete, or create a temporary Vector to
2208 /// store pointers to the DoubleVector objects.
2209 /// This wrapper is meant to make life easier for the user by avoiding calls
2210 /// to new/delete AND without creating a temporary vector of pointers to
2211 /// DoubleVectors.
2212 /// If we had C++ 11, this would be so much nicer since we can use smart
2213 /// pointers which will delete themselves, so we do not have to remember
2214 /// to delete!
2215 //===========================================================================
2217 Vector<DoubleVector>& out_vector)
2218 {
2219 const unsigned n_out_vector = out_vector.size();
2220
2221 Vector<DoubleVector*> out_vector_pt(n_out_vector, 0);
2222
2223 for (unsigned i = 0; i < n_out_vector; i++)
2224 {
2225 out_vector_pt[i] = &out_vector[i];
2226 }
2227
2229 out_vector_pt);
2230
2231 } // function split_distribution_vector
2232
2233 } // namespace DoubleVectorHelpers
2234
2235} // namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:888
void multiply(const DoubleVector &x, DoubleVector &soln) const
Multiply the matrix by the vector x: soln=Ax.
Definition: matrices.cc:1782
bool built() const
access function to the Built flag - indicates whether the matrix has been build - i....
Definition: matrices.h:1210
bool distributed() const
distribution is serial or distributed
unsigned nrow() const
access function to the number of global rows.
bool distribution_built() const
if the communicator_pt is null then the distribution is not setup then false is returned,...
unsigned nrow_local() const
access function for the num of local rows on this processor.
unsigned first_row() const
access function for the first row on this processor
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
A vector in the mathematical sense, initially developed for linear algebra type applications....
Definition: double_vector.h:58
void initialise(const double &v)
initialise the whole vector with value v
void operator+=(const DoubleVector &v)
+= operator with another vector
void output_local_values(std::ostream &outfile, const int &output_precision=-1) const
output the local contents of the vector
double max() const
returns the maximum coefficient
bool operator==(const DoubleVector &v)
== operator
void operator*=(const double &d)
multiply by a double
void build(const DoubleVector &old_vector)
Just copys the argument DoubleVector.
void redistribute(const LinearAlgebraDistribution *const &dist_pt)
The contents of the vector are redistributed to match the new distribution. In a non-MPI rebuild this...
double norm() const
compute the 2 norm of this vector
bool Internal_values
Boolean flag to indicate whether the vector's data (values_pt) is owned by this vector.
void output(std::ostream &outfile, const int &output_precision=-1) const
output the global contents of the vector
bool built() const
void operator-=(const DoubleVector &v)
-= operator with another vector
void operator/=(const double &d)
divide by a double
bool Built
indicates that the vector has been built and is usable
double * Values_pt
the local vector
double dot(const DoubleVector &vec) const
compute the dot product of this vector with the vector vec.
double & operator[](int i)
[] access function to the (local) values of this vector
double * values_pt()
access function to the underlying values
void clear()
wipes the DoubleVector
void output_local_values_with_offset(std::ostream &outfile, const int &output_precision=-1) const
output the local contents of the vector
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
bool distributed() const
access function to the distributed - indicates whether the distribution is serial or distributed
unsigned first_row() const
access function for the first row on this processor. If not distributed then this is just zero.
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
bool built() const
if the communicator_pt is null then the distribution is not setup then false is returned,...
unsigned nrow() const
access function to the number of global rows.
unsigned nrow_local() const
access function for the num of local rows on this processor. If no MPI then Nrow is returned.
unsigned rank_of_global_row(const unsigned i) const
return the processor rank of the global row number i
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:54
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
void split(const DoubleVector &in_vector, Vector< DoubleVector * > &out_vector_pt)
Split a DoubleVector into the out DoubleVectors. Let vec_A be the in Vector, and let vec_B and vec_C ...
void concatenate_without_communication(const Vector< DoubleVector * > &in_vector_pt, DoubleVector &out_vector)
Concatenate DoubleVectors. Takes a Vector of DoubleVectors. If the out vector is built,...
void split_without_communication(const DoubleVector &in_vector, Vector< DoubleVector * > &out_vector_pt)
Split a DoubleVector into the out DoubleVectors. Data stays on its current processor,...
void concatenate(const Vector< DoubleVector * > &in_vector_pt, DoubleVector &out_vector)
Concatenate DoubleVectors. Takes a Vector of DoubleVectors. If the out vector is built,...
void concatenate(const Vector< LinearAlgebraDistribution * > &in_distribution_pt, LinearAlgebraDistribution &out_distribution)
Takes a vector of LinearAlgebraDistribution objects and concatenates them such that the nrow_local of...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
std::ostream & operator<<(std::ostream &out, const DoubleVector &v)
output operator