Changeset 762 for trunk/yat/utility


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/yat/utility
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • 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.