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

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

modified includes in tests

  • 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 301 2005-04-30 13:39:27Z 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; utility::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.