source: trunk/src/matrix.cc @ 192

Last change on this file since 192 was 192, checked in by Peter, 17 years ago

modified outstream operator to be compatible with istream constructor

  • 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 192 2004-10-26 18:04:52Z peter $
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  matrix matrix::mul_elements(const matrix& other)
117  {
118    matrix res(other);
119    gsl_matrix_mul_elements( res.m_, other.m_ );
120    return res;
121  }
122
123
124
125  vector matrix::TEMP_col_return(size_t column) const
126  {
127    vector vec(rows());
128    for (size_t i=0; i<rows(); ++i)
129      vec(i)=operator()(i,column);
130    return vec;
131  }
132
133
134
135  vector matrix::operator[](size_t row) const
136  {
137    vector vec(columns());
138    for (size_t i=0; i<columns(); ++i)
139      vec(i)=operator()(row,i);
140    return vec;
141  }
142
143
144
145  matrix matrix::operator*( const matrix &other ) const
146  {
147    matrix res(rows(), other.columns());
148    gsl_linalg_matmult(m_,other.m_,res.m_);
149    return res;
150  }
151
152
153
154  vector matrix::operator*( const vector &other ) const
155  {
156    gsl_vector* gslvec=gsl_vector_alloc(rows());
157    gsl_blas_dgemv(CblasNoTrans, 1.0, m_, other.gsl_vector_pointer(), 0.0,
158                   gslvec );
159    vector res(gslvec);
160    return res;
161  }
162
163
164
165  matrix matrix::operator+( const matrix &other ) const
166  {
167    matrix res( *this );
168    gsl_matrix_add( res.m_, other.m_ );
169    return res;
170  }
171
172
173
174  matrix matrix::operator-( const matrix &other ) const
175  {
176    matrix res( *this );
177    gsl_matrix_sub( res.m_, other.m_ );
178    return res;
179  }
180
181
182 
183  matrix& matrix::operator=( const matrix& other )
184  {
185    if ( this != &other ) {
186      if ( m_ )
187        gsl_matrix_free( m_ );
188      m_ = other.gsl_matrix_copy();
189    }
190    return *this;
191  }
192
193
194
195  matrix matrix::transpose(void) const
196  {
197    matrix res(columns(),rows());
198    gsl_matrix_transpose_memcpy(res.m_,m_);
199    return res;
200  }
201
202
203
204  std::ostream& operator<<(std::ostream& s, const matrix& a)
205  {
206    using namespace std;
207    s.setf( ios::dec );
208    s.precision(12);
209    // Peter, end with an "\n" to make it compatible with stream constructor
210    for(size_t i=0, j=0; i<a.rows(); i++) {
211      for (j=0; j<a.columns()-1; j++) 
212        s << a(i,j) << "\t";
213      s << a(i,j) << "\n";
214    }
215    return s;
216  }
217
218
219}} // of namespace gslapi and namespace thep
Note: See TracBrowser for help on using the repository browser.