double_multi_vector.h

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented, 
//LIC// multi-physics finite-element library, available 
//LIC// at http://www.oomph-lib.org.
//LIC// 
//LIC//    Version 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
//LIC// 
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC// 
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//LIC// Lesser General Public License for more details.
//LIC// 
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301  USA.
//LIC// 
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC// 
//LIC//====================================================================
#ifndef OOMPH_DOUBLE_MULTI_VECTOR_CLASS_HEADER
#define OOMPH_DOUBLE_MULTI_VECTOR_CLASS_HEADER

// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
  #include <oomph-lib-config.h>
#endif

#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif

//oomph headers
#include "double_vector.h"

//Trilinos headerss
#ifdef OOMPH_HAS_TRILINOS
#include "Teuchos_Range1D.hpp"
#endif

namespace oomph
{

//=============================================================================
/// \short A multi vector in the mathematical sense, initially developed for
/// linear algebra type applications.
/// If MPI then this multi 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*)
//=============================================================================
 class DoubleMultiVector : public DistributableLinearAlgebraObject
{                                                        
 
 public :                     
  
  /// \short Constructor for an uninitialized DoubleMultiVector
  DoubleMultiVector()
   : Values(0), Nvector(0), Internal_values(true), Built(false) {}
 
 /// \short Constructor. Assembles a DoubleMultiVector consisting of 
 /// n_vector vectors, each with a prescribed distribution.
 /// Additionally every entry can be set (with argument v - 
 /// defaults to 0).
 DoubleMultiVector(const unsigned &n_vector,
                   const LinearAlgebraDistribution* const &dist_pt, 
                   const double& v=0.0)
  : Values(0), Nvector(n_vector), Internal_values(true), Built(false)
  {
   this->build(n_vector,dist_pt,v);
   this->setup_doublevector_representation();
  }

 /// \short Constructor. Assembles a DoubleMultiVector consisting of 
 /// n_vector vectors, each with a prescribed distribution. 
 /// Additionally every entry can be set (with argument v - 
 /// defaults to 0).
 DoubleMultiVector(const unsigned &n_vector,
                   const LinearAlgebraDistribution& dist, 
                   const double& v=0.0)
  : Values(0), Nvector(n_vector), Internal_values(true), Built(false)
  {
   this->build(n_vector,dist,v);
   this->setup_doublevector_representation();
  }
   
 /// \short Constructor. Build a multivector using the same distribution
 /// of the input vector with n_vector columns and initialised to the value
 /// v
 DoubleMultiVector(const unsigned &n_vector,
                   const DoubleMultiVector &old_vector,
                   const double &initial_value=0.0) :
  Values(0), Nvector(n_vector), Internal_values(true), Built(false)
  {
   this->build(n_vector,old_vector.distribution_pt(),initial_value);
   this->setup_doublevector_representation();
  }

 /// \short 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)
 DoubleMultiVector(const DoubleMultiVector&old_vector,
                   const std::vector<int> &index,
                   const bool &deep_copy=true) :
  Values(0), Nvector(0), Internal_values(deep_copy), Built(false)
  {
   //Build the storage based on the size of index
   unsigned n_vector = index.size();
   if(deep_copy)
    {
     //Create an entirely new data structure
     this->build(n_vector,old_vector.distribution_pt());
     //Now (deep) copy the data across
     const unsigned n_row_local = this->nrow_local();
     for(unsigned v=0;v<n_vector;v++)
      {
       for(unsigned i=0;i<n_row_local;i++)
        {
         Values[v][i] = old_vector(index[v],i);
        }
      }
    }
   //Otherwise it's a shallow copy
   else
    {
     this->shallow_build(n_vector,old_vector.distribution_pt());
     //Now shallow copy the pointers accross
     for(unsigned v=0;v<n_vector;++v)
      {
       Values[v] = old_vector.values(index[v]);
      }
    }
   this->setup_doublevector_representation();
  }

#ifdef OOMPH_HAS_TRILINOS
 /// \short Constructor that builds a multivector from selected columns
 /// of the input multivector and the provided range. The optional
 /// boolean specifies whether it is a shallow or deep copy. The default
 /// is that it is a deep copy.
 DoubleMultiVector(const DoubleMultiVector&old_vector,
                   const Teuchos::Range1D &index,
                   const bool &deep_copy=true) :
 Values(0), Nvector(0), Internal_values(deep_copy), Built(false)
  {
   //Build the storage based on the size of index
   unsigned n_vector = index.size();
   if(deep_copy)
    {
     //Create entirely new data structure
     this->build(n_vector,old_vector.distribution_pt());
     //Now (deep) copy the data across
     const unsigned n_row_local = this->nrow_local();
     unsigned range_index = index.lbound();
     for(unsigned v=0;v<n_vector;v++)
      {
       for(unsigned i=0;i<n_row_local;i++)
        {
         Values[v][i] = old_vector(range_index,i);
        }
       ++range_index;
      }
    }
   //Otherwise it's a shallow copy
   else
    {
     this->shallow_build(n_vector,old_vector.distribution_pt());
     //Shallow copy the pointers accross
     unsigned range_index = index.lbound();
     for(unsigned v=0;v<n_vector;v++)
      {
       Values[v] = old_vector.values(range_index);
       ++range_index;
      }
    }
   this->setup_doublevector_representation();
  }
#endif
 
