source: trunk/c++_tools/gslapi/matrix.cc @ 593

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

Fixed std includes to compile with g++ 4.1.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 6.0 KB
Line 
1// $Id: matrix.cc 593 2006-08-25 12:59:21Z markus $
2
3/*
4  Copyright (C) 2003 Daniel Dalevi, Peter Johansson
5  Copyright (C) 2004 Jari Häkkinen, Peter Johansson
6  Copyright (C) 2005, 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
7
8  This file is part of the thep c++ tools library,
9                                http://lev.thep.lu.se/trac/c++_tools
10
11  The c++ tools library is free software; you can redistribute it
12  and/or modify it under the terms of the GNU General Public License
13  as published by the Free Software Foundation; either version 2 of
14  the License, or (at your option) any later version.
15
16  The c++ tools library is distributed in the hope that it will be
17  useful, but WITHOUT ANY WARRANTY; without even the implied warranty
18  of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
19  General Public License for more details.
20
21  You should have received a copy of the GNU General Public License
22  along with this program; if not, write to the Free Software
23  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
24  02111-1307, USA.
25*/
26
27#include <c++_tools/gslapi/matrix.h>
28
29#include <c++_tools/gslapi/vector.h>
30#include <c++_tools/utility/stl_utility.h>
31#include <c++_tools/utility/utility.h>
32
33#include <cmath>
34#include <sstream>
35#include <vector>
36
37#include <gsl/gsl_blas.h>
38
39
40namespace theplu {
41namespace gslapi {
42
43
44
45  matrix::matrix(matrix& m, size_t offset_row, size_t offset_column,
46                 size_t n_row, size_t n_column)
47  {
48    // Jari, exception handling needed here. Failure in setting up a
49    // proper gsl_matrix_view is communicated by NULL pointer in the
50    // view structure (cf. GSL manual). How about GSL error state?
51    view_ = new gsl_matrix_view(gsl_matrix_submatrix(m.m_,
52                                                     offset_row,offset_column,
53                                                     n_row,n_column));
54    m_ = &(view_->matrix);
55  }
56
57
58  // Constructor that gets data from istream
59  matrix::matrix(std::istream& is, char sep) 
60    throw (utility::IO_error,std::exception)
61    : view_(NULL)
62  {
63    // Markus to Jari, somewhere we should check that quiet_NaNs are supported
64    // std::numeric_limits<double>::has_quiet_NaN has to be true.
65    // Also in vector
66
67    // read the data file and store in stl vectors (dynamically
68    // expandable)
69    std::vector<std::vector<double> > data_matrix;
70    u_int nof_columns=0;
71    u_int nof_rows = 0;
72    std::string line;
73    while(getline(is, line, '\n')){
74      // Ignoring empty lines
75      if (!line.size()) {
76        continue;
77      }
78      nof_rows++;
79      std::vector<double> v;
80      std::string element;
81      std::stringstream ss(line);
82     
83      bool ok=true;
84      while(ok) {
85        if(sep=='\0')
86          ok=(ss>>element);
87        else
88          ok=getline(ss, element, sep);
89        if(!ok)
90          break;
91       
92        if(utility::is_double(element)) {
93          v.push_back(atof(element.c_str()));
94        }
95        else if (!element.size() || utility::is_nan(element)) {
96          v.push_back(std::numeric_limits<double>::quiet_NaN());
97        }
98        else {
99          // Jari, this should be communicated with as an exception.
100          std::cerr << "Warning: '" << element
101                    << "' is not accepted as a matrix element." << std::endl;
102        }
103      }           
104      if(sep!='\0' && line[line.size()-1]==sep) // add NaN for final separator
105          v.push_back(std::numeric_limits<double>::quiet_NaN());
106      if (!nof_columns)
107        nof_columns=v.size();
108      else if (v.size()!=nof_columns) {
109        std::ostringstream s;
110        s << "matrix::matrix(std::istream&, char) data file error: "
111          << "line " << nof_rows << " has " << v.size()
112          << " columns; expected " << nof_columns << " columns.";
113        throw utility::IO_error(s.str());
114      }
115      data_matrix.push_back(v);
116    }
117
118    // manipulate the state of the stream to be good
119    is.clear(std::ios::goodbit);
120    // convert the data to a gsl matrix
121    m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
122    for(u_int i=0;i<nof_rows;i++)
123      for(u_int j=0;j<nof_columns;j++)
124        gsl_matrix_set( m_, i, j, data_matrix[i][j] );
125  }
126
127
128  matrix::~matrix(void)
129  {
130    if (view_)
131      delete view_;
132    else if (m_)
133      gsl_matrix_free(m_);
134    m_=NULL;
135  }
136
137
138
139  bool matrix::equal(const matrix& other, const double d) const
140  {
141    if (columns()!=other.columns() || rows()!=other.rows())
142      return false;
143    for (size_t i=0; i<rows(); i++)
144      for (size_t j=0; j<columns(); j++)
145        // The two last condition checks are needed for NaN detection
146        if (fabs( (*this)(i,j)-other(i,j) ) > d ||
147            (*this)(i,j)!=(*this)(i,j) || other(i,j)!=other(i,j))
148          return false;
149    return true;
150  }
151
152
153
154  gsl_matrix* matrix::create_gsl_matrix_copy(void) const
155  {
156    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
157    gsl_matrix_memcpy(m,m_);  // Jari, a GSL return value is ignored here
158    return m;
159  }
160
161
162
163  matrix matrix::nan(void) const
164  {
165    matrix m(rows(),columns(),1.0);
166    for (size_t i=0; i<rows(); i++)
167      for (size_t j=0; j<columns(); j++) 
168        if (std::isnan(operator()(i,j)))
169          m(i,j)=0;
170    return m;
171  }
172
173
174
175  // Jari, checkout GSL transpose support in GSL manual 8.4.9
176  void matrix::transpose(void)
177  {
178    if (columns()==rows())
179      gsl_matrix_transpose(m_);
180    else {
181      gsl_matrix* transposed = gsl_matrix_alloc(columns(),rows());
182      gsl_matrix_transpose_memcpy(transposed,m_);
183      gsl_matrix_free(m_);
184      m_=transposed;
185    }
186  }
187
188
189
190  const vector matrix::operator*(const vector&) const
191  {
192    std::cerr << "Not implemented:" << std::endl
193              << "   const vector matrix::operator*(const vector&) const"
194              << std::endl;
195    return vector(0);
196  }
197
198
199
200  const matrix& matrix::operator=( const matrix& other )
201  {
202    if ( this != &other ) {
203      if (view_)
204        delete view_;
205      else if (m_)
206        gsl_matrix_free(m_);
207      m_ = other.create_gsl_matrix_copy();
208    }
209    return *this;
210  }
211
212
213
214  const matrix& matrix::operator*=(const matrix& other)
215  {
216    gsl_matrix* result = gsl_matrix_alloc(rows(),other.columns());
217    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m_, other.m_, 0.0, result);
218    gsl_matrix_free(m_);
219    m_=result;
220    return *this;
221  }
222
223
224
225  std::ostream& operator<<(std::ostream& s, const matrix& m)
226  {
227    s.setf(std::ios::dec);
228    s.precision(12);
229    for(size_t i=0, j=0; i<m.rows(); i++)
230      for (j=0; j<m.columns(); j++) {
231        s << m(i,j);
232        if (j<m.columns()-1)
233          s << s.fill();
234        else if (i<m.rows()-1)
235          s << "\n";
236      }
237    return s;
238  }
239
240
241}} // of namespace gslapi and namespace thep
Note: See TracBrowser for help on using the repository browser.