]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix bug #78: BlockMatrixBase::set/add had race conditions when used from multiple...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 6 Aug 2013 21:44:57 +0000 (21:44 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 6 Aug 2013 21:44:57 +0000 (21:44 +0000)
git-svn-id: https://svn.dealii.org/trunk@30238 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/lac/block_matrix_base.h
tests/lac/block_matrices_02.cc [new file with mode: 0644]
tests/lac/block_matrices_02/cmp/generic [new file with mode: 0644]
tests/lac/block_matrices_03.cc [new file with mode: 0644]
tests/lac/block_matrices_03/cmp/generic [new file with mode: 0644]

index be39ececbd3a0024ac17ec160c13f4e485b8f26b..e87472eac821aaf5b760856ae4671c140e3beea5 100644 (file)
@@ -44,6 +44,14 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li>
+  Fixed: The various block matrix classes are all derived from BlockMatrixBase
+  which had race conditions when the set() or add() functions were called from
+  different threads. This is now fixed.
+  <br>
+  (Wolfgang Bangerth, 2013/08/05)
+  </li>
+
   <li>
   Fixed: various fixes with assignment and reinit of PETScWrappers::MPI::Vector.
   <br>
index 273ad266bead7b78e6e4c1c8d053881b8fefa7d0..6ed20d407ef63845a6723c98792b8da07a9a43bc 100644 (file)
@@ -20,6 +20,7 @@
 
 #include <deal.II/base/config.h>
 #include <deal.II/base/table.h>
+#include <deal.II/base/thread_management.h>
 #include <deal.II/base/utilities.h>
 #include <deal.II/base/smartpointer.h>
 #include <deal.II/base/memory_consumption.h>
@@ -1246,30 +1247,72 @@ protected:
 
 
 private:
-  /**
-   * Temporary vector for counting the
-   * elements written into the
-   * individual blocks when doing a
-   * collective add or set.
-   */
-  std::vector<size_type> counter_within_block;
 
   /**
-   * Temporary vector for column
-   * indices on each block when writing
-   * local to global data on each
-   * sparse matrix.
+   * A structure containing some fields used by the
+   * set() and add() functions that is used to pre-sort
+   * the input fields. Since one can reasonably expect
+   * to call set() and add() from multiple threads at once
+   * as long as the matrix indices that are touched are
+   * disjoint, these temporary data fields need to be
+   * guarded by a mutex; the structure therefore contains such
+   * a mutex as a member variable.
    */
-  std::vector<std::vector<size_type> > column_indices;
+  struct TemporaryData
+  {
+    /**
+     * Temporary vector for counting the
+     * elements written into the
+     * individual blocks when doing a
+     * collective add or set.
+     */
+    std::vector<size_type> counter_within_block;
+
+    /**
+     * Temporary vector for column
+     * indices on each block when writing
+     * local to global data on each
+     * sparse matrix.
+     */
+    std::vector<std::vector<size_type> > column_indices;
+
+    /**
+     * Temporary vector for storing the
+     * local values (they need to be
+     * reordered when writing local to
+     * global).
+     */
+    std::vector<std::vector<double> > column_values;
+
+    /**
+     * A mutex variable used to guard access to the member
+     * variables of this structure;
+     */
+    Threads::Mutex mutex;
+
+    /**
+     * Copy operator. This is needed because the default copy
+     * operator of this class is deleted (since Threads::Mutex is
+     * not copyable) and hence the default copy operator of the
+     * enclosing class is also deleted.
+     *
+     * The implementation here simply does nothing -- TemporaryData
+     * objects are just scratch objects that are resized at the
+     * beginning of their use, so there is no point actually copying
+     * anything.
+     */
+    TemporaryData & operator = (const TemporaryData &)
+    {}
+  };
 
   /**
-   * Temporary vector for storing the
-   * local values (they need to be
-   * reordered when writing local to
-   * global).
+   * A set of scratch arrays that can be used by the add()
+   * and set() functions that take pointers to data to
+   * pre-sort indices before use. Access from multiple threads
+   * is synchronized via the mutex variable that is part of the
+   * structure.
    */
