Changeset 4141

Timestamp:

Feb 1, 2022, 6:11:36 AM (20 months ago)

Author:

Peter

Message:

closes #977. Implement view into gsl_vector, gsl_matrix, and arrays.

Location:

trunk

Files:

• test/matrix_view.cc (modified) (2 diffs)
• test/view.cc (modified) (2 diffs)
• yat/utility/Matrix.h (modified) (1 diff)
• yat/utility/MatrixConstView.cc (modified) (1 diff)
• yat/utility/MatrixConstView.h (modified) (1 diff)
• yat/utility/MatrixView.cc (modified) (1 diff)
• yat/utility/MatrixView.h (modified) (1 diff)
• yat/utility/VectorConstView.cc (modified) (2 diffs)
• yat/utility/VectorConstView.h (modified) (1 diff)
• yat/utility/VectorView.cc (modified) (1 diff)
• yat/utility/VectorView.h (modified) (1 diff)

trunk/test/matrix_view.cc

@@ -219 +219 @@
   test_iterators<MatrixView>(suite, m);
   MatrixView view = make_view(m);
+
 }
 
@@ -242 +243 @@
 
   MatrixConstView view = make_const_view(m);
-}
+
+  std::vector<double> vec(100,0);
+  for (size_t i=0; i<vec.size(); ++i)
+    vec[i] = i;
+  utility::Matrix correct(4, 5);
+  for (size_t i=0; i<correct.rows(); ++i)
+    for (size_t j=0; j<correct.columns(); ++j)
+      correct(i,j) = 5*i + j;
+  MatrixConstView view2(&vec.front(), 4, 5);
+  if (!suite.equal_matrix(view2, correct)) {
+    suite.add(false);
+    suite.err() << "error: view2 not equal correct\n";
+    suite.err() << "view2:\n" << view2 << "\ncorrect:\n" << correct << "\n";
+  }
+
+  size_t tda = 10;
+  MatrixConstView view3(&vec.front(), 4, 5, tda);
+  for (size_t i=0; i<correct.rows(); ++i)
+    for (size_t j=0; j<correct.columns(); ++j)
+      correct(i,j) = tda*i + j;
+  if (!suite.equal_matrix(view3, correct)) {
+    suite.add(false);
+    suite.err() << "error: view3 not equal correct\n";
+    suite.err() << "view3:\n" << view3 << "\ncorrect:\n" << correct << "\n";
+  }
+
+  MatrixConstView view4(correct.gsl_matrix_p());
+  if (!suite.equal_matrix(view4, correct)) {
+    suite.add(false);
+    suite.err() << "error: view4 not equal correct\n";
+    suite.err() << "view4:\n" << view4 << "\ncorrect:\n" << correct << "\n";
+  }
+}

trunk/test/view.cc

@@ -40 +40 @@
 
 int main(int argc, char* argv[])
