Changeset 762


Ignore:
Timestamp:
Feb 20, 2007, 8:18:48 PM (15 years ago)
Author:
Jari Häkkinen
Message:

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

Location:
trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/test/matrix_test.cc

    r680 r762  
    5858  *error << "Testing matrix class" << std::endl;
    5959  bool ok = true;
     60  utility::matrix unit3x3(3,3);
     61  for (size_t i=0; i<unit3x3.rows(); ++i)
     62    unit3x3(i,i)=1;
    6063
    6164  *error << "\tcopy constructor and operator!=" << std::endl;
     
    7881  std::remove("data/tmp_test_matrix.txt");
    7982
    80   *error << "\toperator*(double)" << std::endl;
     83  *error << "\toperator*=(double)" << std::endl;
    8184  utility::matrix m4(3,3,1);
    8285  m4 *= 9;
     
    233236  delete m_nan;
    234237
     238  *error << "\toperator*=(matrix&)" << std::endl;
     239  utility::matrix m6(unit3x3);
     240  m6 *= m;
     241  if (m6!=m) {
     242    ok=false;
     243    *error << "error operator*=(matrix) 1" << std::endl;
     244  }
     245  m6 *= unit3x3;
     246  if (m6!=m) {
     247    ok=false;
     248    *error << "error operator*=(matrix) 2" << std::endl;
     249  }
     250
    235251  if (error!=&std::cerr)
    236252    delete error;
  • trunk/yat/utility/matrix.cc

    r759 r762  
    4141
    4242  matrix::matrix(void)
    43     : m_(NULL), view_(NULL)
     43    : blas_result_(NULL), m_(NULL), view_(NULL)
    4444  {
    4545  }
     
    4747
    4848  matrix::matrix(const size_t& r, const size_t& c, double init_value)
    49     : m_(gsl_matrix_alloc(r,c)), view_(NULL)
     49    : blas_result_(NULL), m_(gsl_matrix_alloc(r,c)), view_(NULL)
    5050  {
    5151    if (!m_)
     
    5656
    5757  matrix::matrix(const matrix& o)
    58     : m_(o.create_gsl_matrix_copy()), view_(NULL)
     58    : blas_result_(NULL), m_(o.create_gsl_matrix_copy()), view_(NULL)
    5959  {
    6060  }
     
    6363  matrix::matrix(matrix& m, size_t offset_row, size_t offset_column,
    6464                 size_t n_row, size_t n_column)
     65    : blas_result_(NULL)
    6566  {
    6667    view_ = new gsl_matrix_view(gsl_matrix_submatrix(m.m_,
     
    7677  matrix::matrix(std::istream& is, char sep)
    7778    throw (utility::IO_error,std::exception)
    78     : view_(NULL)
     79    : blas_result_(NULL), view_(NULL)
    7980  {
    8081    // read the data file and store in stl vectors (dynamically
     
    148149    if (view_)
    149150      delete view_;
    150     else if (m_)
    151       gsl_matrix_free(m_);
     151    else {
     152      if (blas_result_)
     153        gsl_matrix_free(blas_result_);
     154      if (m_)
     155        gsl_matrix_free(m_);
     156    }
     157    blas_result_=NULL;
    152158    m_=NULL;
    153159  }
     
    371377      gsl_matrix_free(m_);
    372378      m_=transposed;
     379      if (blas_result_) {
     380        gsl_matrix_free(blas_result_);
     381        blas_result_=NULL;
     382      }
    373383    }
    374384  }
     
    441451  const matrix& matrix::operator*=(const matrix& other)
    442452  {
    443     gsl_matrix* result = gsl_matrix_alloc(rows(),other.columns());
    444     if (!result)
    445       throw utility::GSL_error("matrix::operator*= failed to allocate memory");
    446     gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, result);
    447     gsl_matrix_free(m_);
    448     m_=result;
     453    if (!blas_result_) {
     454      blas_result_ = gsl_matrix_alloc(rows(),other.columns());
     455      if (!blas_result_)
     456        throw utility::GSL_error("matrix::operator*= failed to allocate memory");
     457    }
     458    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0,
     459                   blas_result_);
     460    gsl_matrix* tmp=m_;
     461    m_=blas_result_;
     462    blas_result_=tmp;
    449463    return *this;
    450464  }
  • trunk/yat/utility/matrix.h

    r759 r762  
    6767  /// the omitted functionality could be included.
    6868  ///
    69   /// @todo Maybe it would be smart to create temporary objects need
    70   /// for BLAS in constructors. As it is now temporary objects are
    71   /// called before BLAS functionality iss used, cf. const matrix&
    72   /// matrix::operator*=(const matrix& other)
    73   ///
    7469  class matrix
    7570  {
     
    323318       \brief Swap row \a i and column \a j.
    324319
     320       The matrix must be square.
     321
    325322       \throw GSL_error if either index is out of bounds, or if matrix
    326323       is not square.
     
    457454    gsl_matrix* create_gsl_matrix_copy(void) const;
    458455
     456    // blas_result_ is used to temporarily store result in BLAS calls
     457    // and when not NULL it should always have the same size as
     458    // m_. Memory is not allocated for blas_result_ until it is
     459    // needed.
     460    gsl_matrix* blas_result_;
    459461    gsl_matrix* m_;
    460462    gsl_matrix_view* view_;
Note: See TracChangeset for help on using the changeset viewer.