]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add read/write operations to LinearAlgebra::Vector.
authorBruno Turcksin <bruno.turcksin@gmail.com>
Wed, 11 Nov 2015 14:39:40 +0000 (08:39 -0600)
committerBruno Turcksin <bruno.turcksin@gmail.com>
Wed, 11 Nov 2015 14:39:40 +0000 (08:39 -0600)
include/deal.II/base/index_set.h
include/deal.II/lac/la_vector.h
include/deal.II/lac/la_vector.templates.h
include/deal.II/lac/read_write_vector.templates.h
include/deal.II/lac/vector_space_vector.h
tests/lac/la_vector_add_and_dot.cc
tests/lac/la_vector_add_and_dot.output
tests/lac/la_vector_norms.cc
tests/lac/la_vector_output.cc [new file with mode: 0644]
tests/lac/la_vector_output.output [new file with mode: 0644]

index b0dc5a4016fa05a4ad2d81f99cf00eb9d3195eee..0ce6ef03b1034309fc846f29f7a7f0390e35a1e9 100644 (file)
@@ -19,6 +19,7 @@
 #include <deal.II/base/config.h>
 #include <deal.II/base/utilities.h>
 #include <deal.II/base/exceptions.h>
+#include <boost/serialization/vector.hpp>
 #include <vector>
 #include <algorithm>
 
index 0e0ea0acb7dde5fcee18a33f8d3f4af5ba67bafd..519fe6114ff0d411d91ed82d6b914fa5183e1957 100644 (file)
@@ -23,6 +23,8 @@
 #include <deal.II/base/index_set.h>
 #include <deal.II/lac/read_write_vector.h>
 #include <deal.II/lac/vector_space_vector.h>
+#include <boost/serialization/array.hpp>
+#include <boost/serialization/split_member.hpp>
 
 #include <cstdio>
 #include <iostream>
@@ -88,6 +90,23 @@ namespace LinearAlgebra
      */
     virtual ~Vector();
 
+    /**
+     * Copies the data of the input vector @p in_vector.
+     */
+    Vector<Number> &operator= (const Vector<Number> &in_vector);
+
+    /**
+     * Copies the data of the input vector @p in_vector.
+     */
+    template <typename Number2>
+    Vector<Number> &operator= (const Vector<Number2> &in_vector);
+
+    /**
+     * Sets all elements of the vector to the scalar @p s. This operation is
+     * only allowed if @p s is equal to zero.
+     */
+    Vector<Number> &operator= (const Number s);
+
     /**
      * Multiply the entire vector by a fixed factor.
      */
@@ -113,6 +132,11 @@ namespace LinearAlgebra
      */
     virtual Number operator* (const VectorSpaceVector<Number> &V);
 
+    /**
+     * Add @p a to all components. Note that @p a is a scalar not a vector.
+     */
+    virtual void add(const Number a);
+
     /**
      * Simple addition of a multiple of a vector, i.e. <tt>*this += a*V</tt>.
      */
@@ -204,6 +228,25 @@ namespace LinearAlgebra
                        const unsigned int precision=3,
                        const bool scientific=true,
                        const bool across=true) const;
+    /**
+     * Write the vector en bloc to a file. This is done in a binary mode, so the
+     * output is neither readable by humans nor (probably) by other computers
+     * using a different operating system or number format.
+     */
+    void block_write (std::ostream &out) const;
+
+    /**
+     * Read a vector en block from a file. This is done using the inverse
+     * operations to the above function, so it is reasonably fast because the
+     * bitstream is not interpreted.
+     *
+     * The vector is resized if necessary.
+     *
+     * A primitive form of error checking is performed which will recognize the
+     * bluntest attempts to interpret some data as a vector stored bitwise to a
+     * file, but not more.
+     */
+    void block_read (std::istream &in);
 
     /**
      * Returns the memory consumption of this class in bytes.
@@ -233,6 +276,15 @@ namespace LinearAlgebra
     typename VectorSpaceVector<Number>::real_type l2_norm_squared_recursive(
       unsigned int i,
       unsigned int j);
+
+    /**
+     * Serialize the data of this object using boost. This function is necessary
+     * to use boost::archive::text_iarchive and boost::archive::text_oarchive.
+     */
+    template <typename Archive>
+    void serialize(Archive &ar, const unsigned int version);
+
+    friend class boost::serialization::access;
   };
 
   /*@}*/
@@ -279,6 +331,56 @@ namespace LinearAlgebra
 
 
 
