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

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

cosmetic change of printout

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 4.8 KB
Line 
1// $Id: matrix.cc 444 2005-12-15 16:20:40Z 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: '" << vec_str[i]
77          //                    << "' is not an integer." << 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  // Jari, checkout GSL transpose support in GSL manual 8.4.9
140  void matrix::transpose(void)
141  {
142    if (columns()==rows())
143      gsl_matrix_transpose(m_);
144    else {
145      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
146      gsl_matrix_transpose_memcpy(transposed,m_);
147      gsl_matrix_free(m_);
148      m_=transposed;
149    }
150  }
151
152
153
154  const vector matrix::operator*(const vector&) const
155  {
156    std::cerr << "Not implemented:" << std::endl
157              << "   const vector matrix::operator*(const vector&) const"
158              << std::endl;
159    return vector(0);
160  }
161
162
163
164  const matrix& matrix::operator=( const matrix& other )
165  {
166    if ( this != &other ) {
167      if (view_)
168        delete view_;
169      else if (m_)
170        gsl_matrix_free(m_);
171      m_ = other.create_gsl_matrix_copy();
172    }
173    return *this;
174  }
175
176
177
178  const matrix& matrix::operator*=(const matrix& other)
179  {
180    gsl_matrix* result = gsl_matrix_alloc(rows(),other.columns());
181    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, result);
182    gsl_matrix_free(m_);
183    m_=result;
184    return *this;
185  }
186
187
188
189  std::ostream& operator<<(std::ostream& s, const matrix& m)
190  {
191    s.setf(std::ios::dec);
192    s.precision(12);
193    for(size_t i=0, j=0; i<m.rows(); i++)
194      for (j=0; j<m.columns(); j++) {
195        s << m(i,j);
196        if (j<m.columns()-1)
197          s << "\t";
198        else if (i<m.rows()-1)
199          s << "\n";
200      }
201    return s;
202  }
203
204
205}} // of namespace gslapi and namespace thep
Note: See TracBrowser for help on using the repository browser.