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