+  template <typename Number>
+  inline
+  Vector<Number> &Vector<Number>::operator= (const Vector<Number> &in_vector)
+  {
+    this->reinit(in_vector.size(),true);
+    std::copy(in_vector.begin(), in_vector.end(), this->begin());
+
+    return *this;
+  }
+
+
+
+  template <typename Number>
+  template <typename Number2>
+  inline
+  Vector<Number> &Vector<Number>::operator= (const Vector<Number2> &in_vector)
+  {
+    this->reinit(in_vector.size(in_vector.size()),true);
+    std::copy(in_vector.begin(), in_vector.end(), this->begin());
+
+    return *this;
+  }
+
+
+
+  template <typename Number>
+  inline
+  Vector<Number> &Vector<Number>::operator= (const Number s)
+  {
+    Assert(s==static_cast<Number>(0), ExcMessage("Only 0 can be assigned to a vector."));
+    (void) s;
+
+    std::fill(this->begin(),this->end(),Number());
+
+    return *this;
+  }
+
+
+
+  template <typename Number>
+  inline
+  void Vector<Number>::add(const Number a)
+  {
+    AssertIsFinite(a);
+    for (unsigned int i=0; i<this->size(); ++i)
+      this->val[i] += a;
+  }
+
+
+
   template <typename Number>
   inline
   typename Vector<Number>::size_type Vector<Number>::size() const
@@ -309,6 +411,22 @@ namespace LinearAlgebra
 
 
 
+  template <typename Number>
+  template <typename Archive>
+  inline
+  void Vector<Number>::serialize(Archive &ar, const unsigned int)
+  {
+    unsigned int current_size = this->size();
+    ar &static_cast<Subscriptor &>(*this);
+    ar &this->stored_elements;
+    // If necessary, resize the vector during a read operation
+    if (this->size()!=current_size)
+      this->reinit(this->size());
+    ar &boost::serialization::make_array(this->val,this->size());
+  }
+
+
+
   template <typename Number>
   inline
   std::size_t Vector<Number>::memory_consumption() const
index f436504d48dabd603906271c1fd4dac225d569b4..92f8d0bbe751cdeadca07f7f5510c14548e36662 100644 (file)
@@ -17,6 +17,8 @@
 #define dealii__la_vector_templates_h
 
 #include <deal.II/lac/la_vector.h>
+#include <iostream>
+#include <iomanip>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -346,6 +348,67 @@ namespace LinearAlgebra
 
     return norm;
   }
+
+
+
+  template <typename Number>
+  void Vector<Number>::block_write(std::ostream &out) const
+  {
+    AssertThrow(out, ExcIO());
+
+    // Other version of the following
+    //  out << size() << std::endl << '[';
+    // Reason: operator<< seems to use some resources  that lead to problems in
+    // a multithreaded environment.
+    const size_type sz = this->size();
+    char buf[16];
+#ifdef DEAL_II_WITH_64BIT_INDICES
+    std::sprintf(buf, "%llu", sz);
+#else
+    std::sprintf(buf, "%u", sz);
+#endif
+    std::strcat(buf, "\n[");
+
+    out.write(buf, std::strlen(buf));
+    out.write(reinterpret_cast<const char *>(this->begin()),
+              reinterpret_cast<const char *>(this->end())
+              - reinterpret_cast<const char *>(this->begin()));
+
+    // out << ']';
+    const char outro = ']';
+    out.write(&outro, 1);
+
+    AssertThrow(out, ExcIO());
+  }
+
+
+
+  template <typename Number>
+  void Vector<Number>::block_read(std::istream &in)
+  {
+    AssertThrow(in, ExcIO());
+
+    size_type sz;
+
+    char buf[16];
+
+    in.getline(buf,16,'\n');
+    sz = std::atoi(buf);
+    this->reinit(sz,true);
+
+    char c;
+    // in >> c;
+    in.read(&c,1);
+    AssertThrow(c=='[', ExcIO());
+
+    in.read(reinterpret_cast<char *>(this->begin()),
+            reinterpret_cast<const char *>(this->end())
+            - reinterpret_cast<const char *>(this->begin()));
+
+    // in >> c;
+    in.read(&c,1);
+    AssertThrow(c=']', ExcIO());
+  }
 }
 
 DEAL_II_NAMESPACE_CLOSE
index 35e2c301d6ddfff964389fd539fca7ef1a643605..9e334de71dd915cc45e04bcc1abd4ed96f304a11 100644 (file)
@@ -61,6 +61,7 @@ namespace LinearAlgebra
     resize_val(size);
 
     stored_elements = complete_index_set(size);
+    stored_elements.compress();
 
     // set entries to zero if so requested
     if (omit_zeroing_entries == false)
index 31ba1b1c7263f1b18793b79aff9f2983612a5d28..1f09437bac1937b8648bc59222190734db367d5c 100644 (file)
@@ -77,6 +77,11 @@ namespace LinearAlgebra
      */
     virtual Number operator* (const VectorSpaceVector<Number> &V) = 0;
 
+    /**
+     * Add @p a to all components. Note that @p a is a scalar not a vector.
+     */
+    virtual void add(const Number a) = 0;
+
     /**
      * Simple addition of a multiple of a vector, i.e. <tt>*this += a*V</tt>.
      */
