source: trunk/lib/gslapi/matrix.cc @ 434

Last change on this file since 434 was 434, checked in by Markus Ringnér, 17 years ago

Modified istream constructors to support both the old format and CSV files. Also support has been added for NaN:s

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