-{  
+{
   using namespace theplu::yat::utility;
   theplu::yat::test::Suite suite(argc, argv);
@@ -132 +132 @@
       suite.add(false);
   }
-  
   test_rvalue(suite);
+
+  // test view from double* or gsl_vector*
+  std::vector<double> vec(100,0);
+  VectorView view2(&vec.front(), 5);
+  VectorConstView cview2(&vec.front(), 5);
+  if (view2.size() != 5) {
+    suite.err() << "error: view2 incorrect dimensions\n";
+    suite.add(false);
+  }
+  else {
+    for (size_t i=0; i<view2.size(); ++i) {
+      double x = 12*i;
+      view2(i) = x;
+      if (vec[i] != x) {
+        suite.add(false);
+      }
+    }
+  }
+  if (view2 != cview2) {
+    suite.err() << "error: view2 != cview2\n";
+    suite.add(false);
+  }
+
+  size_t stride = 3;
+  VectorView view3(&vec.front(), 5, stride);
+  VectorConstView cview3(&vec.front(), 5, stride);
+  if (view3.size() != 5) {
+    suite.err() << "error: view3 incorrect dimensions\n";
+    suite.add(false);
+  }
+  else {
+    for (size_t i=0; i<view3.size(); ++i) {
+      double x = 12*i;
+      view3(i) = x;
+      if (vec[stride*i] != x) {
+        suite.add(false);
+      }
+    }
+  }
+  if (view3 != cview3) {
+    suite.err() << "error: view3 != cview3\n";
+    suite.add(false);
+  }
+
+  Vector v(5);
+  gsl_vector* p = v.gsl_vector_p();
+  const gsl_vector* cp = v.gsl_vector_p();
+  VectorView view4(p);
+  VectorConstView cview4(cp);
+  view4(0) = 3.12;
+  view4(3) = 4.28;
+  if (view4 != v) {
+    suite.add(false);
+    suite.err() << "error: v not equal view4\n";
+  }
+  if (view4 != cview4) {
+    suite.err() << "error: view4 != cview4\n";
+    suite.add(false);
+  }
 
   return suite.return_value();

trunk/yat/utility/Matrix.h

@@ -100 +100 @@
   /// most functionality defined for GSL matrices can be achieved with
   /// this interface class. Most notable GSL functionality not
-  /// supported are no binary file support and views on arrays,
-  /// utility::vectors, gsl_vectors, diagonals, subdiagonals, and
-  /// superdiagonals.
+  /// supported are no binary file support.
   ///
   class Matrix : public MatrixMutable

trunk/yat/utility/MatrixConstView.cc

@@ -59 +59 @@
 
 
+  MatrixConstView::MatrixConstView(const gsl_matrix* m)
+    : view_(nullptr)
+  {
+    copy(m);
+  }
+
+
+  MatrixConstView::MatrixConstView(const double* p, size_t rows, size_t cols)
+    : view_(nullptr)
+  {
+    assert(p);
+    assert(rows);
+    assert(cols);
+    view_ =
+      new gsl_matrix_const_view(gsl_matrix_const_view_array(p, rows, cols));
+  }
+
+
+  MatrixConstView::MatrixConstView(const double* p, size_t rows,
+                                   size_t cols, size_t tda)
+    : view_(nullptr)
+  {
+    assert(p);
+    assert(rows);
+    assert(cols);
+    assert(tda >= cols);
+    view_ =
+      new gsl_matrix_const_view(gsl_matrix_const_view_array_with_tda(p, rows,
+                                                                     cols,
+                                                                     tda));
+  }
+
+
   MatrixConstView::~MatrixConstView(void)
   {

trunk/yat/utility/MatrixConstView.h

@@ -69 +69 @@
 
     /**
+       Create a view into \c m.
+     */
+    explicit MatrixConstView(const gsl_matrix* m);
+
+    /**
+       Create a view with \c rows rows and \c cols columns. Data are
+       stored in row-major order such that p[1] correspond to
+       view(0,1).
+     */
+    MatrixConstView(const double* p, size_t rows, size_t cols);
+
+    /**
+       Create a view with \c rows rows and \c cols columns. The
+       physical row dimension is \c tda, i.e., the first element on
+       the second row is p[tda].
+     */
+    MatrixConstView(const double* p, size_t rows, size_t cols, size_t tda);
+
+    /**
        \brief Destructor
      */

trunk/yat/utility/MatrixView.cc

@@ -57 +57 @@
   {
     copy(other.gsl_matrix_p(), row_offset, col_offset, rows, columns);
+  }
+
+
+  MatrixView::MatrixView(gsl_matrix* m)
+  {
+    assert(0);
+  }
+
+
+  MatrixView::MatrixView(double* p, size_t rows, size_t cols)
+  {
+    assert(0);
+  }
+
+
+  MatrixView::MatrixView(double* p, size_t rows, size_t cols, size_t tda)
+  {
+    assert(0);
   }
 

trunk/yat/utility/MatrixView.h

@@ -85 +85 @@
 
     /**
+       Create a view into \c m.
+     */
+    explicit MatrixView(gsl_matrix* m);
+
+    /**
+       Create a view with \c rows rows and \c cols columns. Data are
+       stored in row-major order such that p[1] correspond to
+       view(0,1).
+     */
+    MatrixView(double* p, size_t rows, size_t cols);
+
+    /**
+       Create a view with \c rows rows and \c cols columns. The
+       physical row dimension is \c tda, i.e., the first element on
+       the second row is p[tda].
+     */
+    MatrixView(double* p, size_t rows, size_t cols, size_t tda);
+
+    /**
        \brief Destructor
      */

trunk/yat/utility/VectorConstView.cc

@@ -29 +29 @@
 #include "MatrixBase.h"
 #include "VectorView.h"
+
+#include <cassert>
 
 namespace theplu {
@@ -78 +80 @@
 
 
+  VectorConstView::VectorConstView(const gsl_vector* v)
+  {
+    const_view_ =
+      new gsl_vector_const_view(gsl_vector_const_subvector(v, 0, v->size));
+    const_vec_ = &(const_view_->vector);
+  }
+
+
+  VectorConstView::VectorConstView(const double* v, size_t size, size_t stride)
+  {
+    const_view_ =
+      new gsl_vector_const_view(gsl_vector_const_view_array_with_stride(v,
+                                                                        stride,
+                                                                        size));
+    const_vec_ = &(const_view_->vector);
+  }
+
+
   VectorConstView::~VectorConstView(void)
   {

trunk/yat/utility/VectorConstView.h

@@ -113 +113 @@
 
     /**
+       Create a view into \c v
+
+       \since New in yat 0.20
+     */
+    explicit VectorConstView(const gsl_vector* v);
+
+    /**
+       Create a view with size \c size and stride \c stride such that
+       view(0) = v[0] and view(1) = v[stride].
+
+       \since New in yat 0.20
+     */
+    VectorConstView(const double* v, size_t size, size_t stride=1);
+
+    /**
        The destructor.
     */

trunk/yat/utility/VectorView.cc

@@ -83 +83 @@
   {
     copy(p.vec_);
+  }
+
+
+  VectorView::VectorView(gsl_vector* v)
+  {
+    copy(v);
+  }
+
+
+  VectorView::VectorView(double* v, size_t size, size_t stride)
+  {
+    assert(v);
+    assert(size);
+    assert(stride);
+    view_ =
+      new gsl_vector_view(gsl_vector_view_array_with_stride(v, stride, size));
+    const_vec_ = vec_ = &(view_->vector);
   }
 

trunk/yat/utility/VectorView.h

@@ -137 +137 @@
     VectorView(MatrixMutable& m, size_t i, bool row=true);
 
+    /**
+       Create a view into \c v
+
+       \since New in yat 0.20
+     */
+    explicit VectorView(gsl_vector* v);
+
+    /**
+       Create a view with size \c size and stride \c stride such that
+       view(0) = v[0] and view(1) = v[stride].
+
+       \since New in yat 0.20
+     */
+    VectorView(double* v, size_t size, size_t stride=1);
+
     ///
     /// The destructor.