index 4441163888dadda4c08025ebfc1d26eb8c19933b..07b68029d49849550bfb8e5e4d508b1708ae2412 100644 (file)
@@ -33,11 +33,13 @@ void check ()
     {
       const unsigned int size = 17 + test*1101;
       LinearAlgebra::Vector<number> v1 (size), v2(size), v3(size), check(size);
+      // Check that the assignment works
+      v1 = 0.;
       for (unsigned int i=0; i<size; ++i)
         {
-          v1(i) = 0.1 + 0.005 * i;
-          v2(i) = -5.2 + 0.18 * i;
-          v3(i) = 3.14159 + 2.7183/(1.+i);
+          v1[i] = 0.1 + 0.005 * i;
+          v2[i] = -5.2 + 0.18 * i;
+          v3[i] = 3.14159 + 2.7183/(1.+i);
         }
       check = v1;
       const number factor = 0.01432;
@@ -55,6 +57,13 @@ void check ()
           for (unsigned int i=0; i<size; ++i)
             deallog << check[i] << " ";
           deallog << std::endl;
+
+          const number constant = 1.;
+          v1.add(constant);
+          deallog << "Vector add constant:    ";
+          for (unsigned int i=0; i<size; ++i)
+            deallog << v1[i] << " ";
+          deallog << std::endl;
         }
 
       deallog << "Add and dot should be " << prod/static_cast<number>(size)
index 329e9f8d540501750a5c2189f5f576482908ef28..0c6f95da7351f19f08b4e9e469c8c8fb56e6b68f 100644 (file)
@@ -6,6 +6,7 @@ DEAL::Add and dot should be 39.61, is 39.61
 DEAL::Add and dot should be 52.71, is 52.71
 DEAL::Vector add reference:   0.03   0.03   0.04   0.05   0.06   0.06   0.07   0.08   0.09   0.09   0.10   0.11   0.12   0.12   0.13   0.14   0.15   
 DEAL::Vector check reference: 0.03   0.03   0.04   0.05   0.06   0.06   0.07   0.08   0.09   0.09   0.10   0.11   0.12   0.12   0.13   0.14   0.15   
+DEAL::Vector add constant:    1.03   1.03   1.04   1.05   1.06   1.06   1.07   1.08   1.09   1.09   1.10   1.11   1.12   1.12   1.13   1.14   1.15   
 DEAL::Add and dot should be 0.30, is 0.30
 DEAL::Add and dot should be 13.40, is 13.40
 DEAL::Add and dot should be 26.50, is 26.50
index 7f2faad3897838380fc7ea9a621b27dda38c2abb..e451006c2c8f7d519ee730c1a4fada70d2d25783 100644 (file)
@@ -40,7 +40,7 @@ void check_norms ()
         skip += 17;
       LinearAlgebra::Vector<number> vec(size);
       for (unsigned int i=0; i<size; ++i)
-        vec(i) = i+1;
+        vec[i] = i+1;
 
       const number l1_norm = vec.l1_norm();
       AssertThrow (std::abs(l1_norm-0.5*size*(size+1)) < acc*0.5*size*(size+1),
@@ -64,7 +64,7 @@ void check_norms ()
 template <typename number>
 void check_complex_norms ()
 {
-  const number acc = 1e1*std::numeric_limits<number>::epsilon();
+  const number acc = 1e2*std::numeric_limits<number>::epsilon();
   unsigned int skip = 73;
   for (unsigned int size=1; size<100000; size+=skip)
     {
diff --git a/tests/lac/la_vector_output.cc b/tests/lac/la_vector_output.cc
new file mode 100644 (file)
index 0000000..0aeb7d8
--- /dev/null
@@ -0,0 +1,88 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2012 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/la_vector.h>
+#include <fstream>
+#include <string>
+#include <boost/archive/text_oarchive.hpp>
+#include <boost/archive/text_iarchive.hpp>
+
+
+void test()
+{
+  std::string filename("test.txt");
+  const unsigned int size(10);
+  LinearAlgebra::Vector<double> vec(size);
+  for (unsigned int i=0; i<size; ++i)
+    vec[i] = i;
+
+  // Check block_write and block_read
+  std::ofstream file_out(filename);
+  // Write the vector in the file
+  vec.block_write(file_out);
+  file_out.close();
+  // Clear the vector
+  vec.reinit(0);
+  // Load the vector form the file
+  std::ifstream file_in(filename);
+  vec.block_read(file_in);
+  file_in.close();
+  // Check the vector
+  double eps=1e-12;
+  for (unsigned int i=0; i<size; ++i)
+    AssertThrow(std::abs(vec[i]-(double)i)<eps,
+                ExcMessage("Value in the vector has been changed by block_write or block_read"));
+
+
+  // save data to archive
+  {
+    std::ofstream file_out2(filename);
+    boost::archive::text_oarchive oa(file_out2);
+    oa<<vec;
+    // archive and stream closed when destructors are called
+  }
+
+  // Clear the vector
+  vec.reinit(0);
+  {
+    std::ifstream file_in2(filename);
+    boost::archive::text_iarchive ia(file_in2);
+    ia >> vec;
+  }
+  // Check the vector
+  for (unsigned int i=0; i<size; ++i)
+    AssertThrow(std::abs(vec[i]-(double)i)<eps,
+                ExcMessage("Value in the vector has been changed by boost archive"));
+}
+
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog << std::fixed;
+  deallog << std::setprecision(2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  test();
+
+  deallog << "OK" << std::endl;
+}
+
+
diff --git a/tests/lac/la_vector_output.output b/tests/lac/la_vector_output.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.