Changeset 762 for trunk/yat

Timestamp:

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

Author:

Jari Häkkinen

Message:

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

Location:

trunk/yat/utility

Files:

• matrix.cc (modified) (8 diffs)
• matrix.h (modified) (3 diffs)

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_;