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