Nemo  2.3.46
tmatrix.h
Go to the documentation of this file.
1 
30 #ifndef _T_MATRIX_H
31 #define _T_MATRIX_H
32 
33 #ifdef HAS_GSL
34 #include <gsl/gsl_vector.h>
35 #include <gsl/gsl_matrix.h>
36 #include <gsl/gsl_permutation.h>
37 #include <gsl/gsl_linalg.h>
38 #endif
39 
40 #include <string.h>
41 #include <sstream>
42 #include <string>
43 #include <assert.h>
44 #include "output.h"
45 
46 using namespace std;
47 
49 class TMatrix {
50 private:
51 
52  unsigned int _rows, _cols, _length;
53 
54  double* _val;
55 
56 public:
57 
58  TMatrix () : _rows(0), _cols(0), _length(0), _val(0) { }
60  TMatrix (const TMatrix& mat) : _rows(0), _cols(0), _length(0), _val(0)
61  {
62  copy(mat);
63  }
65  TMatrix ( unsigned int rows, unsigned int cols ) : _rows(0), _cols(0), _length(0), _val(0)
66  {
67  _length = rows*cols;
68  _val = new double [_length];
69  _rows = rows; _cols = cols;
70  }
71 
72  ~TMatrix () {if(_val != NULL) delete [] _val;}
73 
75  void copy (const TMatrix& mat)
76  {
77  _rows = mat._rows;
78  _cols = mat._cols;
79  _length = _rows*_cols;
80  if(_val) delete [] _val;
81  _val = new double [_length];
82  memcpy(_val,mat._val,_length*sizeof(double));
83  }
85  void set (unsigned int i, unsigned int j, double val) {
86  if( i*j < _length)
87  _val[i*_cols + j] = val;
88  else
89  error("TMatrix::set overflow!\n");
90  }
92  void assign (double val)
93  {
94  for(unsigned int i = 0; i < _length; ++i) _val[i] = val;
95  }
97  void reset (unsigned int rows, unsigned int cols) {
98  if(_val != NULL) delete [] _val;
99  _length = rows * cols;
100  _val = new double [_length];
101  memset(_val, 0, _length*sizeof(double));
102  _rows = rows; _cols = cols;
103  }
105  void reset (unsigned int rows, unsigned int cols, double* array) {
106  reset(rows, cols);
107  memcpy(_val, array, _length * sizeof(double));
108  }
110  void reset ( )
111  {
112  _rows = 0;
113  _cols = 0;
114  _length = 0;
115  if(_val != NULL) delete [] _val;
116  _val = NULL;
117  }
118 
120  double get (unsigned int i, unsigned int j) {
121  if( !((i+1)*(j+1) > _length) )
122  return _val[i*_cols + j];
123  else
124  fatal("TMatrix::get overflow!\n");
125  return 0;
126  }
128  double* get () const {return _val;}
129  double* getValArray () const {return _val;}
134  unsigned int get_dims (unsigned int* dims) {
135  if(dims != NULL) { dims[0] = _rows; dims[1] = _cols; }
136  return _length;
137  }
139  unsigned int getNbRows ( ) {return _rows;}
140  unsigned int nrows () {return _rows;}
142  unsigned int getNbCols ( ) {return _cols;}
143  unsigned int ncols () {return _cols;}
145  unsigned int length ( ) {return _length;}
151  void getColumnView (unsigned int col, unsigned int n, double* array)
152  {
153  if(col > _cols-1) {
154  error("TMatrix::getColumnView: not that many columns in matrix\n");
155  return;
156  }
157  if(n != _rows) {
158  error("TMatrix::getColumnView: array size not equal to number of rows in matrix\n");
159  return;
160  }
161  for(unsigned int i = 0; i < _rows; ++i)
162  array[i] = _val[i*_cols + col];
163  }
169  void getRowView (unsigned int row, unsigned int n, double* array)
170  {
171  if(row > _rows-1) {
172  error("TMatrix::getRowView: not that many rows in matrix\n");
173  return;
174  }
175  if(n != _cols) {
176  error("TMatrix::getRowView: array size not equal to number of columns in matrix\n");
177  return;
178  }
179  for(unsigned int i = 0, stride = row*_cols; i < _cols; ++i)
180  array[i] = _val[stride + i];
181  }
183  void plus (unsigned int i, unsigned int j, double value)
184  {
185  if( i*j < _length)
186  _val[i*_cols + j] += value;
187  else
188  error("TMatrix::set overflow!\n");
189  }
191  void minus (unsigned int i, unsigned int j, double value)
192  {
193  if( i*j < _length)
194  _val[i*_cols + j] -= value;
195  else
196  error("TMatrix::set overflow!\n");
197  }
199  void multi (unsigned int i, unsigned int j, double value)
200  {
201  if( i*j < _length)
202  _val[i*_cols + j] *= value;
203  else
204  error("TMatrix::set overflow!\n");
205  }
207  void divide (unsigned int i, unsigned int j, double value)
208  {
209  if( i*j < _length)
210  _val[i*_cols + j] /= value;
211  else
212  error("TMatrix::set overflow!\n");
213  }
215  void transpose ()
216  {
217  TMatrix tmp(_cols, _rows);
218 
219  for(unsigned int i = 0; i < _rows; i++)
220  for(unsigned int j = 0; j < _cols; j++)
221  tmp.set(j, i, get(i, j));
222 
223  reset(_cols, _rows, tmp.get());
224  }
226  double colSum (unsigned int col)
227  {
228  assert(col < _cols); //has to be [0, _cols-1]
229  double sum = 0;
230  for (unsigned int i = 0; i < _rows; ++i) {
231  sum += _val[i*_cols + col];
232  }
233  return sum;
234  }
236  double rowSum (unsigned int row)
237  {
238  assert(row < _rows); // has to be [0, _rows-1]
239  double sum = 0;
240  for (unsigned int i = 0; i < _cols; ++i) {
241  sum += _val[row*_cols + i];
242  }
243  return sum;
244  }
245 
246  void show_up()
247  {
248  message("TMatrix dimensions: rows = %i, columns = %i\n",_rows,_cols);
249  for(unsigned int i = 0; i < _rows; i++) {
250  for(unsigned int j = 0; j < _cols; j++)
251  message("%.3f ",_val[i*_cols + j]);
252  message("\n");
253  }
254  }
256  string to_string ()
257  {
258  ostringstream OUT;
259 
260  OUT << "{";
261  for(unsigned int i = 0; i < _rows; i++) {
262  OUT<<"{";
263  for(unsigned int j = 0; j < _cols-1; j++) {
264  OUT<<_val[i*_cols + j]<<",";
265  }
266  OUT<<_val[i*_cols + _cols-1]<<"}";
267  }
268  OUT << "}";
269 
270  return OUT.str();
271  }
272 
273 #ifdef HAS_GSL
274 
277  void get_gsl_matrix(gsl_matrix* M)
278  {
279  if(M->size1 != _rows)
280  fatal("TMatrix::get_gsl_matrix row size of input matrix doesn't match!\n");
281  if(M->size2 != _cols)
282  fatal("TMatrix::get_gsl_matrix col size of input matrix doesn't match!\n");
283  for(unsigned int i = 0; i < _rows; ++i)
284  for(unsigned int j = 0; j < _cols; ++j)
285  gsl_matrix_set(M,i,j,_val[i*_cols + j]);
286  }
290  void set_from_gsl_matrix (gsl_matrix* M)
291  {
292  if(_val != 0) delete [] _val;
293  _rows = M->size1;
294  _cols = M->size2;
295  _length = _rows*_cols;
296  _val = new double [_length];
297  memcpy(_val, M->data, _length*sizeof(double));
298  }
300  void inverse ()
301  {
302  if(_rows != _cols) {
303  error("TMatrix::inverse matrix must be square to invert!\n");
304  return;
305  }
306 
307  gsl_matrix *S = gsl_matrix_alloc(_rows,_cols);
308  get_gsl_matrix(S);
309  gsl_matrix *LU = gsl_matrix_alloc(_rows,_cols);
310  gsl_permutation *P = gsl_permutation_alloc(_rows);
311  int snum;
312 
313  gsl_matrix_memcpy(LU,S);
314  gsl_linalg_LU_decomp(LU,P,&snum);
315  gsl_linalg_LU_invert(LU,P,S);
316 
317  set_from_gsl_matrix(S);
318 
319  gsl_matrix_free(S);
320  gsl_matrix_free(LU);
321  gsl_permutation_free(P);
322  }
323 #endif
324 };
325 
326 #endif
327 
void message(const char *message,...)
Definition: output.cc:40
unsigned int _rows
Definition: tmatrix.h:52
void plus(unsigned int i, unsigned int j, double value)
Adds a value to an element of the matrix.
Definition: tmatrix.h:183
unsigned int ncols()
Definition: tmatrix.h:143
unsigned int _cols
Definition: tmatrix.h:52
void assign(double val)
Assigns a value to all element of the matrix.
Definition: tmatrix.h:92
void minus(unsigned int i, unsigned int j, double value)
Substracts a value from an element of the matrix.
Definition: tmatrix.h:191
unsigned int getNbCols()
Gives the number of columns.
Definition: tmatrix.h:142
unsigned int getNbRows()
Gives the number of rows.
Definition: tmatrix.h:139
unsigned int get_dims(unsigned int *dims)
Accessor to the matrix dimensions.
Definition: tmatrix.h:134
unsigned int length()
Returns the number of elements in the matrix.
Definition: tmatrix.h:145
TMatrix(const TMatrix &mat)
copy constructor.
Definition: tmatrix.h:60
double get(unsigned int i, unsigned int j)
Accessor to element at row i and column j.
Definition: tmatrix.h:120
void fatal(const char *str,...)
Definition: output.cc:90
double * _val
Definition: tmatrix.h:54
int error(const char *str,...)
Definition: output.cc:73
~TMatrix()
Definition: tmatrix.h:72
void multi(unsigned int i, unsigned int j, double value)
Multiply an element of the matrix by a value.
Definition: tmatrix.h:199
double rowSum(unsigned int row)
Sum all elements in a row.
Definition: tmatrix.h:236
void set(unsigned int i, unsigned int j, double val)
Sets element at row i and column j to value val.
Definition: tmatrix.h:85
void transpose()
Transpose the matrix, swaps columns for rows.
Definition: tmatrix.h:215
TMatrix(unsigned int rows, unsigned int cols)
Creates an array of doubles of size = rows*cols.
Definition: tmatrix.h:65
void getRowView(unsigned int row, unsigned int n, double *array)
Gives access to a row of the matrix.
Definition: tmatrix.h:169
A class to handle matrix in params, coerces matrix into a vector of same total size.
Definition: tmatrix.h:49
void show_up()
Definition: tmatrix.h:246
void divide(unsigned int i, unsigned int j, double value)
Divide an element of the matrix by a value.
Definition: tmatrix.h:207
double colSum(unsigned int col)
Sum all elements in a column.
Definition: tmatrix.h:226
void reset(unsigned int rows, unsigned int cols)
Re-allocate the existing matrix with assigned rows and cols dimensions.
Definition: tmatrix.h:97
string to_string()
Writes the matrix into a string in Nemo's matrix input format.
Definition: tmatrix.h:256
void reset()
Reset members to zero state.
Definition: tmatrix.h:110
void copy(const TMatrix &mat)
Copy a matrix.
Definition: tmatrix.h:75
void reset(unsigned int rows, unsigned int cols, double *array)
Reset the existing matrix to the new dimensions and copies the array.
Definition: tmatrix.h:105
void getColumnView(unsigned int col, unsigned int n, double *array)
Gives access to a column of the matrix.
Definition: tmatrix.h:151
TMatrix()
Definition: tmatrix.h:58
unsigned int nrows()
Definition: tmatrix.h:140
double * getValArray() const
Definition: tmatrix.h:129

Generated for Nemo v2.3.0 by  doxygen 1.8.8 --
Catalogued on GSR