source: trunk/src/matrix.cc @ 240

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

Changed mul_elements interface, and fixed serious bug in it.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 3.9 KB
Line 
1// $Id: matrix.cc 240 2005-02-22 11:13:00Z jari $
2
3#include <iostream>
4#include <sstream>
5#include <string>
6#include <vector>
7
8#include <gsl/gsl_blas.h>
9#include <gsl/gsl_linalg.h>
10
11#include "matrix.h"
12#include "vector.h"
13
14
15namespace theplu {
16namespace gslapi {
17
18
19
20  matrix::matrix(void)
21    : m_(NULL)
22  {
23  }
24
25
26
27  matrix::matrix(const size_t& r, const size_t& c, double init_value)
28  {
29    m_ = gsl_matrix_alloc(r,c);
30    set_all(init_value);
31  }
32
33
34
35  matrix::matrix(const matrix& other)
36  {
37    m_ = other.gsl_matrix_copy();
38  }
39
40
41
42  // Constructor that gets data from istream
43  matrix::matrix(std::istream& is)
44  {
45    using namespace std;
46
47    // read the data file and store in stl vectors (dynamically
48    // expandable)
49    std::vector<std::vector<double> > data_matrix;
50    u_int nof_columns=0;
51    u_int nof_rows = 0;
52    string s;
53    for (nof_rows = 0; getline(is, s, '\n'); nof_rows++) {
54      istringstream line(s);
55      std::vector<double> v;
56      string tmp_string;
57      while (line >> tmp_string) {
58        if(!is.good()) {
59          cerr << "matrix::matrix(std::istream&): "
60               << "error reading data file!" << endl;
61          exit(1);
62        }
63        double t=atof(tmp_string.c_str());
64        // Peter, should check that it actually was a double!!!!!
65        v.push_back(t);
66      }
67   
68      // Ignoring empty lines
69      if(!v.size()){
70        nof_rows--;
71        continue;
72      }
73     
74      if(nof_columns==0)
75        nof_columns=v.size();
76      else if(v.size()!=nof_columns) {
77        cerr << "matrix::matrix(std::istream&) data file error: "
78             << "line" << nof_rows+1 << " has " << v.size() 
79             << " columns; expected " << nof_columns
80             << " columns"
81             << endl; 
82        exit(1);
83      }
84      data_matrix.push_back(v);
85    }
86    // manipulate the state of the stream to be good
87    is.clear(std::ios::goodbit);
88    // convert the data to a gsl matrix
89    m_ = gsl_matrix_alloc ( nof_rows, nof_columns );
90    for(u_int i=0;i<nof_rows;i++) 
91      for(u_int j=0;j<nof_columns;j++) 
92        gsl_matrix_set( m_, i, j, data_matrix[i][j] );
93  }
94
95
96
97  matrix::~matrix(void)
98  {
99    if (m_) {
100      gsl_matrix_free(m_);
101      m_=NULL;
102    }
103  }
104
105
106
107  gsl_matrix* matrix::gsl_matrix_copy(void) const
108  {
109    gsl_matrix* m = gsl_matrix_alloc(rows(),columns());
110    gsl_matrix_memcpy(m,m_);
111    return m;
112  }
113
114
115
116  void matrix::mul_elements(const matrix& other)
117  {
118    gsl_matrix_mul_elements( m_, other.m_ );
119  }
120
121
122
123  vector matrix::TEMP_col_return(size_t column) const
124  {
125    vector vec(rows());
126    for (size_t i=0; i<rows(); ++i)
127      vec(i)=operator()(i,column);
128    return vec;
129  }
130
131
132
133  vector matrix::operator[](size_t row) const
134  {
135    vector vec(columns());
136    for (size_t i=0; i<columns(); ++i)
137      vec(i)=operator()(row,i);
138    return vec;
139  }
140
141
142
143  matrix matrix::operator*( const matrix &other ) const
144  {
145    matrix res(rows(), other.columns());
146    gsl_linalg_matmult(m_,other.m_,res.m_);
147    return res;
148  }
149
150
151
152  vector matrix::operator*( const vector &other ) const
153  {
154    gsl_vector* gslvec=gsl_vector_alloc(rows());
155    gsl_blas_dgemv(CblasNoTrans, 1.0, m_, other.gsl_vector_pointer(), 0.0,
156                   gslvec );
157    vector res(gslvec);
158    return res;
159  }
160
161
162
163  matrix matrix::operator+( const matrix &other ) const
164  {
165    matrix res( *this );
166    gsl_matrix_add( res.m_, other.m_ );
167    return res;
168  }
169
170
171
172  matrix matrix::operator-( const matrix &other ) const
173  {
174    matrix res( *this );
175    gsl_matrix_sub( res.m_, other.m_ );
176    return res;
177  }
178
179
180 
181  matrix& matrix::operator=( const matrix& other )
182  {
183    if ( this != &other ) {
184      if ( m_ )
185        gsl_matrix_free( m_ );
186      m_ = other.gsl_matrix_copy();
187    }
188    return *this;
189  }
190
191
192
193  matrix matrix::transpose(void) const
194  {
195    matrix res(columns(),rows());
196    gsl_matrix_transpose_memcpy(res.m_,m_);
197    return res;
198  }
199
200
201
202  std::ostream& operator<<(std::ostream& s, const matrix& a)
203  {
204    using namespace std;
205    s.setf( ios::dec );
206    s.precision(12);
207    // Peter, end with an "\n" to make it compatible with stream constructor
208    for(size_t i=0, j=0; i<a.rows(); i++) {
209      for (j=0; j<a.columns()-1; j++) 
210        s << a(i,j) << "\t";
211      s << a(i,j) << "\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.