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

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

file structure modifications. NOTE, this revision is not working, please wait for the next...

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