 /// Copy constructor
 DoubleMultiVector(const DoubleMultiVector& new_vector)
  : DistributableLinearAlgebraObject(),
  Values(0), Nvector(0), Internal_values(true), Built(false)
  {
   this->build(new_vector);
   this->setup_doublevector_representation();
  }


 /// Destructor - just calls this->clear() to delete the distribution and data
 ~DoubleMultiVector()
  {
   this->clear();
  }                 
 
 /// assignment operator (deep copy)
 void operator=(const DoubleMultiVector& old_vector)
  {
   this->build(old_vector);
   this->setup_doublevector_representation();
  }

 /// Return the number of vectors
 unsigned nvector() const {return Nvector;}

 /// \short Provide a (shallow) copy of the old vector
 void shallow_build(const DoubleMultiVector& old_vector)
  {
   //Only bother if the old_vector is not the same as current vector
   if (!(*this == old_vector))
    {
     // the vector does not own the internal data
     Internal_values = false;
     
     //Copy the number of vectors
     Nvector = old_vector.nvector();
     //Allocate storage for pointers to the values
     this->shallow_build(Nvector,old_vector.distribution_pt());
     
     // copy all the pointers accross
     if (this->distribution_built())
      {
       for(unsigned v=0;v<Nvector;++v)
        {
         Values[v] = old_vector.values(v);
        }
      }
    }
  }

 /// \short Build the storage for pointers to vectors with a given
 /// distribution, but do not populate the pointers
 void shallow_build(const unsigned &n_vector,
                    const LinearAlgebraDistribution &dist)
  {
   this->shallow_build(n_vector,&dist);
  }


 /// \short Build the storage for pointers to vectors with a given
 /// distribution, but do not populate the pointers
 void shallow_build(const unsigned &n_vector,
                    const LinearAlgebraDistribution* const &dist_pt)
  {
   // clean the memory
   this->clear();

   //The vector does not own the data
   Internal_values = false;

   // Set the distribution
   this->build_distribution(dist_pt);

   // update the values
   if (dist_pt->built())
    {
     //Set the number of vectors
     Nvector = n_vector;
     //Build the array of pointers to each vector's data
     Values = new double*[n_vector];
     Built=true;
    }
   else
    {
     Built=false;
    }
  }

 
 /// \short Provides a (deep) copy of the old_vector
 void build(const DoubleMultiVector& old_vector)
  {
   //Only bother if the old_vector is not the same as current vector
   if (!(*this == old_vector))
    {
     // the vector owns the internal data
     Internal_values = true;
     
     //Copy the number of vectors
     Nvector = old_vector.nvector();
     // reset the distribution and resize the data
     this->build(Nvector,old_vector.distribution_pt(),0.0);
     
     // copy the data
     if (this->distribution_built())
      {
       unsigned n_row_local = this->nrow_local();
       double** const old_vector_values = old_vector.values();
       for (unsigned i = 0; i < n_row_local; i++)
        {
         for(unsigned v=0;v<Nvector;v++)
          {
           Values[v][i] = old_vector_values[v][i];
          }
        }
      }
    }
  }

