Changeset 762 for trunk/yat/utility/matrix.cc

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.

File:

• trunk/yat/utility/matrix.cc (modified) (8 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;
   }