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

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

Made output from matrix and vector coherrent.

  • 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 418 2005-12-01 23:52:12Z 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    else if (m_)
78      gsl_matrix_free(m_);
79    m_=NULL;
80  }
81
82
83
84  bool matrix::equal(const matrix& other, const double d) const
85  {
86    if (columns()!=other.columns() || rows()!=other.rows())
87      return false;
88    for (size_t i=0; i<rows(); i++)
89      for (size_t j=0; j<columns(); j++)
90        // The two last condition checks are needed for NaN detection
91        if (fabs( (*this)(i,j)-other(i,j) ) > d ||
92            (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j))
93          return false;
94    return true;
95  }
96
97
98
99  gsl_matrix* matrix::create_gsl_matrix_copy(void) const
100  {
101    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
102    gsl_matrix_memcpy(m,m_);  // Jari, a GSL return value is ignored here
103    return m;
104  }
105
106
107
108  // Jari, checkout GSL transpose support in GSL manual 8.4.9
109  void matrix::transpose(void)
110  {
111    if (columns()==rows())
112      gsl_matrix_transpose(m_);
113    else {
114      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
115      gsl_matrix_transpose_memcpy(transposed,m_);
116      gsl_matrix_free(m_);
117      m_=transposed;
118    }
119  }
120
121
122
123  const matrix& matrix::operator=( const matrix& other )
124  {
125    if ( this != &other ) {
126      if (view_)
127        delete view_;
128      else if (m_)
129        gsl_matrix_free(m_);
130      m_ = other.create_gsl_matrix_copy();
131    }
132    return *this;
133  }
134
135
136
137  const matrix& matrix::operator*=(const matrix& other)
138  {
139    gsl_matrix* result = gsl_matrix_alloc(rows(),other.columns());
140    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, result);
141    gsl_matrix_free(m_);
142    m_=result;
143    return *this;
144  }
145
146
147
148  std::ostream& operator<<(std::ostream& s, const matrix& m)
149  {
150    s.setf(std::ios::dec);
151    s.precision(12);
152    for(size_t i=0, j=0; i<m.rows(); i++)
153      for (j=0; j<m.columns(); j++) {
154        s << m(i,j);
155        if (j<m.columns()-1)
156          s << "\t";
157        else if (i<m.rows()-1)
158          s << "\n";
159      }
160    return s;
161  }
162
163
164}} // of namespace gslapi and namespace thep
Note: See TracBrowser for help on using the repository browser.