Changeset 762

Timestamp:

Feb 20, 2007, 8:18:48 PM (17 years ago)

Author:

Jari Häkkinen

Message:

Fixes #76. Creating/recreating BLAS support matrix only when needed.

Location:

trunk

Files:

• test/matrix_test.cc (modified) (3 diffs)
• yat/utility/matrix.cc (modified) (8 diffs)
• yat/utility/matrix.h (modified) (3 diffs)

trunk/test/matrix_test.cc

@@ -58 +58 @@
   *error << "Testing matrix class" << std::endl;
   bool ok = true;
+  utility::matrix unit3x3(3,3);
+  for (size_t i=0; i<unit3x3.rows(); ++i)
+    unit3x3(i,i)=1;
 
   *error << "\tcopy constructor and operator!=" << std::endl;
@@ -78 +81 @@
   std::remove("data/tmp_test_matrix.txt");
 
-  *error << "\toperator*(double)" << std::endl;
+  *error << "\toperator*=(double)" << std::endl;
   utility::matrix m4(3,3,1);
   m4 *= 9;
@@ -233 +236 @@
   delete m_nan;
 
+  *error << "\toperator*=(matrix&)" << std::endl;
+  utility::matrix m6(unit3x3);
+  m6 *= m;
+  if (m6!=m) {
+    ok=false;
+    *error << "error operator*=(matrix) 1" << std::endl;
+  }
+  m6 *= unit3x3;
+  if (m6!=m) {
+    ok=false;
+    *error << "error operator*=(matrix) 2" << std::endl;
+  }
+
   if (error!=&std::cerr)
     delete error;

trunk/yat/utility/matrix.cc

@@ -41 +41 @@
 
   matrix::matrix(void)
-    : m_(NULL), view_(NULL)
+    : blas_result_(NULL), m_(NULL), view_(NULL)
   {
   }
@@ -47 +47 @@
 
   matrix::matrix(const size_t& r, const size_t& c, double init_value)
-    : m_(gsl_matrix_alloc(r,c)), view_(NULL)
+    : blas_result_(NULL), m_(gsl_matrix_alloc(r,c)), view_(NULL)
   {
     if (!m_)
@@ -56 +56 @@
 
   matrix::matrix(const matrix& o)
-    : m_(o.create_gsl_matrix_copy()), view_(NULL)
+    : blas_result_(NULL), m_(o.create_gsl_matrix_copy()), view_(NULL)
   {
   }
@@ -63 +63 @@
   matrix::matrix(matrix& m, size_t offset_row, size_t offset_column,
                  size_t n_row, size_t n_column)
+    : blas_result_(NULL)
   {
     view_ = new gsl_matrix_view(gsl_matrix_submatrix(m.m_,
@@ -76 +77 @@
   matrix::matrix(std::istream& is, char sep) 
     throw (utility::IO_error,std::exception)
-    : view_(NULL)
+    : blas_result_(NULL), view_(NULL)
   {
     // read the data file and store in stl vectors (dynamically
@@ -148 +149 @@
     if (view_)
       delete view_;
-    else if (m_)
-      gsl_matrix_free(m_);
+    else {
+      if (blas_result_)
+        gsl_matrix_free(blas_result_);
+      if (m_)
+        gsl_matrix_free(m_);
+    }
+    blas_result_=NULL;
     m_=NULL;
   }
@@ -371 +377 @@
       gsl_matrix_free(m_);
       m_=transposed;
+      if (blas_result_) {
+        gsl_matrix_free(blas_result_);
+        blas_result_=NULL;
+      }
     }
   }
@@ -441 +451 @@
   const matrix& matrix::operator*=(const matrix& other)
   {
-    gsl_matrix* result = gsl_matrix_alloc(rows(),other.columns());
-    if (!result)
-      throw utility::GSL_error("matrix::operator*= failed to allocate memory");
-    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, result);
-    gsl_matrix_free(m_);
-    m_=result;
+    if (!blas_result_) {
+      blas_result_ = gsl_matrix_alloc(rows(),other.columns());
+      if (!blas_result_)
+        throw utility::GSL_error("matrix::operator*= failed to allocate memory");
+    }
+    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0,
+                   blas_result_);
+    gsl_matrix* tmp=m_;
+    m_=blas_result_;
+    blas_result_=tmp;
     return *this;
   }

trunk/yat/utility/matrix.h

@@ -67 +67 @@
   /// the omitted functionality could be included.
   ///
-  /// @todo Maybe it would be smart to create temporary objects need
-  /// for BLAS in constructors. As it is now temporary objects are
-  /// called before BLAS functionality iss used, cf. const matrix&
-  /// matrix::operator*=(const matrix& other)
-  ///
   class matrix
   {
@@ -323 +318 @@
        \brief Swap row \a i and column \a j.
 
+       The matrix must be square.
+
        \throw GSL_error if either index is out of bounds, or if matrix
        is not square.
@@ -457 +454 @@
     gsl_matrix* create_gsl_matrix_copy(void) const;
 
+    // blas_result_ is used to temporarily store result in BLAS calls
+    // and when not NULL it should always have the same size as
+    // m_. Memory is not allocated for blas_result_ until it is
+    // needed.
+    gsl_matrix* blas_result_;
     gsl_matrix* m_;
     gsl_matrix_view* view_;