 /// \short Assembles a DoubleMultiVector 
 /// with n_vector vectors, a distribution dist, if v is specified 
 /// each element is set to v, otherwise each element is set to 0.0
 void build(const unsigned &n_vector,
            const LinearAlgebraDistribution& dist, 
            const double& initial_value=0.0)
 {
  this->build(n_vector,&dist,initial_value);
 }

 /// \short 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
 void build(const unsigned &n_vector,
            const LinearAlgebraDistribution* const &dist_pt, 
            const double& initial_value=0.0)
  {
   // clean the memory
   this->clear();

   // the vector owns the internal data
   Internal_values = true;

   // Set the distribution
   this->build_distribution(dist_pt);

   // update the values
   if (dist_pt->built())
    {
     //Set the number of vectors
     Nvector = n_vector;
     //Build the array of pointers to each vector's data
     Values = new double*[n_vector];
     //Now build the contiguous array of real data
     const unsigned n_row_local = this->nrow_local();
     double *values = new double[n_row_local*n_vector];
     // set the data
     for(unsigned v=0;v<n_vector;v++)
      {
       Values[v] = &values[v*n_row_local];
       for (unsigned i = 0; i < n_row_local; i++)
        {
         Values[v][i] = initial_value;
        }
      }
     Built=true;
    }
   else
    {
     Built=false;
    }
  }

 /// \short initialise the whole vector with value v
 void initialise(const double& initial_value)
  {
   if (Built)
    {
     const unsigned n_vector = this->Nvector;
     const unsigned n_row_local = this->nrow_local();
     
     // set the residuals
     for(unsigned v=0;v<n_vector;v++)
      {
       for (unsigned i=0; i<n_row_local; i++)
        {
         Values[v][i] = initial_value;
        }
      }
    }
  }

 /// \short initialise the vector with coefficient from the vector v.
 /// Note: The vector v must be of length 
 //void initialise(const Vector<double> v);

 /// \short wipes the DoubleVector
 void clear() 
  {
   //Return if nothing to do
   if(Values==0) {return;}

   //If we are in charge of the data then delete it
   if (Internal_values)
    {
     //Delete the double storage arrays at once 
     //(they were allocated as a contiguous block)
     delete[] Values[0];
    }

   //Always Delete the pointers to the arrays
   delete[] Values;

   //Then set the pointer (to a pointer) to null
   Values = 0;
   this->clear_distribution();
   Built=false;
  }

 // indicates whether this DoubleVector is built
 bool built() const { return Built; }

 /// \short Allows are external data to be used by this vector. 
 /// WARNING: The size of the external data must correspond to the 
 /// LinearAlgebraDistribution dist_pt argument.
 /// 1. When a rebuild method is called new internal values are created.
 /// 2. It is not possible to redistribute(...) a vector with external
 /// values .
 /// 3. External values are only deleted by this vector if
 /// delete_external_values = true.
 /*void set_external_values(const LinearAlgebraDistribution* const& dist_pt,
                          double* external_values,
                          bool delete_external_values)
  {
   // clean the memory
   this->clear();

   // Set the distribution
   this->build_distribution(dist_pt);

   // set the external values
   set_external_values(external_values,delete_external_values);
   }*/
 
 /// \short Allows are external data to be used by this vector. 
 /// WARNING: The size of the external data must correspond to the 
 /// distribution of this vector.
 /// 1. When a rebuild method is called new internal values are created.
 /// 2. It is not possible to redistribute(...) a vector with external
 /// values .
 /// 3. External values are only deleted by this vector if
 /// delete_external_values = true.
 /*void set_external_values(double* external_values, 
                          bool delete_external_values)
  {
#ifdef PARANOID
   // check that this distribution is setup
   if (!this->distribution_built())
    {
    // if this vector does not own the double* values then it cannot be
    // distributed.
    // note: this is not stictly necessary - would just need to be careful 
    // with delete[] below.
     std::ostringstream error_message;    
     error_message << "The distribution of the vector must be setup before "
                   << "external values can be set"; 
     throw OomphLibError(error_message.str(),
     OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
   if (Internal_values)
    {
     delete[] Values_pt;
    }
   Values_pt = external_values;
   Internal_values = delete_external_values;
   }*/

