double_multi_vector.h
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: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 #ifndef OOMPH_DOUBLE_MULTI_VECTOR_CLASS_HEADER
31 #define OOMPH_DOUBLE_MULTI_VECTOR_CLASS_HEADER
32 
33 // Config header generated by autoconfig
34 #ifdef HAVE_CONFIG_H
35  #include <oomph-lib-config.h>
36 #endif
37 
38 #ifdef OOMPH_HAS_MPI
39 #include "mpi.h"
40 #endif
41 
42 //oomph headers
43 #include "double_vector.h"
44 
45 //Trilinos headerss
46 #ifdef OOMPH_HAS_TRILINOS
47 #include "Teuchos_Range1D.hpp"
48 #endif
49 
50 namespace oomph
51 {
52 
53 //=============================================================================
54 /// \short A multi vector in the mathematical sense, initially developed for
55 /// linear algebra type applications.
56 /// If MPI then this multi vector can be distributed - its distribution is
57 /// described by the LinearAlgebraDistribution object at Distribution_pt.
58 /// Data is stored in a C-style pointer vector (double*)
59 //=============================================================================
61 {
62 
63  public :
64 
65  /// \short Constructor for an uninitialized DoubleMultiVector
67  : Values(0), Nvector(0), Internal_values(true), Built(false) {}
68 
69  /// \short Constructor. Assembles a DoubleMultiVector consisting of
70  /// n_vector vectors, each with a prescribed distribution.
71  /// Additionally every entry can be set (with argument v -
72  /// defaults to 0).
73  DoubleMultiVector(const unsigned &n_vector,
74  const LinearAlgebraDistribution* const &dist_pt,
75  const double& v=0.0)
76  : Values(0), Nvector(n_vector), Internal_values(true), Built(false)
77  {
78  this->build(n_vector,dist_pt,v);
80  }
81 
82  /// \short Constructor. Assembles a DoubleMultiVector consisting of
83  /// n_vector vectors, each with a prescribed distribution.
84  /// Additionally every entry can be set (with argument v -
85  /// defaults to 0).
86  DoubleMultiVector(const unsigned &n_vector,
87  const LinearAlgebraDistribution& dist,
88  const double& v=0.0)
89  : Values(0), Nvector(n_vector), Internal_values(true), Built(false)
90  {
91  this->build(n_vector,dist,v);
93  }
94 
95  /// \short Constructor. Build a multivector using the same distribution
96  /// of the input vector with n_vector columns and initialised to the value
97  /// v
98  DoubleMultiVector(const unsigned &n_vector,
99  const DoubleMultiVector &old_vector,
100  const double &initial_value=0.0) :
101  Values(0), Nvector(n_vector), Internal_values(true), Built(false)
102  {
103  this->build(n_vector,old_vector.distribution_pt(),initial_value);
105  }
106 
107  /// \short Constructor that builds a multivector from selected columns
108  /// of the input multivector. The boolean controls whether it is a
109  /// shallow or deep copy (default deep)
111  const std::vector<int> &index,
112  const bool &deep_copy=true) :
113  Values(0), Nvector(0), Internal_values(deep_copy), Built(false)
114  {
115  //Build the storage based on the size of index
116  unsigned n_vector = index.size();
117  if(deep_copy)
118  {
119  //Create an entirely new data structure
120  this->build(n_vector,old_vector.distribution_pt());
121  //Now (deep) copy the data across
122  const unsigned n_row_local = this->nrow_local();
123  for(unsigned v=0;v<n_vector;v++)
124  {
125  for(unsigned i=0;i<n_row_local;i++)
126  {
127  Values[v][i] = old_vector(index[v],i);
128  }
129  }
130  }
131  //Otherwise it's a shallow copy
132  else
133  {
134  this->shallow_build(n_vector,old_vector.distribution_pt());
135  //Now shallow copy the pointers accross
136  for(unsigned v=0;v<n_vector;++v)
137  {
138  Values[v] = old_vector.values(index[v]);
139  }
140  }
142  }
143 
144 #ifdef OOMPH_HAS_TRILINOS
145  /// \short Constructor that builds a multivector from selected columns
146  /// of the input multivector and the provided range. The optional
147  /// boolean specifies whether it is a shallow or deep copy. The default
148  /// is that it is a deep copy.
150  const Teuchos::Range1D &index,
151  const bool &deep_copy=true) :
152  Values(0), Nvector(0), Internal_values(deep_copy), Built(false)
153  {
154  //Build the storage based on the size of index
155  unsigned n_vector = index.size();
156  if(deep_copy)
157  {
158  //Create entirely new data structure
159  this->build(n_vector,old_vector.distribution_pt());
160  //Now (deep) copy the data across
161  const unsigned n_row_local = this->nrow_local();
162  unsigned range_index = index.lbound();
163  for(unsigned v=0;v<n_vector;v++)
164  {
165  for(unsigned i=0;i<n_row_local;i++)
166  {
167  Values[v][i] = old_vector(range_index,i);
168  }
169  ++range_index;
170  }
171  }
172  //Otherwise it's a shallow copy
173  else
174  {
175  this->shallow_build(n_vector,old_vector.distribution_pt());
176  //Shallow copy the pointers accross
177  unsigned range_index = index.lbound();
178  for(unsigned v=0;v<n_vector;v++)
179  {
180  Values[v] = old_vector.values(range_index);
181  ++range_index;
182  }
183  }
185  }
186 #endif
187 
188  /// Copy constructor
191  Values(0), Nvector(0), Internal_values(true), Built(false)
192  {
193  this->build(new_vector);
195  }
196 
197 
198  /// Destructor - just calls this->clear() to delete the distribution and data
200  {
201  this->clear();
202  }
203 
204  /// assignment operator (deep copy)
205  void operator=(const DoubleMultiVector& old_vector)
206  {
207  this->build(old_vector);
209  }
210 
211  /// Return the number of vectors
212  unsigned nvector() const {return Nvector;}
213 
214  /// \short Provide a (shallow) copy of the old vector
215  void shallow_build(const DoubleMultiVector& old_vector)
216  {
217  //Only bother if the old_vector is not the same as current vector
218  if (!(*this == old_vector))
219  {
220  // the vector does not own the internal data
221  Internal_values = false;
222 
223  //Copy the number of vectors
224  Nvector = old_vector.nvector();
225  //Allocate storage for pointers to the values
226  this->shallow_build(Nvector,old_vector.distribution_pt());
227 
228  // copy all the pointers accross
229  if (this->distribution_built())
230  {
231  for(unsigned v=0;v<Nvector;++v)
232  {
233  Values[v] = old_vector.values(v);
234  }
235  }
236  }
237  }
238 
239  /// \short Build the storage for pointers to vectors with a given
240  /// distribution, but do not populate the pointers
241  void shallow_build(const unsigned &n_vector,
242  const LinearAlgebraDistribution &dist)
243  {
244  this->shallow_build(n_vector,&dist);
245  }
246 
247 
248  /// \short Build the storage for pointers to vectors with a given
249  /// distribution, but do not populate the pointers
250  void shallow_build(const unsigned &n_vector,
251  const LinearAlgebraDistribution* const &dist_pt)
252  {
253  // clean the memory
254  this->clear();
255 
256  //The vector does not own the data
257  Internal_values = false;
258 
259  // Set the distribution
260  this->build_distribution(dist_pt);
261 
262  // update the values
263  if (dist_pt->built())
264  {
265  //Set the number of vectors
266  Nvector = n_vector;
267  //Build the array of pointers to each vector's data
268  Values = new double*[n_vector];
269  Built=true;
270  }
271  else
272  {
273  Built=false;
274  }
275  }
276 
277 
278  /// \short Provides a (deep) copy of the old_vector
279  void build(const DoubleMultiVector& old_vector)
280  {
281  //Only bother if the old_vector is not the same as current vector
282  if (!(*this == old_vector))
283  {
284  // the vector owns the internal data
285  Internal_values = true;
286 
287  //Copy the number of vectors
288  Nvector = old_vector.nvector();
289  // reset the distribution and resize the data
290  this->build(Nvector,old_vector.distribution_pt(),0.0);
291 
292  // copy the data
293  if (this->distribution_built())
294  {
295  unsigned n_row_local = this->nrow_local();
296  double** const old_vector_values = old_vector.values();
297  for (unsigned i = 0; i < n_row_local; i++)
298  {
299  for(unsigned v=0;v<Nvector;v++)
300  {
301  Values[v][i] = old_vector_values[v][i];
302  }
303  }
304  }
305  }
306  }
307 
308  /// \short Assembles a DoubleMultiVector
309  /// with n_vector vectors, a distribution dist, if v is specified
310  /// each element is set to v, otherwise each element is set to 0.0
311  void build(const unsigned &n_vector,
312  const LinearAlgebraDistribution& dist,
313  const double& initial_value=0.0)
314  {
315  this->build(n_vector,&dist,initial_value);
316  }
317 
318  /// \short Assembles a DoubleMultiVector with n_vector vectors, each with a
319  /// distribution dist, if v is specified
320  /// each element is set to v, otherwise each element is set to 0.0
321  void build(const unsigned &n_vector,
322  const LinearAlgebraDistribution* const &dist_pt,
323  const double& initial_value=0.0)
324  {
325  // clean the memory
326  this->clear();
327 
328  // the vector owns the internal data
329  Internal_values = true;
330 
331  // Set the distribution
332  this->build_distribution(dist_pt);
333 
334  // update the values
335  if (dist_pt->built())
336  {
337  //Set the number of vectors
338  Nvector = n_vector;
339  //Build the array of pointers to each vector's data
340  Values = new double*[n_vector];
341  //Now build the contiguous array of real data
342  const unsigned n_row_local = this->nrow_local();
343  double *values = new double[n_row_local*n_vector];
344  // set the data
345  for(unsigned v=0;v<n_vector;v++)
346  {
347  Values[v] = &values[v*n_row_local];
348  for (unsigned i = 0; i < n_row_local; i++)
349  {
350  Values[v][i] = initial_value;
351  }
352  }
353  Built=true;
354  }
355  else
356  {
357  Built=false;
358  }
359  }
360 
361  /// \short initialise the whole vector with value v
362  void initialise(const double& initial_value)
363  {
364  if (Built)
365  {
366  const unsigned n_vector = this->Nvector;
367  const unsigned n_row_local = this->nrow_local();
368 
369  // set the residuals
370  for(unsigned v=0;v<n_vector;v++)
371  {
372  for (unsigned i=0; i<n_row_local; i++)
373  {
374  Values[v][i] = initial_value;
375  }
376  }
377  }
378  }
379 
380  /// \short initialise the vector with coefficient from the vector v.
381  /// Note: The vector v must be of length
382  //void initialise(const Vector<double> v);
383 
384  /// \short wipes the DoubleVector
385  void clear()
386  {
387  //Return if nothing to do
388  if(Values==0) {return;}
389 
390  //If we are in charge of the data then delete it
391  if (Internal_values)
392  {
393  //Delete the double storage arrays at once
394  //(they were allocated as a contiguous block)
395  delete[] Values[0];
396  }
397 
398  //Always Delete the pointers to the arrays
399  delete[] Values;
400 
401  //Then set the pointer (to a pointer) to null
402  Values = 0;
403  this->clear_distribution();
404  Built=false;
405  }
406 
407  // indicates whether this DoubleVector is built
408  bool built() const { return Built; }
409 
410  /// \short Allows are external data to be used by this vector.
411  /// WARNING: The size of the external data must correspond to the
412  /// LinearAlgebraDistribution dist_pt argument.
413  /// 1. When a rebuild method is called new internal values are created.
414  /// 2. It is not possible to redistribute(...) a vector with external
415  /// values .
416  /// 3. External values are only deleted by this vector if
417  /// delete_external_values = true.
418  /*void set_external_values(const LinearAlgebraDistribution* const& dist_pt,
419  double* external_values,
420  bool delete_external_values)
421  {
422  // clean the memory
423  this->clear();
424 
425  // Set the distribution
426  this->build_distribution(dist_pt);
427 
428  // set the external values
429  set_external_values(external_values,delete_external_values);
430  }*/
431 
432  /// \short Allows are external data to be used by this vector.
433  /// WARNING: The size of the external data must correspond to the
434  /// distribution of this vector.
435  /// 1. When a rebuild method is called new internal values are created.
436  /// 2. It is not possible to redistribute(...) a vector with external
437  /// values .
438  /// 3. External values are only deleted by this vector if
439  /// delete_external_values = true.
440  /*void set_external_values(double* external_values,
441  bool delete_external_values)
442  {
443 #ifdef PARANOID
444  // check that this distribution is setup
445  if (!this->distribution_built())
446  {
447  // if this vector does not own the double* values then it cannot be
448  // distributed.
449  // note: this is not stictly necessary - would just need to be careful
450  // with delete[] below.
451  std::ostringstream error_message;
452  error_message << "The distribution of the vector must be setup before "
453  << "external values can be set";
454  throw OomphLibError(error_message.str(),
455  OOMPH_CURRENT_FUNCTION,
456  OOMPH_EXCEPTION_LOCATION);
457  }
458 #endif
459  if (Internal_values)
460  {
461  delete[] Values_pt;
462  }
463  Values_pt = external_values;
464  Internal_values = delete_external_values;
465  }*/
466 
467  /// \short The contents of the vector are redistributed to match the new
468  /// distribution. In a non-MPI rebuild this method works, but does nothing.
469  /// \b NOTE 1: The current distribution and the new distribution must have
470  /// the same number of global rows.
471  /// \b NOTE 2: The current distribution and the new distribution must have
472  /// the same Communicator.
473  void redistribute(const LinearAlgebraDistribution* const& dist_pt);
474 
475  /// \short [] access function to the (local) values of the v-th vector
476  double& operator()(int v, int i) const
477  {
478 #ifdef RANGE_CHECKING
479  std::ostringstream error_message;
480  bool error=false;
481  if(v > int(Nvector))
482  {
483  error_message << "Range Error: Vector " << v
484  << "is not in the range (0,"
485  << Nvector-1 << ")";
486  error = true;
487  }
488 
489  if (i >= int(this->nrow_local()))
490  {
491  error_message << "Range Error: " << i
492  << " is not in the range (0,"
493  << this->nrow_local()-1 << ")";
494  error = true;
495  }
496 
497  if(error)
498  {
499  throw OomphLibError(error_message.str(),
500  OOMPH_CURRENT_FUNCTION,
501  OOMPH_EXCEPTION_LOCATION);
502  }
503 #endif
504  return Values[v][i];
505  }
506 
507  /// \short == operator
508  bool operator==(const DoubleMultiVector& vec)
509  {
510  // if vec is not setup return false
511  if (vec.built() && !this->built())
512  {
513  return false;
514  }
515  else if (!vec.built() && this->built())
516  {
517  return false;
518  }
519  else if (!vec.built() && !this->built())
520  {
521  return true;
522  }
523  else
524  {
525  double** const v_values = vec.values();
526  const unsigned n_row_local = this->nrow_local();
527  const unsigned n_vector = this->nvector();
528  for(unsigned v=0;v<n_vector;++v)
529  {
530  for(unsigned i=0;i<n_row_local;++i)
531  {
532  if (Values[v][i] != v_values[v][i])
533  {
534  return false;
535  }
536  }
537  }
538  return true;
539  }
540  }
541 
542  /// \short += operator
544  {
545 #ifdef PARANOID
546  // PARANOID check that this vector is setup
547  if (!this->built())
548  {
549  std::ostringstream error_message;
550  error_message << "This vector must be setup.";
551  throw OomphLibError(error_message.str(),
552  OOMPH_CURRENT_FUNCTION,
553  OOMPH_EXCEPTION_LOCATION);
554  }
555  // PARANOID check that the vector v is setup
556  if (!vec.built())
557  {
558  std::ostringstream error_message;
559  error_message << "The vector v must be setup.";
560  throw OomphLibError(error_message.str(),
561  OOMPH_CURRENT_FUNCTION,
562  OOMPH_EXCEPTION_LOCATION);
563  }
564  // PARANOID check that the vectors have the same distribution
565  if (!(*vec.distribution_pt() == *this->distribution_pt()))
566  {
567  std::ostringstream error_message;
568  error_message << "The vector v and this vector must have the same "
569  << "distribution.";
570  throw OomphLibError(error_message.str(),
571  OOMPH_CURRENT_FUNCTION,
572  OOMPH_EXCEPTION_LOCATION);
573  }
574 #endif//
575 
576 
577  double** v_values = vec.values();
578  const unsigned n_vector = this->nvector();
579  const unsigned n_row_local = this->nrow_local();
580  for(unsigned v=0; v<n_vector;++v)
581  {
582  for(unsigned i=0; i<n_row_local;++i)
583  {
584  Values[v][i] += v_values[v][i];
585  }
586  }
587  }
588 
589  /// -= operator
591  {
592 #ifdef PARANOID
593  // PARANOID check that this vector is setup
594  if (!this->distribution_built())
595  {
596  std::ostringstream error_message;
597  error_message << "This vector must be setup.";
598  throw OomphLibError(error_message.str(),
599  OOMPH_CURRENT_FUNCTION,
600  OOMPH_EXCEPTION_LOCATION);
601  }
602  // PARANOID check that the vector v is setup
603  if (!vec.built())
604  {
605  std::ostringstream error_message;
606  error_message << "The vector v must be setup.";
607  throw OomphLibError(error_message.str(),
608  OOMPH_CURRENT_FUNCTION,
609  OOMPH_EXCEPTION_LOCATION);
610  }
611  // PARANOID check that the vectors have the same distribution
612  if (!(*vec.distribution_pt() == *this->distribution_pt()))
613  {
614  std::ostringstream error_message;
615  error_message << "The vector v and this vector must have the same "
616  << "distribution.";
617  throw OomphLibError(error_message.str(),
618  OOMPH_CURRENT_FUNCTION,
619  OOMPH_EXCEPTION_LOCATION);
620  }
621 #endif
622 
623  double** v_values = vec.values();
624  const unsigned n_vector = this->nvector();
625  const unsigned n_row_local = this->nrow_local();
626  for(unsigned v=0;v<n_vector;++v)
627  {
628  for(unsigned i=0;i<n_row_local;++i)
629  {
630  Values[v][i] -= v_values[v][i];
631  }
632  }
633  }
634 
635  ///Multiply by a scalar
636  void operator *=(const double& scalar_value)
637  {
638  #ifdef PARANOID
639  // PARANOID check that this vector is setup
640  if (!this->distribution_built())
641  {
642  std::ostringstream error_message;
643  error_message << "This vector must be setup.";
644  throw OomphLibError(error_message.str(),
645  OOMPH_CURRENT_FUNCTION,
646  OOMPH_EXCEPTION_LOCATION);
647  }
648 #endif
649  const unsigned n_vector = this->nvector();
650  const unsigned n_row_local = this->nrow_local();
651  for(unsigned v=0;v<n_vector;++v)
652  {
653  for(unsigned i=0;i<n_row_local;++i)
654  {
655  Values[v][i] *= scalar_value;
656  }
657  }
658  }
659 
660 
661  /// access function to the underlying values
662  double** values() {return Values;}
663 
664  /// \short access function to the underlying values (const version)
665  double** values() const {return Values;}
666 
667  /// access function to the i-th vector's data
668  double* values(const unsigned &i) {return Values[i];}
669 
670  /// access function to the i-th vector's data (const version)
671  double* values(const unsigned &i) const {return Values[i];}
672 
673  /// access to the DoubleVector representatoin
674  DoubleVector &doublevector(const unsigned &i)
675  {return Internal_doublevector[i];}
676 
677  /// access to the DoubleVector representation (const version)
678  const DoubleVector &doublevector(const unsigned &i) const
679  {return Internal_doublevector[i];}
680 
681  /// output the contents of the vector
682  void output(std::ostream &outfile) const
683  {
684  // temp pointer to values
685  double** temp;
686 
687  // Number of vectors
688  unsigned n_vector = this->nvector();
689 
690  // number of global row
691  unsigned nrow = this->nrow();
692 
693 #ifdef OOMPH_HAS_MPI
694 
695  // number of local rows
696  int nrow_local = this->nrow_local();
697 
698  // gather from all processors
699  if (this->distributed() &&
700  this->distribution_pt()->communicator_pt()->nproc() > 1)
701  {
702  // number of processors
703  int nproc = this->distribution_pt()->communicator_pt()->nproc();
704 
705  // number of gobal row
706  unsigned nrow = this->nrow();
707 
708  // get the vector of first_row s and nrow_local s
709  int* dist_first_row = new int[nproc];
710  int* dist_nrow_local = new int[nproc];
711  for (int p = 0; p < nproc; p++)
712  {
713  dist_first_row[p] = this->first_row(p);
714  dist_nrow_local[p] = this->nrow_local(p);
715  }
716 
717  // gather
718  temp = new double*[n_vector];
719  double* temp_value = new double[nrow*n_vector];
720  for(unsigned v=0;v<n_vector;v++)
721  {
722  temp[v] = &temp_value[v*nrow];
723  }
724 
725  //Now do an all gather for each vector
726  //Possibly costly in terms of extra communication, but it's only output!
727  for(unsigned v=0;v<n_vector;++v)
728  {
729  MPI_Allgatherv(Values[v],nrow_local,MPI_DOUBLE,
730  temp[v],dist_nrow_local,dist_first_row,
731  MPI_DOUBLE,
732  this->distribution_pt()->communicator_pt()->mpi_comm());
733  }
734 
735  // clean up
736  delete[] dist_first_row;
737  delete[] dist_nrow_local;
738  }
739  else
740  {
741  temp = Values;
742  }
743 #else
744  temp = Values;
745 #endif
746 
747  // output
748  for (unsigned i = 0; i < nrow; i++)
749  {
750  outfile << i << " ";
751  for(unsigned v = 0; v<n_vector; v++)
752  {
753  outfile << temp[v][i] << " ";
754  }
755  outfile << "\n";
756  }
757 
758  // clean up if required
759 #ifdef OOMPH_HAS_MPI
760  if (this->distributed() &&
761  this->distribution_pt()->communicator_pt()->nproc() > 1)
762  {
763  delete[] temp[0];
764  delete[] temp;
765  }
766 #endif
767 
768  }
769 
770  /// output the contents of the vector
771  void output(std::string filename)
772  {
773  // Open file
774  std::ofstream some_file;
775  some_file.open(filename.c_str());
776  output(some_file);
777  some_file.close();
778  }
779 
780 
781  /// compute the 2 norm of this vector
782  void dot(const DoubleMultiVector& vec, std::vector<double> &result) const
783  {
784 #ifdef PARANOID
785  // paranoid check that the vector is setup
786  if (!this->built())
787  {
788  std::ostringstream error_message;
789  error_message << "This vector must be setup.";
790  throw OomphLibError(error_message.str(),
791  OOMPH_CURRENT_FUNCTION,
792  OOMPH_EXCEPTION_LOCATION);
793  }
794  if (!vec.built())
795  {
796  std::ostringstream error_message;
797  error_message << "The input vector be setup.";
798  throw OomphLibError(error_message.str(),
799  OOMPH_CURRENT_FUNCTION,
800  OOMPH_EXCEPTION_LOCATION);
801  }
802  if (*this->distribution_pt() != *vec.distribution_pt())
803  {
804  std::ostringstream error_message;
805  error_message << "The distribution of this vector and the vector vec "
806  << "must be the same."
807  << "\n\n this: " << *this->distribution_pt()
808  << "\n vec: " << *vec.distribution_pt();
809  throw OomphLibError(error_message.str(),
810  OOMPH_CURRENT_FUNCTION,
811  OOMPH_EXCEPTION_LOCATION);
812  }
813 #endif
814 
815  // compute the local norm
816  unsigned nrow_local = this->nrow_local();
817  int n_vector = this->nvector();
818  double n[n_vector];
819  for(int v=0;v<n_vector;v++)
820  {
821  //Initialise
822  n[v] = 0.0;
823  const double* vec_values_pt = vec.values(v);
824  for (unsigned i = 0; i < nrow_local; i++)
825  {
826  n[v] += Values[v][i]*vec_values_pt[i];
827  }
828  }
829 
830  // if this vector is distributed and on multiple processors then gather
831 #ifdef OOMPH_HAS_MPI
832  double n2[n_vector]; for(int v=0;v<n_vector;v++) {n2[v] = n[v];}
833 
834  if (this->distributed() &&
835  this->distribution_pt()->communicator_pt()->nproc() > 1)
836  {
837  MPI_Allreduce(&n,&n2,n_vector,MPI_DOUBLE,MPI_SUM,
838  this->distribution_pt()->communicator_pt()->mpi_comm());
839  }
840  for(int v=0;v<n_vector;v++) {n[v] = n2[v];}
841 #endif
842 
843  result.resize(n_vector);
844  for(int v=0;v<n_vector;v++) {result[v] = n[v];}
845  }
846 
847  /// compute the 2 norm of this vector
848  void norm(std::vector<double> &result) const
849  {
850 #ifdef PARANOID
851  // paranoid check that the vector is setup
852  if (!this->built())
853  {
854  std::ostringstream error_message;
855  error_message << "This vector must be setup.";
856  throw OomphLibError(error_message.str(),
857  OOMPH_CURRENT_FUNCTION,
858  OOMPH_EXCEPTION_LOCATION);
859  }
860 #endif
861 
862  // compute the local norm
863  unsigned nrow_local = this->nrow_local();
864  int n_vector = this->nvector();
865  double n[n_vector];
866  for(int v=0;v<n_vector;v++)
867  {
868  n[v] = 0.0;
869  for (unsigned i = 0; i < nrow_local; i++)
870  {
871  n[v] += Values[v][i]*Values[v][i];
872  }
873  }
874 
875  // if this vector is distributed and on multiple processors then gather
876 #ifdef OOMPH_HAS_MPI
877  double n2[n_vector]; for(int v=0;v<n_vector;v++) {n2[v] = n[v];}
878  if (this->distributed() &&
879  this->distribution_pt()->communicator_pt()->nproc() > 1)
880  {
881  MPI_Allreduce(&n,&n2,n_vector,MPI_DOUBLE,MPI_SUM,
882  this->distribution_pt()->communicator_pt()->mpi_comm());
883  }
884  for(int v=0;v<n_vector;v++) {n[v] = n2[v];}
885 #endif
886 
887  //Now sqrt the norm and fill in result
888  result.resize(n_vector);
889  for(int v=0;v<n_vector;v++) {result[v] = sqrt(n[v]);}
890  }
891 
892  /// compute the A-norm using the matrix at matrix_pt
893  /*double norm(const CRDoubleMatrix* matrix_pt) const
894  {
895 #ifdef PARANOID
896  // paranoid check that the vector is setup
897  if (!this->built())
898  {
899  std::ostringstream error_message;
900  error_message << "This vector must be setup.";
901  throw OomphLibError(error_message.str(),
902  OOMPH_CURRENT_FUNCTION,
903  OOMPH_EXCEPTION_LOCATION);
904  }
905  if (!matrix_pt->built())
906  {
907  std::ostringstream error_message;
908  error_message << "The input matrix be built.";
909  throw OomphLibError(error_message.str(),
910  OOMPH_CURRENT_FUNCTION,
911  OOMPH_EXCEPTION_LOCATION);
912  }
913  if (*this->distribution_pt() != *matrix_pt->distribution_pt())
914  {
915  std::ostringstream error_message;
916  error_message << "The distribution of this vector and the matrix at "
917  << "matrix_pt must be the same";
918  throw OomphLibError(error_message.str(),
919  OOMPH_CURRENT_FUNCTION,
920  OOMPH_EXCEPTION_LOCATION);
921  }
922 #endif
923 
924  // compute the matrix norm
925  DoubleVector x(this->distribution_pt(),0.0);
926  matrix_pt->multiply(*this,x);
927  return sqrt(this->dot(x));
928 
929  }*/
930 
931  private :
932 
933  /// Setup the doublevector representation
935  {
936  const unsigned n_vector = this->nvector();
937  Internal_doublevector.resize(n_vector);
938  //Loop over all and set the external values of the DoubleVectors
939  for(unsigned v=0;v<n_vector;v++)
940  {
941  Internal_doublevector[v].set_external_values(this->distribution_pt(),
942  this->values(v),false);
943  }
944  }
945 
946  /// the local data, need a pointer to a pointer so that the
947  /// individual vectors can be extracted
948  double** Values;
949 
950  /// The number of vectors
951  unsigned Nvector;
952 
953  /// \short Boolean flag to indicate whether the vector's data (values_pt)
954  /// is owned by this vector.
956 
957  /// indicates that the vector has been built and is usable
958  bool Built;
959 
960  /// Need a vector of DoubleVectors to interface with our linear solvers
962 
963 }; //end of DoubleVector
964 
965 } // end of oomph namespace
966 
967 
968 
969 #endif
unsigned nrow_local() const
access function for the num of local rows on this processor.
void shallow_build(const unsigned &n_vector, const LinearAlgebraDistribution *const &dist_pt)
Build the storage for pointers to vectors with a given distribution, but do not populate the pointers...
unsigned Nvector
The number of vectors.
void operator+=(DoubleMultiVector vec)
+= operator
unsigned first_row() const
access function for the first row on this processor
bool operator==(const DoubleMultiVector &vec)
== operator
void redistribute(const LinearAlgebraDistribution *const &dist_pt)
Allows are external data to be used by this vector. WARNING: The size of the external data must corre...
double ** values()
access function to the underlying values
DoubleMultiVector()
Constructor for an uninitialized DoubleMultiVector.
cstr elem_len * i
Definition: cfortran.h:607
void shallow_build(const unsigned &n_vector, const LinearAlgebraDistribution &dist)
Build the storage for pointers to vectors with a given distribution, but do not populate the pointers...
DoubleMultiVector(const unsigned &n_vector, const LinearAlgebraDistribution *const &dist_pt, const double &v=0.0)
Constructor. Assembles a DoubleMultiVector consisting of n_vector vectors, each with a prescribed dis...
unsigned nrow() const
access function to the number of global rows.
double & operator()(int v, int i) const
[] access function to the (local) values of the v-th vector
void dot(const DoubleMultiVector &vec, std::vector< double > &result) const
compute the 2 norm of this vector
DoubleMultiVector(const DoubleMultiVector &new_vector)
Copy constructor.
bool distributed() const
distribution is serial or distributed
void clear()
initialise the vector with coefficient from the vector v. Note: The vector v must be of length ...
void setup_doublevector_representation()
compute the A-norm using the matrix at matrix_pt
void operator*=(const double &scalar_value)
Multiply by a scalar.
void build(const unsigned &n_vector, const LinearAlgebraDistribution *const &dist_pt, const double &initial_value=0.0)
Assembles a DoubleMultiVector with n_vector vectors, each with a distribution dist, if v is specified each element is set to v, otherwise each element is set to 0.0.
DoubleMultiVector(const DoubleMultiVector &old_vector, const Teuchos::Range1D &index, const bool &deep_copy=true)
Constructor that builds a multivector from selected columns of the input multivector and the provided...
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
bool Built
indicates that the vector has been built and is usable
void output(std::ostream &outfile) const
output the contents of the vector
DoubleVector & doublevector(const unsigned &i)
access to the DoubleVector representatoin
DoubleMultiVector(const DoubleMultiVector &old_vector, const std::vector< int > &index, const bool &deep_copy=true)
Constructor that builds a multivector from selected columns of the input multivector. The boolean controls whether it is a shallow or deep copy (default deep)
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
DoubleMultiVector(const unsigned &n_vector, const DoubleMultiVector &old_vector, const double &initial_value=0.0)
Constructor. Build a multivector using the same distribution of the input vector with n_vector column...
void initialise(const double &initial_value)
initialise the whole vector with value v
void deep_copy(const CRDoubleMatrix *const in_matrix_pt, CRDoubleMatrix &out_matrix)
Create a deep copy of the matrix pointed to by in_matrix_pt.
Definition: matrices.h:3244
Base class for any linear algebra object that is distributable. Just contains storage for the LinearA...
double * values(const unsigned &i) const
access function to the i-th vector's data (const version)
void output(std::string filename)
output the contents of the vector
~DoubleMultiVector()
Destructor - just calls this->clear() to delete the distribution and data.
double ** values() const
access function to the underlying values (const version)
void build(const unsigned &n_vector, const LinearAlgebraDistribution &dist, const double &initial_value=0.0)
Assembles a DoubleMultiVector with n_vector vectors, a distribution dist, if v is specified each elem...
void clear_distribution()
clear the distribution of this distributable linear algebra object
DoubleMultiVector(const unsigned &n_vector, const LinearAlgebraDistribution &dist, const double &v=0.0)
Constructor. Assembles a DoubleMultiVector consisting of n_vector vectors, each with a prescribed dis...
bool Internal_values
Boolean flag to indicate whether the vector's data (values_pt) is owned by this vector.
void shallow_build(const DoubleMultiVector &old_vector)
Provide a (shallow) copy of the old vector.
void norm(std::vector< double > &result) const
compute the 2 norm of this vector
void build(const DoubleMultiVector &old_vector)
Provides a (deep) copy of the old_vector.
A multi vector in the mathematical sense, initially developed for linear algebra type applications...
void operator=(const DoubleMultiVector &old_vector)
assignment operator (deep copy)
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
Vector< DoubleVector > Internal_doublevector
Need a vector of DoubleVectors to interface with our linear solvers.
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
const DoubleVector & doublevector(const unsigned &i) const
access to the DoubleVector representation (const version)
void operator-=(DoubleMultiVector vec)
-= operator
unsigned nvector() const
Return the number of vectors.
double * values(const unsigned &i)
access function to the i-th vector's data