source: branches/better_matrix_class/lib/gslapi/matrix.cc @ 415

Last change on this file since 415 was 415, checked in by Jari Häkkinen, 17 years ago

Removed some gslapi copy intensive operators.
Made gslapi assignment operators ignore views, and change them
into normal vectors if assigned.
Cleaning up of code aswell.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.7 KB
Line 
1// $Id: matrix.cc 415 2005-12-01 15:52:36Z jari $
2
3#include <c++_tools/gslapi/matrix.h>
4
5#include <c++_tools/gslapi/vector.h>
6#include <c++_tools/utility/stl_utility.h>
7
8#include <cmath>
9#include <sstream>
10#include <vector>
11
12#include <gsl/gsl_blas.h>
13
14
15namespace theplu {
16namespace gslapi {
17
18
19
20  matrix::matrix(matrix& m, size_t offset_row, size_t offset_column,
21                 size_t n_row, size_t n_column)
22  {
23    // Jari, exception handling needed here. Failure in setting up a
24    // proper gsl_matrix_view is communicated by NULL pointer in the
25    // view structure (cf. GSL manual). How about GSL error state?
26    view_ = new gsl_matrix_view(gsl_matrix_submatrix(m.m_,
27                                                     offset_row,offset_column,
28                                                     n_row,n_column));
29    m_ = &(view_->matrix);
30  }
31
32
33
34  // Constructor that gets data from istream
35  matrix::matrix(std::istream& is) throw (utility::IO_error,std::exception)
36    : view_(NULL)
37  {
38    // read the data file and store in stl vectors (dynamically
39    // expandable)
40    std::vector<std::vector<double> > data_matrix;
41    u_int nof_columns=0;
42    u_int nof_rows = 0;
43    std::vector<double> v;
44    for (nof_rows = 0; utility::read_to_double(is, v); nof_rows++) {
45      // Ignoring empty lines
46      if (!v.size()) {
47        nof_rows--;
48        continue;
49      }
50      if (!nof_columns)
51        nof_columns=v.size();
52      else if (v.size()!=nof_columns) {
53        std::ostringstream s;
54        s << "matrix::matrix(std::istream&) data file error: "
55          << "line" << nof_rows+1 << " has " << v.size()
56          << " columns; expected " << nof_columns << " columns.";
57        throw utility::IO_error(s.str());
58      }
59      data_matrix.push_back(v);
60    }
61
62    // manipulate the state of the stream to be good
63    is.clear(std::ios::goodbit);
64    // convert the data to a gsl matrix
65    m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
66    for(u_int i=0;i<nof_rows;i++)
67      for(u_int j=0;j<nof_columns;j++)
68        gsl_matrix_set( m_, i, j, data_matrix[i][j] );
69  }
70
71
72
73  matrix::~matrix(void)
74  {
75    if (view_)
76      delete view_;
77    if (m_) {
78      gsl_matrix_free(m_);
79      m_=NULL;
80    }
81  }
82
83
84
85  bool matrix::equal(const matrix& other, const double d) const
86  {
87    if (columns()!=other.columns() || rows()!=other.rows())
88      return false;
89    for (size_t i=0; i<rows(); i++)
90      for (size_t j=0; j<columns(); j++)
91        // The two last condition checks are needed for NaN detection
92        if (fabs( (*this)(i,j)-other(i,j) ) > d ||
93            (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j))
94          return false;
95    return true;
96  }
97
98
99
100  gsl_matrix* matrix::create_gsl_matrix_copy(void) const
101  {
102    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
103    gsl_matrix_memcpy(m,m_);  // Jari, a GSL return value is ignored here
104    return m;
105  }
106
107
108
109  // Jari, checkout GSL transpose support in GSL manual 8.4.9
110  void matrix::transpose(void)
111  {
112    if (columns()==rows())
113      gsl_matrix_transpose(m_);
114    else {
115      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
116      gsl_matrix_transpose_memcpy(transposed,m_);
117      gsl_matrix_free(m_);
118      m_=transposed;
119    }
120  }
121
122
123
124  const matrix& matrix::operator=( const matrix& other )
125  {
126    if ( this != &other ) {
127      if (view_)
128        delete view_;
129      else if (m_)
130        gsl_matrix_free(m_);
131      m_ = other.create_gsl_matrix_copy();
132    }
133    return *this;
134  }
135
136
137
138  const matrix& matrix::operator*=(const matrix& other)
139  {
140    gsl_matrix* result = gsl_matrix_alloc(rows(),other.columns());
141    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, result);
142    gsl_matrix_free(m_);
143    m_=result;
144    return *this;
145  }
146
147
148
149  std::ostream& operator<<(std::ostream& s, const matrix& m)
150  {
151    using namespace std;
152    s.setf( ios::dec );
153    s.precision(12);
154    for(size_t i=0, j=0; i<m.rows(); i++)
155      for (j=0; j<m.columns(); j++) {
156        s << m(i,j);
157        if (j<m.columns()-1)
158          s << "\t";
159        else if (i<m.rows()-1)
160          s << "\n";
161      }
162    return s;
163  }
164
165
166}} // of namespace gslapi and namespace thep
Note: See TracBrowser for help on using the repository browser.