Essentia  2.1-beta5-dev
tnt_cmat.h
Go to the documentation of this file.
1 /*
2 *
3 * Template Numerical Toolkit (TNT)
4 *
5 * Mathematical and Computational Sciences Division
6 * National Institute of Technology,
7 * Gaithersburg, MD USA
8 *
9 *
10 * This software was developed at the National Institute of Standards and
11 * Technology (NIST) by employees of the Federal Government in the course
12 * of their official duties. Pursuant to title 17 Section 105 of the
13 * United States Code, this software is not subject to copyright protection
14 * and is in the public domain. NIST assumes no responsibility whatsoever for
15 * its use by other parties, and makes no guarantees, expressed or implied,
16 * about its quality, reliability, or any other characteristic.
17 *
18 */
19 
20 
21 // C compatible matrix: row-oriented, 0-based [i][j] and 1-based (i,j) indexing
22 //
23 
24 #ifndef TNT_CMAT_H
25 #define TNT_CMAT_H
26 
27 #include "tnt_subscript.h"
28 #include "tnt_vec.h"
29 #include <cstdlib>
30 #include <cassert>
31 #include <iostream>
32 #include <sstream>
33 
34 namespace TNT
35 {
36 
37 
38 template <class T>
39 class Matrix
40 {
41 
42 
43  public:
44 
46  typedef T value_type;
47  typedef T element_type;
48  typedef T* pointer;
49  typedef T* iterator;
50  typedef T& reference;
51  typedef const T* const_iterator;
52  typedef const T& const_reference;
53 
54  Subscript lbound() const { return 1;}
55 
56  protected:
59  Subscript mn_; // total size
60  T* v_;
61  T** row_;
62  T* vm1_ ; // these point to the same data, but are 1-based
63  T** rowm1_;
64 
65  // internal helper function to create the array
66  // of row pointers
67 
69  {
70  mn_ = M*N;
71  m_ = M;
72  n_ = N;
73 
74  v_ = new T[mn_];
75  row_ = new T*[M];
76  rowm1_ = new T*[M];
77 
78  assert(v_ != NULL);
79  assert(row_ != NULL);
80  assert(rowm1_ != NULL);
81 
82  T* p = v_;
83  vm1_ = v_ - 1;
84  for (Subscript i=0; i<M; i++)
85  {
86  row_[i] = p;
87  rowm1_[i] = p-1;
88  p += N ;
89 
90  }
91 
92  rowm1_ -- ; // compensate for 1-based offset
93  }
94 
95  void copy(const T* v)
96  {
97  Subscript N = m_ * n_;
98  Subscript i;
99 
100 #ifdef TNT_UNROLL_LOOPS
101  Subscript Nmod4 = N & 3;
102  Subscript N4 = N - Nmod4;
103 
104  for (i=0; i<N4; i+=4)
105  {
106  v_[i] = v[i];
107  v_[i+1] = v[i+1];
108  v_[i+2] = v[i+2];
109  v_[i+3] = v[i+3];
110  }
111 
112  for (i=N4; i< N; i++)
113  v_[i] = v[i];
114 #else
115 
116  for (i=0; i< N; i++)
117  v_[i] = v[i];
118 #endif
119  }
120 
121  void set(const T& val)
122  {
123  Subscript N = m_ * n_;
124  Subscript i;
125 
126 #ifdef TNT_UNROLL_LOOPS
127  Subscript Nmod4 = N & 3;
128  Subscript N4 = N - Nmod4;
129 
130  for (i=0; i<N4; i+=4)
131  {
132  v_[i] = val;
133  v_[i+1] = val;
134  v_[i+2] = val;
135  v_[i+3] = val;
136  }
137 
138  for (i=N4; i< N; i++)
139  v_[i] = val;
140 #else
141 
142  for (i=0; i< N; i++)
143  v_[i] = val;
144 
145 #endif
146  }
147 
148 
149 
150  void destroy()
151  {
152  /* do nothing, if no memory has been previously allocated */
153  if (v_ == NULL) return ;
154 
155  /* if we are here, then matrix was previously allocated */
156  if (v_ != NULL) delete [] (v_);
157  if (row_ != NULL) delete [] (row_);
158 
159  /* return rowm1_ back to original value */
160  rowm1_ ++;
161  if (rowm1_ != NULL ) delete [] (rowm1_);
162  }
163 
164 
165  public:
166 
167  operator T**(){ return row_; }
168  operator T**() const { return row_; }
169 
170 
171  Subscript size() const { return mn_; }
172 
173  // constructors
174 
175  Matrix() : m_(0), n_(0), mn_(0), v_(0), row_(0), vm1_(0), rowm1_(0) {};
176 
177  Matrix(const Matrix<T> &A)
178  {
179  initialize(A.m_, A.n_);
180  copy(A.v_);
181  }
182 
183  Matrix(Subscript M, Subscript N, const T& value = T())
184  {
185  initialize(M,N);
186  set(value);
187  }
188 
189  Matrix(Subscript M, Subscript N, const T* v)
190  {
191  initialize(M,N);
192  copy(v);
193  }
194 
195  Matrix(Subscript M, Subscript N, const char *s)
196  {
197  initialize(M,N);
198  //std::istrstream ins(s);
199  std::istringstream ins(s);
200 
201  Subscript i, j;
202 
203  for (i=0; i<M; i++)
204  for (j=0; j<N; j++)
205  ins >> row_[i][j];
206  }
207 
208  // destructor
209  //
211  {
212  destroy();
213  }
214 
215 
216  // reallocating
217  //
219  {
220  if (num_rows() == M && num_cols() == N)
221  return *this;
222 
223  destroy();
224  initialize(M,N);
225 
226  return *this;
227  }
228 
229 
230 
231 
232  // assignments
233  //
235  {
236  if (v_ == A.v_)
237  return *this;
238 
239  if (m_ == A.m_ && n_ == A.n_) // no need to re-alloc
240  copy(A.v_);
241 
242  else
243  {
244  destroy();
245  initialize(A.m_, A.n_);
246  copy(A.v_);
247  }
248 
249  return *this;
250  }
251 
252  Matrix<T>& operator=(const T& scalar)
253  {
254  set(scalar);
255  return *this;
256  }
257 
258 
260  {
261 #ifdef TNT_BOUNDS_CHECK
262  assert( d >= 1);
263  assert( d <= 2);
264 #endif
265  return (d==1) ? m_ : ((d==2) ? n_ : 0);
266  }
267 
268  Subscript num_rows() const { return m_; }
269  Subscript num_cols() const { return n_; }
270 
271 
272 
273 
274  inline T* operator[](Subscript i)
275  {
276 #ifdef TNT_BOUNDS_CHECK
277  assert(0<=i);
278  assert(i < m_) ;
279 #endif
280  return row_[i];
281  }
282 
283  inline const T* operator[](Subscript i) const
284  {
285 #ifdef TNT_BOUNDS_CHECK
286  assert(0<=i);
287  assert(i < m_) ;
288 #endif
289  return row_[i];
290  }
291 
292  inline reference operator()(Subscript i)
293  {
294 #ifdef TNT_BOUNDS_CHECK
295  assert(1<=i);
296  assert(i <= mn_) ;
297 #endif
298  return vm1_[i];
299  }
300 
301  inline const_reference operator()(Subscript i) const
302  {
303 #ifdef TNT_BOUNDS_CHECK
304  assert(1<=i);
305  assert(i <= mn_) ;
306 #endif
307  return vm1_[i];
308  }
309 
310 
311 
312  inline reference operator()(Subscript i, Subscript j)
313  {
314 #ifdef TNT_BOUNDS_CHECK
315  assert(1<=i);
316  assert(i <= m_) ;
317  assert(1<=j);
318  assert(j <= n_);
319 #endif
320  return rowm1_[i][j];
321  }
322 
323 
324 
325  inline const_reference operator() (Subscript i, Subscript j) const
326  {
327 #ifdef TNT_BOUNDS_CHECK
328  assert(1<=i);
329  assert(i <= m_) ;
330  assert(1<=j);
331  assert(j <= n_);
332 #endif
333  return rowm1_[i][j];
334  }
335 
336 
337 
338 
339 };
340 
341 
342 /* *************************** I/O ********************************/
343 
344 template <class T>
345 std::ostream& operator<<(std::ostream &s, const Matrix<T> &A)
346 {
347  Subscript M=A.num_rows();
348  Subscript N=A.num_cols();
349 
350  s << M << " " << N << "\n";
351 
352  for (Subscript i=0; i<M; i++)
353  {
354  for (Subscript j=0; j<N; j++)
355  {
356  s << A[i][j] << " ";
357  }
358  s << "\n";
359  }
360 
361 
362  return s;
363 }
364 
365 template <class T>
366 std::istream& operator>>(std::istream &s, Matrix<T> &A)
367 {
368 
369  Subscript M, N;
370 
371  s >> M >> N;
372 
373  if ( !(M == A.num_rows() && N == A.num_cols() ))
374  {
375  A.newsize(M,N);
376  }
377 
378 
379  for (Subscript i=0; i<M; i++)
380  for (Subscript j=0; j<N; j++)
381  {
382  s >> A[i][j];
383  }
384 
385 
386  return s;
387 }
388 
389 // *******************[ basic matrix algorithms ]***************************
390 
391 
392 template <class T>
394  const Matrix<T> &B)
395 {
396  Subscript M = A.num_rows();
397  Subscript N = A.num_cols();
398 
399  assert(M==B.num_rows());
400  assert(N==B.num_cols());
401 
402  Matrix<T> tmp(M,N);
403  Subscript i,j;
404 
405  for (i=0; i<M; i++)
406  for (j=0; j<N; j++)
407  tmp[i][j] = A[i][j] + B[i][j];
408 
409  return tmp;
410 }
411 
412 template <class T>
414  const Matrix<T> &B)
415 {
416  Subscript M = A.num_rows();
417  Subscript N = A.num_cols();
418 
419  assert(M==B.num_rows());
420  assert(N==B.num_cols());
421 
422  Matrix<T> tmp(M,N);
423  Subscript i,j;
424 
425  for (i=0; i<M; i++)
426  for (j=0; j<N; j++)
427  tmp[i][j] = A[i][j] - B[i][j];
428 
429  return tmp;
430 }
431 
432 template <class T>
434  const Matrix<T> &B)
435 {
436  Subscript M = A.num_rows();
437  Subscript N = A.num_cols();
438 
439  assert(M==B.num_rows());
440  assert(N==B.num_cols());
441 
442  Matrix<T> tmp(M,N);
443  Subscript i,j;
444 
445  for (i=0; i<M; i++)
446  for (j=0; j<N; j++)
447  tmp[i][j] = A[i][j] * B[i][j];
448 
449  return tmp;
450 }
451 
452 
453 template <class T>
455 {
456  Subscript M = A.num_rows();
457  Subscript N = A.num_cols();
458 
459  Matrix<T> S(N,M);
460  Subscript i, j;
461 
462  for (i=0; i<M; i++)
463  for (j=0; j<N; j++)
464  S[j][i] = A[i][j];
465 
466  return S;
467 }
468 
469 
470 
471 template <class T>
472 inline Matrix<T> matmult(const Matrix<T> &A,
473  const Matrix<T> &B)
474 {
475 
476 #ifdef TNT_BOUNDS_CHECK
477  assert(A.num_cols() == B.num_rows());
478 #endif
479 
480  Subscript M = A.num_rows();
481  Subscript N = A.num_cols();
482  Subscript K = B.num_cols();
483 
484  Matrix<T> tmp(M,K);
485  T sum;
486 
487  for (Subscript i=0; i<M; i++)
488  for (Subscript k=0; k<K; k++)
489  {
490  sum = 0;
491  for (Subscript j=0; j<N; j++)
492  sum = sum + A[i][j] * B[j][k];
493 
494  tmp[i][k] = sum;
495  }
496 
497  return tmp;
498 }
499 
500 template <class T>
501 inline Matrix<T> operator*(const Matrix<T> &A,
502  const Matrix<T> &B)
503 {
504  return matmult(A,B);
505 }
506 
507 template <class T>
508 inline int matmult(Matrix<T>& C, const Matrix<T> &A,
509  const Matrix<T> &B)
510 {
511 
512  assert(A.num_cols() == B.num_rows());
513 
514  Subscript M = A.num_rows();
515  Subscript N = A.num_cols();
516  Subscript K = B.num_cols();
517 
518  C.newsize(M,K);
519 
520  T sum;
521 
522  const T* row_i;
523  const T* col_k;
524 
525  for (Subscript i=0; i<M; i++)
526  for (Subscript k=0; k<K; k++)
527  {
528  row_i = &(A[i][0]);
529  col_k = &(B[0][k]);
530  sum = 0;
531  for (Subscript j=0; j<N; j++)
532  {
533  sum += *row_i * *col_k;
534  row_i++;
535  col_k += K;
536  }
537  C[i][k] = sum;
538  }
539 
540  return 0;
541 }
542 
543 
544 template <class T>
546 {
547 
548 #ifdef TNT_BOUNDS_CHECK
549  assert(A.num_cols() == x.dim());
550 #endif
551 
552  Subscript M = A.num_rows();
553  Subscript N = A.num_cols();
554 
555  Vector<T> tmp(M);
556  T sum;
557 
558  for (Subscript i=0; i<M; i++)
559  {
560  sum = 0;
561  const T* rowi = A[i];
562  for (Subscript j=0; j<N; j++)
563  sum = sum + rowi[j] * x[j];
564 
565  tmp[i] = sum;
566  }
567 
568  return tmp;
569 }
570 
571 template <class T>
572 inline Vector<T> operator*(const Matrix<T> &A, const Vector<T> &x)
573 {
574  return matmult(A,x);
575 }
576 
577 } // namespace TNT
578 
579 #endif
580 // CMAT_H
Matrix< T > transpose(const Matrix< T > &A)
Definition: tnt_cmat.h:454
Definition: tnt_cmat.h:39
T sum(const std::vector< T > &array, int start, int end)
Definition: essentiamath.h:116
T * pointer
Definition: tnt_cmat.h:48
Subscript dim() const
Definition: tnt_vec.h:236
Subscript lbound() const
Definition: tnt_cmat.h:54
Matrix(Subscript M, Subscript N, const T *v)
Definition: tnt_cmat.h:189
reference operator()(Subscript i)
Definition: tnt_cmat.h:292
Subscript n_
Definition: tnt_cmat.h:58
T element_type
Definition: tnt_cmat.h:47
Definition: tnt_array1d.h:35
T value_type
Definition: tnt_cmat.h:46
Subscript dim(Subscript d) const
Definition: tnt_cmat.h:259
Matrix(Subscript M, Subscript N, const T &value=T())
Definition: tnt_cmat.h:183
void initialize(Subscript M, Subscript N)
Definition: tnt_cmat.h:68
Subscript m_
Definition: tnt_cmat.h:57
int Subscript
Definition: tnt_subscript.h:43
Matrix< T > mult_element(const Matrix< T > &A, const Matrix< T > &B)
Definition: tnt_cmat.h:433
const T * operator[](Subscript i) const
Definition: tnt_cmat.h:283
T * operator[](Subscript i)
Definition: tnt_cmat.h:274
Matrix()
Definition: tnt_cmat.h:175
#define NULL
Definition: tnt_i_refvec.h:33
Subscript size_type
Definition: tnt_cmat.h:45
void copy(const T *v)
Definition: tnt_cmat.h:95
T ** rowm1_
Definition: tnt_cmat.h:63
T * iterator
Definition: tnt_cmat.h:49
const_reference operator()(Subscript i) const
Definition: tnt_cmat.h:301
Matrix< T > & operator=(const Matrix< T > &A)
Definition: tnt_cmat.h:234
Matrix(Subscript M, Subscript N, const char *s)
Definition: tnt_cmat.h:195
Matrix(const Matrix< T > &A)
Definition: tnt_cmat.h:177
Array1D< T > operator+(const Array1D< T > &A, const Array1D< T > &B)
Definition: tnt_array1d_utils.h:64
Subscript size() const
Definition: tnt_cmat.h:171
Subscript num_rows() const
Definition: tnt_cmat.h:268
Subscript mn_
Definition: tnt_cmat.h:59
Matrix< T > & operator=(const T &scalar)
Definition: tnt_cmat.h:252
void destroy()
Definition: tnt_cmat.h:150
Array2D< T > matmult(const Array2D< T > &A, const Array2D< T > &B)
Definition: tnt_array2d_utils.h:259
const T * const_iterator
Definition: tnt_cmat.h:51
T & reference
Definition: tnt_cmat.h:50
T * v_
Definition: tnt_cmat.h:60
Array1D< T > operator*(const Array1D< T > &A, const Array1D< T > &B)
Definition: tnt_array1d_utils.h:107
Array1D< T > operator-(const Array1D< T > &A, const Array1D< T > &B)
Definition: tnt_array1d_utils.h:86
const T & const_reference
Definition: tnt_cmat.h:52
T * vm1_
Definition: tnt_cmat.h:62
reference operator()(Subscript i, Subscript j)
Definition: tnt_cmat.h:312
std::istream & operator>>(std::istream &s, Array1D< T > &A)
Definition: tnt_array1d_utils.h:49
~Matrix()
Definition: tnt_cmat.h:210
T ** row_
Definition: tnt_cmat.h:61
Matrix< T > & newsize(Subscript M, Subscript N)
Definition: tnt_cmat.h:218
Definition: tnt_vec.h:42
Subscript num_cols() const
Definition: tnt_cmat.h:269