 /// \short The contents of the vector are redistributed to match the new
 /// distribution. In a non-MPI rebuild this method works, but does nothing. 
 /// \b NOTE 1: The current distribution and the new distribution must have
 /// the same number of global rows.
 /// \b NOTE 2: The current distribution and the new distribution must have
 /// the same Communicator.
 void redistribute(const LinearAlgebraDistribution* const& dist_pt);
   
 /// \short [] access function to the (local) values of the v-th vector
 double& operator()(int v, int i) const
  {
#ifdef RANGE_CHECKING
   std::ostringstream error_message;
   bool error=false;
   if(v > int(Nvector))
    {
     error_message << "Range Error: Vector " << v 
                   << "is not in the range (0,"
                   << Nvector-1 << ")";
     error = true;
    }
   
   if (i >= int(this->nrow_local()))
    {   
     error_message << "Range Error: " << i 
                   << " is not in the range (0," 
                   << this->nrow_local()-1 << ")"; 
     error = true;
    }

   if(error)
    {
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
   return Values[v][i];
  }

 /// \short == operator
 bool operator==(const DoubleMultiVector& vec)
  {
   // if vec is not setup return false 
   if (vec.built() && !this->built())
    {
     return false;
    }
   else if (!vec.built() && this->built())
    {
     return false;
    }
   else if (!vec.built() && !this->built())
    {
     return true;
    }
   else
    {
     double** const v_values = vec.values();
     const unsigned n_row_local = this->nrow_local();
     const unsigned n_vector = this->nvector();
     for(unsigned v=0;v<n_vector;++v)
      {
       for(unsigned i=0;i<n_row_local;++i)
        {
         if (Values[v][i] != v_values[v][i])
          {
           return false;
          }
        }
      }
     return true;
    }
  }

 /// \short += operator
 void operator+=(DoubleMultiVector vec)
  {
#ifdef PARANOID
   // PARANOID check that this vector is setup
   if (!this->built())
    {
     std::ostringstream error_message;
     error_message << "This vector must be setup."; 
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
   // PARANOID check that the vector v is setup
   if (!vec.built())
    {
     std::ostringstream error_message;
     error_message << "The vector v must be setup."; 
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
   // PARANOID check that the vectors have the same distribution
   if (!(*vec.distribution_pt() == *this->distribution_pt()))
    {
     std::ostringstream error_message;
     error_message << "The vector v and this vector must have the same "
                   << "distribution."; 
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif//


   double** v_values = vec.values();
   const unsigned n_vector = this->nvector();
   const unsigned n_row_local = this->nrow_local();
   for(unsigned v=0; v<n_vector;++v)
    {
     for(unsigned i=0; i<n_row_local;++i)
      {
       Values[v][i] += v_values[v][i];
      }
    }
  }

 /// -= operator
 void operator-=(DoubleMultiVector vec)
  {
#ifdef PARANOID
   // PARANOID check that this vector is setup
   if (!this->distribution_built())
    {
     std::ostringstream error_message;
     error_message << "This vector must be setup."; 
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
   // PARANOID check that the vector v is setup
   if (!vec.built())
    {
     std::ostringstream error_message;
     error_message << "The vector v must be setup."; 
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
   // PARANOID check that the vectors have the same distribution
   if (!(*vec.distribution_pt() == *this->distribution_pt()))
    {
     std::ostringstream error_message;
     error_message << "The vector v and this vector must have the same "
                   << "distribution."; 
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif

   double** v_values = vec.values();
   const unsigned n_vector = this->nvector();
   const unsigned n_row_local = this->nrow_local();
   for(unsigned v=0;v<n_vector;++v)
    {
     for(unsigned i=0;i<n_row_local;++i)
      {
       Values[v][i] -= v_values[v][i];
      }
    }
  }

 ///Multiply by a scalar
 void operator *=(const double& scalar_value)
 {
  #ifdef PARANOID
   // PARANOID check that this vector is setup
   if (!this->distribution_built())
    {
     std::ostringstream error_message;
     error_message << "This vector must be setup."; 
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif
   const unsigned n_vector = this->nvector();
   const unsigned n_row_local = this->nrow_local();
   for(unsigned v=0;v<n_vector;++v)
    {
     for(unsigned i=0;i<n_row_local;++i)
      {
       Values[v][i] *= scalar_value;
      }
    }
 }
   

 /// access function to the underlying values
 double** values() {return Values;}

 /// \short access function to the underlying values (const version)
 double** values() const {return Values;}

 /// access function to the i-th vector's data
 double* values(const unsigned &i) {return Values[i];}
 
 /// access function to the i-th vector's data (const version)
 double* values(const unsigned &i) const {return Values[i];}

 /// access to the DoubleVector representatoin
 DoubleVector &doublevector(const unsigned &i) 
  {return Internal_doublevector[i];}

 /// access to the DoubleVector representation (const version)
 const DoubleVector &doublevector(const unsigned &i) const
  {return Internal_doublevector[i];}

 /// output the contents of the vector
 void output(std::ostream &outfile) const
  {
   // temp pointer to values
   double** temp;

   // Number of vectors
   unsigned n_vector = this->nvector();
 
   // number of global row
   unsigned nrow = this->nrow();
  
#ifdef OOMPH_HAS_MPI

   // number of local rows
   int nrow_local = this->nrow_local();

   // gather from all processors
   if (this->distributed() && 
       this->distribution_pt()->communicator_pt()->nproc() > 1)
    {
     // number of processors
     int nproc = this->distribution_pt()->communicator_pt()->nproc();

     // number of gobal row
     unsigned nrow = this->nrow();

     // get the vector of first_row s and nrow_local s
     int* dist_first_row = new int[nproc];
     int* dist_nrow_local =  new int[nproc];
     for (int p = 0; p < nproc; p++)
      {
       dist_first_row[p] = this->first_row(p);
       dist_nrow_local[p] = this->nrow_local(p);
      }

     // gather
     temp = new double*[n_vector];
     double* temp_value = new double[nrow*n_vector];
     for(unsigned v=0;v<n_vector;v++)
      {
       temp[v] = &temp_value[v*nrow];
      }

     //Now do an all gather for each vector
     //Possibly costly in terms of extra communication, but it's only output!
     for(unsigned v=0;v<n_vector;++v)
      {
       MPI_Allgatherv(Values[v],nrow_local,MPI_DOUBLE,
                      temp[v],dist_nrow_local,dist_first_row,
                      MPI_DOUBLE,
                      this->distribution_pt()->communicator_pt()->mpi_comm());
      }

     // clean up
     delete[] dist_first_row;
     delete[] dist_nrow_local;
    }
   else
    {
     temp = Values;
    }
#else
   temp = Values;
#endif

   // output
     for (unsigned i = 0; i < nrow; i++)
      {
       outfile << i << " ";
       for(unsigned v = 0; v<n_vector; v++)
        {
         outfile << temp[v][i] << " ";
        }
       outfile << "\n";
      }
     
     // clean up if required
#ifdef OOMPH_HAS_MPI
   if (this->distributed() && 
       this->distribution_pt()->communicator_pt()->nproc() > 1)
    {
     delete[] temp[0];
     delete[] temp;
    }
#endif

  }

 /// output the contents of the vector
 void output(std::string filename)
  {
    // Open file
    std::ofstream some_file;
    some_file.open(filename.c_str());
    output(some_file);
    some_file.close();
  }


 /// compute the 2 norm of this vector
 void dot(const DoubleMultiVector& vec, std::vector<double> &result) const
  {
#ifdef PARANOID
   // paranoid check that the vector is setup
    if (!this->built())
    {
     std::ostringstream error_message;
     error_message << "This vector must be setup."; 
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
   if (!vec.built())
    {
     std::ostringstream error_message;
     error_message << "The input vector be setup."; 
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
   if (*this->distribution_pt() != *vec.distribution_pt())
    {
     std::ostringstream error_message;    
     error_message << "The distribution of this vector and the vector vec "
                   << "must be the same."
                   << "\n\n  this: " << *this->distribution_pt()
                   << "\n  vec:  " << *vec.distribution_pt();
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif

   // compute the local norm
   unsigned nrow_local = this->nrow_local();
   int n_vector = this->nvector();
   double n[n_vector];
   for(int v=0;v<n_vector;v++)
    {
     //Initialise
     n[v] = 0.0;
     const double* vec_values_pt = vec.values(v);
     for (unsigned i = 0; i < nrow_local; i++)
      {
       n[v] += Values[v][i]*vec_values_pt[i];
      }
    }

   // if this vector is distributed and on multiple processors then gather
#ifdef OOMPH_HAS_MPI
   double n2[n_vector]; for(int v=0;v<n_vector;v++) {n2[v] = n[v];}

   if (this->distributed() && 
       this->distribution_pt()->communicator_pt()->nproc() > 1)
    {
     MPI_Allreduce(&n,&n2,n_vector,MPI_DOUBLE,MPI_SUM,
                   this->distribution_pt()->communicator_pt()->mpi_comm());
    }
   for(int v=0;v<n_vector;v++) {n[v] = n2[v];}
#endif

   result.resize(n_vector);
   for(int v=0;v<n_vector;v++) {result[v] = n[v];}
  }

 /// compute the 2 norm of this vector 
 void norm(std::vector<double> &result) const
  {
#ifdef PARANOID
   // paranoid check that the vector is setup
   if (!this->built())
    {
     std::ostringstream error_message;
     error_message << "This vector must be setup."; 
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif

   // compute the local norm
   unsigned nrow_local = this->nrow_local();
   int n_vector = this->nvector();
   double n[n_vector]; 
   for(int v=0;v<n_vector;v++)
    {
     n[v] = 0.0;
     for (unsigned i = 0; i < nrow_local; i++)
      {
       n[v] += Values[v][i]*Values[v][i];
      }
    }
   
   // if this vector is distributed and on multiple processors then gather
#ifdef OOMPH_HAS_MPI
   double n2[n_vector]; for(int v=0;v<n_vector;v++) {n2[v] = n[v];}
   if (this->distributed() && 
       this->distribution_pt()->communicator_pt()->nproc() > 1)
    {
     MPI_Allreduce(&n,&n2,n_vector,MPI_DOUBLE,MPI_SUM,
                   this->distribution_pt()->communicator_pt()->mpi_comm());
    }
   for(int v=0;v<n_vector;v++) {n[v] = n2[v];}
#endif

   //Now sqrt the norm and fill in result
   result.resize(n_vector);
   for(int v=0;v<n_vector;v++) {result[v] = sqrt(n[v]);}
  }

 /// compute the A-norm using the matrix at matrix_pt
 /*double norm(const CRDoubleMatrix* matrix_pt) const
  {
#ifdef PARANOID
   // paranoid check that the vector is setup
   if (!this->built())
    {
     std::ostringstream error_message;
     error_message << "This vector must be setup."; 
     throw OomphLibError(error_message.str(),
     OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
   if (!matrix_pt->built())
    {
     std::ostringstream error_message;
     error_message << "The input matrix be built."; 
     throw OomphLibError(error_message.str(),
     OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
   if (*this->distribution_pt() != *matrix_pt->distribution_pt())
    {
     std::ostringstream error_message;    
     error_message << "The distribution of this vector and the matrix at "
                   << "matrix_pt must be the same";
     throw OomphLibError(error_message.str(),
     OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif

   // compute the matrix norm
   DoubleVector x(this->distribution_pt(),0.0);
   matrix_pt->multiply(*this,x);
   return sqrt(this->dot(x));

   }*/
 
 private :

 /// Setup the doublevector representation
 void setup_doublevector_representation()
  {
   const unsigned n_vector = this->nvector();
   Internal_doublevector.resize(n_vector);
   //Loop over all and set the external values of the DoubleVectors
   for(unsigned v=0;v<n_vector;v++)
    {
     Internal_doublevector[v].set_external_values(this->distribution_pt(),
                                                  this->values(v),false);
    }
  }
 
 /// the local data, need a pointer to a pointer so that the
 /// individual vectors can be extracted
 double** Values;

 /// The number of vectors
 unsigned Nvector;

 /// \short Boolean flag to indicate whether the vector's data (values_pt) 
 /// is owned by this vector.
 bool Internal_values;

 /// indicates that the vector has been built and is usable
 bool Built;

 /// Need a vector of DoubleVectors to interface with our linear solvers
 Vector<DoubleVector> Internal_doublevector;

}; //end of DoubleVector                 

} // end of oomph namespace



#endif