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-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26#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
46namespace 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
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],
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
unsigned nrow() const
access function to the number of global rows.
bool distribution_built() const
if the communicator_pt is null then the distribution is not setup then false is returned,...
unsigned nrow_local() const
access function for the num of local rows on this processor.
unsigned first_row() const
access function for the first row on this processor
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
void build_distribution(const LinearAlgebraDistribution *const dist_pt)
setup the distribution of this distributable linear algebra object
A 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
const DoubleVector & doublevector(const unsigned &i) const
access to the DoubleVector representation (const version)
void output(std::ostream &outfile) const
output the contents of the vector
double * values(const unsigned &i) const
access function to the i-th vector's data (const version)
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...
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...
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.
double ** values()
access function to the underlying values
~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
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)
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
double & operator()(int v, int i) const
[] access function to the (local) values of the v-th vector
void shallow_build(const DoubleMultiVector &old_vector)
Provide a (shallow) copy of the old vector.
bool operator==(const DoubleMultiVector &vec)
== operator
double ** values() const
access function to the underlying values (const version)
void operator-=(DoubleMultiVector vec)
-= operator
void operator+=(DoubleMultiVector vec)
+= operator
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.
double * values(const unsigned &i)
access function to the i-th vector's data
DoubleMultiVector(const DoubleMultiVector &new_vector)
Copy constructor.
DoubleVector & doublevector(const unsigned &i)
access to the DoubleVector representatoin
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.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...