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

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

Fixed error message when ,atrix contains elements not covertible to double or nan

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.1 KB
Line 
1// $Id: matrix.cc 567 2006-03-21 20:28:51Z 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     
59      bool ok=true;
60      while(ok) {
61        if(sep=='\0')
62          ok=(ss>>element);
63        else
64          ok=getline(ss, element, sep);
65        if(!ok)
66          break;
67       
68        if(utility::is_double(element)) {
69          v.push_back(atof(element.c_str()));
70        }
71        else if (!element.size() || utility::is_nan(element)) {
72          v.push_back(std::numeric_limits<double>::quiet_NaN());
73        }
74        else {
75          // Jari, this should be communicated with as an exception.
76          std::cerr << "Warning: '" << element
77                    << "' is not accepted as a matrix element." << std::endl;
78        }
79      }           
80      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
81          v.push_back(std::numeric_limits<double>::quiet_NaN());
82      if (!nof_columns)
83        nof_columns=v.size();
84      else if (v.size()!=nof_columns) {
85        std::ostringstream s;
86        s << "matrix::matrix(std::istream&, char) data file error: "
87          << "line " << nof_rows << " has " << v.size()
88          << " columns; expected " << nof_columns << " columns.";
89        throw utility::IO_error(s.str());
90      }
91      data_matrix.push_back(v);
92    }
93
94    // manipulate the state of the stream to be good
95    is.clear(std::ios::goodbit);
96    // convert the data to a gsl matrix
97    m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
98    for(u_int i=0;i<nof_rows;i++)
99      for(u_int j=0;j<nof_columns;j++)
100        gsl_matrix_set( m_, i, j, data_matrix[i][j] );
101  }
102
103
104  matrix::~matrix(void)
105  {
106    if (view_)
107      delete view_;
108    else if (m_)
109      gsl_matrix_free(m_);
110    m_=NULL;
111  }
112
113
114
115  bool matrix::equal(const matrix& other, const double d) const
116  {
117    if (columns()!=other.columns() || rows()!=other.rows())
118      return false;
119    for (size_t i=0; i<rows(); i++)
120      for (size_t j=0; j<columns(); j++)
121        // The two last condition checks are needed for NaN detection
122        if (fabs( (*this)(i,j)-other(i,j) ) > d ||
123            (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j))
124          return false;
125    return true;
126  }
127
128
129
130  gsl_matrix* matrix::create_gsl_matrix_copy(void) const
131  {
132    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
133    gsl_matrix_memcpy(m,m_);  // Jari, a GSL return value is ignored here
134    return m;
135  }
136
137
138
139  matrix matrix::nan(void) const
140  {
141    matrix m(rows(),columns(),1.0);
142    for (size_t i=0; i<rows(); i++)
143      for (size_t j=0; j<columns(); j++) 
144        if (std::isnan(operator()(i,j)))
145          m(i,j)=0;
146    return m;
147  }
148
149
150
151  // Jari, checkout GSL transpose support in GSL manual 8.4.9
152  void matrix::transpose(void)
153  {
154    if (columns()==rows())
155      gsl_matrix_transpose(m_);
156    else {
157      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
158      gsl_matrix_transpose_memcpy(transposed,m_);
159      gsl_matrix_free(m_);
160      m_=transposed;
161    }
162  }
163
164
165
166  const vector matrix::operator*(const vector&) const
167  {
168    std::cerr << "Not implemented:" << std::endl
169              << "   const vector matrix::operator*(const vector&) const"
170              << std::endl;
171    return vector(0);
172  }
173
174
175
176  const matrix& matrix::operator=( const matrix& other )
177  {
178    if ( this != &other ) {
179      if (view_)
180        delete view_;
181      else if (m_)
182        gsl_matrix_free(m_);
183      m_ = other.create_gsl_matrix_copy();
184    }
185    return *this;
186  }
187
188
189
190  const matrix& matrix::operator*=(const matrix& other)
191  {
192    gsl_matrix* result = gsl_matrix_alloc(rows(),other.columns());
193    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, result);
194    gsl_matrix_free(m_);
195    m_=result;
196    return *this;
197  }
198
199
200
201  std::ostream& operator<<(std::ostream& s, const matrix& m)
202  {
203    s.setf(std::ios::dec);
204    s.precision(12);
205    for(size_t i=0, j=0; i<m.rows(); i++)
206      for (j=0; j<m.columns(); j++) {
207        s << m(i,j);
208        if (j<m.columns()-1)
209          s << "\t";
210        else if (i<m.rows()-1)
211          s << "\n";
212      }
213    return s;
214  }
215
216
217}} // of namespace gslapi and namespace thep
Note: See TracBrowser for help on using the repository browser.