-  std::vector<std::vector<double> > column_values;
-
+  TemporaryData temporary_data;
 
   /**
    * Make the iterator class a
@@ -1772,9 +1815,10 @@ BlockMatrixBase<MatrixType>::memory_consumption () const
     MemoryConsumption::memory_consumption(row_block_indices)+
     MemoryConsumption::memory_consumption(column_block_indices)+
     MemoryConsumption::memory_consumption(sub_objects)+
-    MemoryConsumption::memory_consumption(counter_within_block)+
-    MemoryConsumption::memory_consumption(column_indices)+
-    MemoryConsumption::memory_consumption(column_values);
+    MemoryConsumption::memory_consumption(temporary_data.counter_within_block)+
+    MemoryConsumption::memory_consumption(temporary_data.column_indices)+
+    MemoryConsumption::memory_consumption(temporary_data.column_values)+
+    MemoryConsumption::memory_consumption(temporary_data.mutex);
 
   for (unsigned int r=0; r<n_block_rows(); ++r)
     for (unsigned int c=0; c<n_block_cols(); ++c)
@@ -1976,12 +2020,16 @@ BlockMatrixBase<MatrixType>::set (const size_type  row,
 {
   prepare_set_operation();
 
+  // lock access to the temporary data structure to
+  // allow multiple threads to call this function concurrently
+  Threads::Mutex::ScopedLock lock (temporary_data.mutex);
+
   // Resize scratch arrays
-  if (column_indices.size() < this->n_block_cols())
+  if (temporary_data.column_indices.size() < this->n_block_cols())
     {
-      column_indices.resize (this->n_block_cols());
-      column_values.resize (this->n_block_cols());
-      counter_within_block.resize (this->n_block_cols());
+      temporary_data.column_indices.resize (this->n_block_cols());
+      temporary_data.column_values.resize (this->n_block_cols());
+      temporary_data.counter_within_block.resize (this->n_block_cols());
     }
 
   // Resize sub-arrays to n_cols. This
@@ -1994,19 +2042,19 @@ BlockMatrixBase<MatrixType>::set (const size_type  row,
   // whether the size of one is large
   // enough before actually going
   // through all of them.
-  if (column_indices[0].size() < n_cols)
+  if (temporary_data.column_indices[0].size() < n_cols)
     {
       for (unsigned int i=0; i<this->n_block_cols(); ++i)
         {
-          column_indices[i].resize(n_cols);
-          column_values[i].resize(n_cols);
+          temporary_data.column_indices[i].resize(n_cols);
+          temporary_data.column_values[i].resize(n_cols);
         }
     }
 
   // Reset the number of added elements
   // in each block to zero.
   for (unsigned int i=0; i<this->n_block_cols(); ++i)
-    counter_within_block[i] = 0;
+    temporary_data.counter_within_block[i] = 0;
 
   // Go through the column indices to
   // find out which portions of the
@@ -2027,10 +2075,10 @@ BlockMatrixBase<MatrixType>::set (const size_type  row,
       const std::pair<unsigned int, size_type>
       col_index = this->column_block_indices.global_to_local(col_indices[j]);
 
-      const size_type local_index = counter_within_block[col_index.first]++;
+      const size_type local_index = temporary_data.counter_within_block[col_index.first]++;
 
-      column_indices[col_index.first][local_index] = col_index.second;
-      column_values[col_index.first][local_index] = value;
+      temporary_data.column_indices[col_index.first][local_index] = col_index.second;
+      temporary_data.column_values[col_index.first][local_index] = value;
     }
 
 #ifdef DEBUG
@@ -2038,7 +2086,7 @@ BlockMatrixBase<MatrixType>::set (const size_type  row,
   // the right length has been obtained.
   size_type length = 0;
   for (unsigned int i=0; i<this->n_block_cols(); ++i)
-    length += counter_within_block[i];
+    length += temporary_data.counter_within_block[i];
   Assert (length <= n_cols, ExcInternalError());
 #endif
 
@@ -2051,14 +2099,14 @@ BlockMatrixBase<MatrixType>::set (const size_type  row,
   row_index = this->row_block_indices.global_to_local (row);
   for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
     {
-      if (counter_within_block[block_col] == 0)
+      if (temporary_data.counter_within_block[block_col] == 0)
         continue;
 
       block(row_index.first, block_col).set
       (row_index.second,
-       counter_within_block[block_col],
-       &column_indices[block_col][0],
-       &column_values[block_col][0],
+       temporary_data.counter_within_block[block_col],
+       &temporary_data.column_indices[block_col][0],
+       &temporary_data.column_values[block_col][0],
        false);
     }
 }
@@ -2221,12 +2269,14 @@ BlockMatrixBase<MatrixType>::add (const size_type  row,
       return;
     }
 
-  // Resize scratch arrays
-  if (column_indices.size() < this->n_block_cols())
+  // Lock scratch arrays, then resize them
+  Threads::Mutex::ScopedLock lock (temporary_data.mutex);
+
+  if (temporary_data.column_indices.size() < this->n_block_cols())
     {
-      column_indices.resize (this->n_block_cols());
-      column_values.resize (this->n_block_cols());
-      counter_within_block.resize (this->n_block_cols());
+      temporary_data.column_indices.resize (this->n_block_cols());
+      temporary_data.column_values.resize (this->n_block_cols());
+      temporary_data.counter_within_block.resize (this->n_block_cols());
     }
 
   // Resize sub-arrays to n_cols. This
@@ -2239,19 +2289,19 @@ BlockMatrixBase<MatrixType>::add (const size_type  row,
   // whether the size of one is large
   // enough before actually going
   // through all of them.
-  if (column_indices[0].size() < n_cols)
+  if (temporary_data.column_indices[0].size() < n_cols)
     {
       for (unsigned int i=0; i<this->n_block_cols(); ++i)
         {
-          column_indices[i].resize(n_cols);
-          column_values[i].resize(n_cols);
+          temporary_data.column_indices[i].resize(n_cols);
+          temporary_data.column_values[i].resize(n_cols);
         }
     }
 
   // Reset the number of added elements
   // in each block to zero.
   for (unsigned int i=0; i<this->n_block_cols(); ++i)
-    counter_within_block[i] = 0;
+    temporary_data.counter_within_block[i] = 0;
 
   // Go through the column indices to
   // find out which portions of the
@@ -2272,10 +2322,10 @@ BlockMatrixBase<MatrixType>::add (const size_type  row,
       const std::pair<unsigned int, size_type>
       col_index = this->column_block_indices.global_to_local(col_indices[j]);
 
-      const size_type local_index = counter_within_block[col_index.first]++;
+      const size_type local_index = temporary_data.counter_within_block[col_index.first]++;
 
-      column_indices[col_index.first][local_index] = col_index.second;
-      column_values[col_index.first][local_index] = value;
+      temporary_data.column_indices[col_index.first][local_index] = col_index.second;
+      temporary_data.column_values[col_index.first][local_index] = value;
     }
 
 #ifdef DEBUG
@@ -2283,7 +2333,7 @@ BlockMatrixBase<MatrixType>::add (const size_type  row,
   // the right length has been obtained.
   size_type length = 0;
   for (unsigned int i=0; i<this->n_block_cols(); ++i)
-    length += counter_within_block[i];
+    length += temporary_data.counter_within_block[i];
   Assert (length <= n_cols, ExcInternalError());
 #endif
 
@@ -2296,14 +2346,14 @@ BlockMatrixBase<MatrixType>::add (const size_type  row,
   row_index = this->row_block_indices.global_to_local (row);
   for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
     {
-      if (counter_within_block[block_col] == 0)
+      if (temporary_data.counter_within_block[block_col] == 0)
         continue;
 
       block(row_index.first, block_col).add
       (row_index.second,
-       counter_within_block[block_col],
-       &column_indices[block_col][0],
-       &column_values[block_col][0],
+       temporary_data.counter_within_block[block_col],
+       &temporary_data.column_indices[block_col][0],
+       &temporary_data.column_values[block_col][0],
        false,
        col_indices_are_sorted);
     }
diff --git a/tests/lac/block_matrices_02.cc b/tests/lac/block_matrices_02.cc
new file mode 100644 (file)
index 0000000..0eece67
--- /dev/null
@@ -0,0 +1,131 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2000 - 2013 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.
+//
+// ---------------------------------------------------------------------
+
+
+// the versions of BlockMatrixBase::set/add that take a pointer to a set of
+// values were not thread safe, even if they were writing into disjoint sets
+// of elements. this test verifies that the set() function is now indeed
+// thread safe under these conditions
+//
+// we test this by calling add() from a multitude of threads at once. without
+// the patch that fixed this, one can elicit a great deal of different memory
+// corruption errors and assertions from this test, depending on luck and
+// phase of the moon :-)
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/block_sparsity_pattern.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_vector.h>
+#include <fstream>
+#include <iomanip>
+#include <algorithm>
+
+
+
+void do_set (const bool even_or_odd,
+            BlockSparseMatrix<double> &bsm)
+{
+  BlockSparseMatrix<double>::size_type col_indices[5];
+  for (unsigned int i=0; i<5 ; ++i)
+    if (even_or_odd)
+      col_indices[i] = 2*i;
+    else
+      col_indices[i] = 2*i+1;
+
+  BlockSparseMatrix<double>::value_type values[5];
+  for (unsigned int i=0; i<5 ; ++i)
+    if (even_or_odd)
+      values[i] = 1;
+    else
+      values[i] = 2;
+
+  bsm.set (0, 5, col_indices, values, false);
+}
+
+
+void test ()
+{
+  std::ofstream logfile("block_matrices_02/output");
+  deallog << std::fixed;
+  deallog << std::setprecision(2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  BlockSparsityPattern bsp(2,2);
+  // set sizes
+  for (unsigned int i=0; i<2; ++i)
+    for (unsigned int j=0; j<2; ++j)
+      bsp.block(i,j).reinit ( 5, 5, 5);
+  bsp.collect_sizes ();
+
+  // make a full matrix
+  for (unsigned int row=0; row<10; ++row)
+    for (unsigned int i=0; i<10; ++i)
+      bsp.add (row, i);
+  bsp.compress ();
+
+  BlockSparseMatrix<double> bsm (bsp);
+
+  Threads::ThreadGroup<> tg;
+  for (unsigned int i=0; i<100; ++i)
+    {
+      tg += Threads::new_thread (&do_set, true, bsm);
+      tg += Threads::new_thread (&do_set, false, bsm);
+    }
+  tg.join_all ();
+
+  bsm.print_formatted (deallog.get_file_stream());
+}
+
+
+
+
+int main ()
+{
+  try
+    {
+      test ();
+    }
+  catch (std::exception &e)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << e.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      // abort
+      return 2;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      // abort
+      return 3;
+    };
+
+
+  return 0;
+}
diff --git a/tests/lac/block_matrices_02/cmp/generic b/tests/lac/block_matrices_02/cmp/generic
new file mode 100644 (file)
index 0000000..3dca993
--- /dev/null
@@ -0,0 +1,25 @@
+
+Component (0,0)
+1.000e+00  2.000e+00  1.000e+00  2.000e+00  1.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+Component (0,1)
+2.000e+00  1.000e+00  2.000e+00  1.000e+00  2.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+Component (1,0)
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+Component (1,1)
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
diff --git a/tests/lac/block_matrices_03.cc b/tests/lac/block_matrices_03.cc
new file mode 100644 (file)
index 0000000..16d9767
--- /dev/null
@@ -0,0 +1,134 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2000 - 2013 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.
+//
+// ---------------------------------------------------------------------
+
+
+// the versions of BlockMatrixBase::set/add that take a pointer to a set of
+// values were not thread safe, even if they were writing into disjoint sets
+// of elements. this test verifies that the add() function is now indeed
+// thread safe under these conditions
+//
+// we test this by calling add() from a multitude of threads at once. without
+// the patch that fixed this, one can elicit a great deal of different memory
+// corruption errors and assertions from this test, depending on luck and
+// phase of the moon :-)
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/block_sparsity_pattern.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_vector.h>
+#include <fstream>
+#include <iomanip>
+#include <algorithm>
+
+
+
+void do_add (const bool even_or_odd,
+            BlockSparseMatrix<double> &bsm)
+{
+  BlockSparseMatrix<double>::size_type col_indices[5];
+  for (unsigned int i=0; i<5 ; ++i)
+    if (even_or_odd)
+      col_indices[i] = 2*i;
+    else
+      col_indices[i] = 2*i+1;
+
+  BlockSparseMatrix<double>::value_type values[5];
+  for (unsigned int i=0; i<5 ; ++i)
+    if (even_or_odd)
+      values[i] = 1;
+    else
+      values[i] = 2;
+
+  bsm.add (0, 5, col_indices, values, false);
+}
+
+
+void test ()
+{
+  std::ofstream logfile("block_matrices_03/output");
+  deallog << std::fixed;
+  deallog << std::setprecision(2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  BlockSparsityPattern bsp(2,2);
+  // set sizes
+  for (unsigned int i=0; i<2; ++i)
+    for (unsigned int j=0; j<2; ++j)
+      bsp.block(i,j).reinit ( 5, 5, 5);
+  bsp.collect_sizes ();
+
+  // make a full matrix
+  for (unsigned int row=0; row<10; ++row)
+    for (unsigned int i=0; i<10; ++i)
+      bsp.add (row, i);
+  bsp.compress ();
+
+  BlockSparseMatrix<double> bsm (bsp);
+
+  Threads::ThreadGroup<> tg;
+  for (unsigned int i=0; i<100; ++i)
+    {
+      tg += Threads::new_thread (&do_add, true, bsm);
+      tg += Threads::new_thread (&do_add, false, bsm);
+    }
+  tg.join_all ();
+
+  // divide whole matrix by 100 to get back to the effect of a single set
+  bsm /= 100;
+  
+  bsm.print_formatted (deallog.get_file_stream());
+}
+
+
+
+
+int main ()
+{
+  try
+    {
+      test ();
+    }
+  catch (std::exception &e)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << e.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      // abort
+      return 2;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      // abort
+      return 3;
+    };
+
+
+  return 0;
+}
diff --git a/tests/lac/block_matrices_03/cmp/generic b/tests/lac/block_matrices_03/cmp/generic
new file mode 100644 (file)
index 0000000..3dca993
--- /dev/null
@@ -0,0 +1,25 @@
+
+Component (0,0)
+1.000e+00  2.000e+00  1.000e+00  2.000e+00  1.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+Component (0,1)
+2.000e+00  1.000e+00  2.000e+00  1.000e+00  2.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+Component (1,0)
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+Component (1,1)
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  
+0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  

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.