IT++ Logo
mat.cpp
Go to the documentation of this file.
1 
29 #include <itpp/base/mat.h>
30 
31 #ifndef _MSC_VER
32 # include <itpp/config.h>
33 #else
34 # include <itpp/config_msvc.h>
35 #endif
36 
37 #if defined (HAVE_BLAS)
38 # include <itpp/base/blas.h>
39 #endif
40 
42 
43 namespace itpp
44 {
45 
46 template<>
47 cmat cmat::hermitian_transpose() const
48 {
49  cmat temp(no_cols, no_rows);
50  for (int i = 0; i < no_rows; i++)
51  for (int j = 0; j < no_cols; j++)
52  temp(j, i) = std::conj(operator()(i,j));
53 
54  return temp;
55 }
56 
57 
58 // -------- Multiplication operator -------------
59 
60 #if defined(HAVE_BLAS)
61 
62 template<>
63 mat& mat::operator*=(const mat &m)
64 {
65  it_assert_debug(no_cols == m.no_rows, "mat::operator*=(): Wrong sizes");
66  mat r(no_rows, m.no_cols); // unnecessary memory??
67  double alpha = 1.0;
68  double beta = 0.0;
69  char trans = 'n';
70  blas::dgemm_(&trans, &trans, &no_rows, &m.no_cols, &no_cols, &alpha, data,
71  &no_rows, m.data, &m.no_rows, &beta, r.data, &r.no_rows);
72  operator=(r); // time consuming
73  return *this;
74 }
75 
76 template<>
77 cmat& cmat::operator*=(const cmat &m)
78 {
79  it_assert_debug(no_cols == m.no_rows, "cmat::operator*=(): Wrong sizes");
80  cmat r(no_rows, m.no_cols); // unnecessary memory??
81  std::complex<double> alpha = std::complex<double>(1.0);
82  std::complex<double> beta = std::complex<double>(0.0);
83  char trans = 'n';
84  blas::zgemm_(&trans, &trans, &no_rows, &m.no_cols, &no_cols, &alpha, data,
85  &no_rows, m.data, &m.no_rows, &beta, r.data, &r.no_rows);
86  operator=(r); // time consuming
87  return *this;
88 }
89 #else
90 template<>
91 mat& mat::operator*=(const mat &m)
92 {
93  it_assert_debug(no_cols == m.no_rows, "Mat<>::operator*=(): Wrong sizes");
94  mat r(no_rows, m.no_cols);
95  int r_pos = 0, pos = 0, m_pos = 0;
96 
97  for (int i = 0; i < r.no_cols; i++) {
98  for (int j = 0; j < r.no_rows; j++) {
99  double tmp = 0.0;
100  pos = 0;
101  for (int k = 0; k < no_cols; k++) {
102  tmp += data[pos+j] * m.data[m_pos+k];
103  pos += no_rows;
104  }
105  r.data[r_pos+j] = tmp;
106  }
107  r_pos += r.no_rows;
108  m_pos += m.no_rows;
109  }
110  operator=(r); // time consuming
111  return *this;
112 }
113 
114 template<>
115 cmat& cmat::operator*=(const cmat &m)
116 {
117  it_assert_debug(no_cols == m.no_rows, "Mat<>::operator*=(): Wrong sizes");
118  cmat r(no_rows, m.no_cols);
119  int r_pos = 0, pos = 0, m_pos = 0;
120 
121  for (int i = 0; i < r.no_cols; i++) {
122  for (int j = 0; j < r.no_rows; j++) {
123  std::complex<double> tmp(0.0);
124  pos = 0;
125  for (int k = 0; k < no_cols; k++) {
126  tmp += data[pos+j] * m.data[m_pos+k];
127  pos += no_rows;
128  }
129  r.data[r_pos+j] = tmp;
130  }
131  r_pos += r.no_rows;
132  m_pos += m.no_rows;
133  }
134  operator=(r); // time consuming
135  return *this;
136 }
137 #endif // HAVE_BLAS
138 
139 
140 #if defined(HAVE_BLAS)
141 template<>
142 mat operator*(const mat &m1, const mat &m2)
143 {
144  it_assert_debug(m1.cols() == m2.rows(), "mat::operator*(): Wrong sizes");
145  int m1_r = m1.rows(); int m1_c = m1.cols();
146  int m2_r = m2.rows(); int m2_c = m2.cols();
147  mat r(m1_r, m2_c);
148  double alpha = 1.0;
149  double beta = 0.0;
150  char trans = 'n';
151  blas::dgemm_(&trans, &trans, &m1_r, &m2_c, &m1_c, &alpha,
152  m1._data(), &m1_r, m2._data(), &m2_r, &beta, r._data(),
153  &m1_r);
154  return r;
155 }
156 
157 template<>
158 cmat operator*(const cmat &m1, const cmat &m2)
159 {
160  it_assert_debug(m1.cols() == m2.rows(), "cmat::operator*(): Wrong sizes");
161  int m1_r = m1.rows(); int m1_c = m1.cols();
162  int m2_r = m2.rows(); int m2_c = m2.cols();
163  cmat r(m1_r, m2_c);
164  std::complex<double> alpha = std::complex<double>(1.0);
165  std::complex<double> beta = std::complex<double>(0.0);
166  char trans = 'n';
167  blas::zgemm_(&trans, &trans, &m1_r, &m2_c, &m1_c, &alpha,
168  m1._data(), &m1_r, m2._data(), &m2_r, &beta, r._data(),
169  &m1_r);
170  return r;
171 }
172 #else
173 template<>
174 mat operator*(const mat &m1, const mat &m2)
175 {
176  it_assert_debug(m1.rows() == m2.cols(),
177  "Mat<>::operator*(): Wrong sizes");
178  mat r(m1.rows(), m2.cols());
179  double *tr = r._data();
180  const double *t1;
181  const double *t2 = m2._data();
182  for (int i = 0; i < r.cols(); i++) {
183  for (int j = 0; j < r.rows(); j++) {
184  double tmp = 0.0;
185  t1 = m1._data() + j;
186  for (int k = m1.cols(); k > 0; k--) {
187  tmp += *(t1) * *(t2++);
188  t1 += m1.rows();
189  }
190  *(tr++) = tmp;
191  t2 -= m2.rows();
192  }
193  t2 += m2.rows();
194  }
195  return r;
196 }
197 
198 template<>
199 cmat operator*(const cmat &m1, const cmat &m2)
200 {
201  it_assert_debug(m1.cols() == m2.rows(),
202  "Mat<>::operator*(): Wrong sizes");
203  cmat r(m1.rows(), m2.cols());
204  std::complex<double> *tr = r._data();
205  const std::complex<double> *t1;
206  const std::complex<double> *t2 = m2._data();
207  for (int i = 0; i < r.cols(); i++) {
208  for (int j = 0; j < r.rows(); j++) {
209  std::complex<double> tmp(0.0);
210  t1 = m1._data() + j;
211  for (int k = m1.cols(); k > 0; k--) {
212  tmp += *(t1) * *(t2++);
213  t1 += m1.rows();
214  }
215  *(tr++) = tmp;
216  t2 -= m2.rows();
217  }
218  t2 += m2.rows();
219  }
220  return r;
221 }
222 #endif // HAVE_BLAS
223 
224 
225 #if defined(HAVE_BLAS)
226 template<>
227 vec operator*(const mat &m, const vec &v)
228 {
229  it_assert_debug(m.cols() == v.size(), "mat::operator*(): Wrong sizes");
230  int m_r = m.rows(); int m_c = m.cols();
231  vec r(m_r);
232  double alpha = 1.0;
233  double beta = 0.0;
234  char trans = 'n';
235  int incr = 1;
236  blas::dgemv_(&trans, &m_r, &m_c, &alpha, m._data(), &m_r,
237  v._data(), &incr, &beta, r._data(), &incr);
238  return r;
239 }
240 
241 template<>
242 cvec operator*(const cmat &m, const cvec &v)
243 {
244  it_assert_debug(m.cols() == v.size(), "cmat::operator*(): Wrong sizes");
245  int m_r = m.rows(); int m_c = m.cols();
246  cvec r(m_r);
247  std::complex<double> alpha = std::complex<double>(1.0);
248  std::complex<double> beta = std::complex<double>(0.0);
249  char trans = 'n';
250  int incr = 1;
251  blas::zgemv_(&trans, &m_r, &m_c, &alpha, m._data(), &m_r,
252  v._data(), &incr, &beta, r._data(), &incr);
253  return r;
254 }
255 #else
256 template<>
257 vec operator*(const mat &m, const vec &v)
258 {
259  it_assert_debug(m.cols() == v.size(),
260  "Mat<>::operator*(): Wrong sizes");
261  vec r(m.rows());
262  for (int i = 0; i < m.rows(); i++) {
263  r(i) = 0.0;
264  int m_pos = 0;
265  for (int k = 0; k < m.cols(); k++) {
266  r(i) += m._data()[m_pos+i] * v(k);
267  m_pos += m.rows();
268  }
269  }
270  return r;
271 }
272 
273 template<>
274 cvec operator*(const cmat &m, const cvec &v)
275 {
276  it_assert_debug(m.cols() == v.size(),
277  "Mat<>::operator*(): Wrong sizes");
278  cvec r(m.rows());
279  for (int i = 0; i < m.rows(); i++) {
280  r(i) = std::complex<double>(0.0);
281  int m_pos = 0;
282  for (int k = 0; k < m.cols(); k++) {
283  r(i) += m._data()[m_pos+i] * v(k);
284  m_pos += m.rows();
285  }
286  }
287  return r;
288 }
289 #endif // HAVE_BLAS
290 
291 
292 //---------------------------------------------------------------------
293 // Instantiations
294 //---------------------------------------------------------------------
295 
296 // class instantiations
297 
298 template class ITPP_EXPORT Mat<double>;
299 template class ITPP_EXPORT Mat<std::complex<double> >;
300 template class ITPP_EXPORT Mat<int>;
301 template class ITPP_EXPORT Mat<short int>;
302 template class ITPP_EXPORT Mat<bin>;
303 
304 // addition operators
305 
306 template ITPP_EXPORT mat operator+(const mat &m1, const mat &m2);
307 template ITPP_EXPORT cmat operator+(const cmat &m1, const cmat &m2);
308 template ITPP_EXPORT imat operator+(const imat &m1, const imat &m2);
309 template ITPP_EXPORT smat operator+(const smat &m1, const smat &m2);
310 template ITPP_EXPORT bmat operator+(const bmat &m1, const bmat &m2);
311 
312 template ITPP_EXPORT mat operator+(const mat &m, double t);
313 template ITPP_EXPORT cmat operator+(const cmat &m, std::complex<double> t);
314 template ITPP_EXPORT imat operator+(const imat &m, int t);
315 template ITPP_EXPORT smat operator+(const smat &m, short t);
316 template ITPP_EXPORT bmat operator+(const bmat &m, bin t);
317 
318 template ITPP_EXPORT mat operator+(double t, const mat &m);
319 template ITPP_EXPORT cmat operator+(std::complex<double> t, const cmat &m);
320 template ITPP_EXPORT imat operator+(int t, const imat &m);
321 template ITPP_EXPORT smat operator+(short t, const smat &m);
322 template ITPP_EXPORT bmat operator+(bin t, const bmat &m);
323 
324 // subraction operators
325 
326 template ITPP_EXPORT mat operator-(const mat &m1, const mat &m2);
327 template ITPP_EXPORT cmat operator-(const cmat &m1, const cmat &m2);
328 template ITPP_EXPORT imat operator-(const imat &m1, const imat &m2);
329 template ITPP_EXPORT smat operator-(const smat &m1, const smat &m2);
330 template ITPP_EXPORT bmat operator-(const bmat &m1, const bmat &m2);
331 
332 template ITPP_EXPORT mat operator-(const mat &m, double t);
333 template ITPP_EXPORT cmat operator-(const cmat &m, std::complex<double> t);
334 template ITPP_EXPORT imat operator-(const imat &m, int t);
335 template ITPP_EXPORT smat operator-(const smat &m, short t);
336 template ITPP_EXPORT bmat operator-(const bmat &m, bin t);
337 
338 template ITPP_EXPORT mat operator-(double t, const mat &m);
339 template ITPP_EXPORT cmat operator-(std::complex<double> t, const cmat &m);
340 template ITPP_EXPORT imat operator-(int t, const imat &m);
341 template ITPP_EXPORT smat operator-(short t, const smat &m);
342 template ITPP_EXPORT bmat operator-(bin t, const bmat &m);
343 
344 // unary minus
345 
346 template ITPP_EXPORT mat operator-(const mat &m);
347 template ITPP_EXPORT cmat operator-(const cmat &m);
348 template ITPP_EXPORT imat operator-(const imat &m);
349 template ITPP_EXPORT smat operator-(const smat &m);
350 template ITPP_EXPORT bmat operator-(const bmat &m);
351 
352 // multiplication operators
353 
354 template ITPP_EXPORT imat operator*(const imat &m1, const imat &m2);
355 template ITPP_EXPORT smat operator*(const smat &m1, const smat &m2);
356 template ITPP_EXPORT bmat operator*(const bmat &m1, const bmat &m2);
357 
358 template ITPP_EXPORT ivec operator*(const imat &m, const ivec &v);
359 template ITPP_EXPORT svec operator*(const smat &m, const svec &v);
360 template ITPP_EXPORT bvec operator*(const bmat &m, const bvec &v);
361 
362 template ITPP_EXPORT mat operator*(const mat &m, double t);
363 template ITPP_EXPORT cmat operator*(const cmat &m, std::complex<double> t);
364 template ITPP_EXPORT imat operator*(const imat &m, int t);
365 template ITPP_EXPORT smat operator*(const smat &m, short t);
366 template ITPP_EXPORT bmat operator*(const bmat &m, bin t);
367 
368 template ITPP_EXPORT mat operator*(double t, const mat &m);
369 template ITPP_EXPORT cmat operator*(std::complex<double> t, const cmat &m);
370 template ITPP_EXPORT imat operator*(int t, const imat &m);
371 template ITPP_EXPORT smat operator*(short t, const smat &m);
372 template ITPP_EXPORT bmat operator*(bin t, const bmat &m);
373 
374 // elementwise multiplication
375 
376 template ITPP_EXPORT mat elem_mult(const mat &m1, const mat &m2);
377 template ITPP_EXPORT cmat elem_mult(const cmat &m1, const cmat &m2);
378 template ITPP_EXPORT imat elem_mult(const imat &m1, const imat &m2);
379 template ITPP_EXPORT smat elem_mult(const smat &m1, const smat &m2);
380 template ITPP_EXPORT bmat elem_mult(const bmat &m1, const bmat &m2);
381 
382 template ITPP_EXPORT void elem_mult_out(const mat &m1, const mat &m2, mat &out);
383 template ITPP_EXPORT void elem_mult_out(const cmat &m1, const cmat &m2, cmat &out);
384 template ITPP_EXPORT void elem_mult_out(const imat &m1, const imat &m2, imat &out);
385 template ITPP_EXPORT void elem_mult_out(const smat &m1, const smat &m2, smat &out);
386 template ITPP_EXPORT void elem_mult_out(const bmat &m1, const bmat &m2, bmat &out);
387 
388 template ITPP_EXPORT void elem_mult_out(const mat &m1, const mat &m2,
389  const mat &m3, mat &out);
390 template ITPP_EXPORT void elem_mult_out(const cmat &m1, const cmat &m2,
391  const cmat &m3, cmat &out);
392 template ITPP_EXPORT void elem_mult_out(const imat &m1, const imat &m2,
393  const imat &m3, imat &out);
394 template ITPP_EXPORT void elem_mult_out(const smat &m1, const smat &m2,
395  const smat &m3, smat &out);
396 template ITPP_EXPORT void elem_mult_out(const bmat &m1, const bmat &m2,
397  const bmat &m3, bmat &out);
398 
399 template ITPP_EXPORT void elem_mult_out(const mat &m1, const mat &m2, const mat &m3,
400  const mat &m4, mat &out);
401 template ITPP_EXPORT void elem_mult_out(const cmat &m1, const cmat &m2,
402  const cmat &m3, const cmat &m4, cmat &out);
403 template ITPP_EXPORT void elem_mult_out(const imat &m1, const imat &m2,
404  const imat &m3, const imat &m4, imat &out);
405 template ITPP_EXPORT void elem_mult_out(const smat &m1, const smat &m2,
406  const smat &m3, const smat &m4, smat &out);
407 template ITPP_EXPORT void elem_mult_out(const bmat &m1, const bmat &m2,
408  const bmat &m3, const bmat &m4, bmat &out);
409 
410 template ITPP_EXPORT void elem_mult_inplace(const mat &m1, mat &m2);
411 template ITPP_EXPORT void elem_mult_inplace(const cmat &m1, cmat &m2);
412 template ITPP_EXPORT void elem_mult_inplace(const imat &m1, imat &m2);
413 template ITPP_EXPORT void elem_mult_inplace(const smat &m1, smat &m2);
414 template ITPP_EXPORT void elem_mult_inplace(const bmat &m1, bmat &m2);
415 
416 template ITPP_EXPORT double elem_mult_sum(const mat &m1, const mat &m2);
417 template ITPP_EXPORT std::complex<double> elem_mult_sum(const cmat &m1, const cmat &m2);
418 template ITPP_EXPORT int elem_mult_sum(const imat &m1, const imat &m2);
419 template ITPP_EXPORT short elem_mult_sum(const smat &m1, const smat &m2);
420 template ITPP_EXPORT bin elem_mult_sum(const bmat &m1, const bmat &m2);
421 
422 // division operator
423 
424 template ITPP_EXPORT mat operator/(double t, const mat &m);
425 template ITPP_EXPORT cmat operator/(std::complex<double> t, const cmat &m);
426 template ITPP_EXPORT imat operator/(int t, const imat &m);
427 template ITPP_EXPORT smat operator/(short t, const smat &m);
428 template ITPP_EXPORT bmat operator/(bin t, const bmat &m);
429 
430 template ITPP_EXPORT mat operator/(const mat &m, double t);
431 template ITPP_EXPORT cmat operator/(const cmat &m, std::complex<double> t);
432 template ITPP_EXPORT imat operator/(const imat &m, int t);
433 template ITPP_EXPORT smat operator/(const smat &m, short t);
434 template ITPP_EXPORT bmat operator/(const bmat &m, bin t);
435 
436 // elementwise division
437 
438 template ITPP_EXPORT mat elem_div(const mat &m1, const mat &m2);
439 template ITPP_EXPORT cmat elem_div(const cmat &m1, const cmat &m2);
440 template ITPP_EXPORT imat elem_div(const imat &m1, const imat &m2);
441 template ITPP_EXPORT smat elem_div(const smat &m1, const smat &m2);
442 template ITPP_EXPORT bmat elem_div(const bmat &m1, const bmat &m2);
443 
444 template ITPP_EXPORT void elem_div_out(const mat &m1, const mat &m2, mat &out);
445 template ITPP_EXPORT void elem_div_out(const cmat &m1, const cmat &m2, cmat &out);
446 template ITPP_EXPORT void elem_div_out(const imat &m1, const imat &m2, imat &out);
447 template ITPP_EXPORT void elem_div_out(const smat &m1, const smat &m2, smat &out);
448 template ITPP_EXPORT void elem_div_out(const bmat &m1, const bmat &m2, bmat &out);
449 
450 template ITPP_EXPORT double elem_div_sum(const mat &m1, const mat &m2);
451 template ITPP_EXPORT std::complex<double> elem_div_sum(const cmat &m1,
452  const cmat &m2);
453 template ITPP_EXPORT int elem_div_sum(const imat &m1, const imat &m2);
454 template ITPP_EXPORT short elem_div_sum(const smat &m1, const smat &m2);
455 template ITPP_EXPORT bin elem_div_sum(const bmat &m1, const bmat &m2);
456 
457 // concatenation
458 
459 template ITPP_EXPORT mat concat_horizontal(const mat &m1, const mat &m2);
460 template ITPP_EXPORT cmat concat_horizontal(const cmat &m1, const cmat &m2);
461 template ITPP_EXPORT imat concat_horizontal(const imat &m1, const imat &m2);
462 template ITPP_EXPORT smat concat_horizontal(const smat &m1, const smat &m2);
463 template ITPP_EXPORT bmat concat_horizontal(const bmat &m1, const bmat &m2);
464 
465 template ITPP_EXPORT mat concat_vertical(const mat &m1, const mat &m2);
466 template ITPP_EXPORT cmat concat_vertical(const cmat &m1, const cmat &m2);
467 template ITPP_EXPORT imat concat_vertical(const imat &m1, const imat &m2);
468 template ITPP_EXPORT smat concat_vertical(const smat &m1, const smat &m2);
469 template ITPP_EXPORT bmat concat_vertical(const bmat &m1, const bmat &m2);
470 
471 // I/O streams
472 
473 template ITPP_EXPORT std::ostream &operator<<(std::ostream &os, const mat &m);
474 template ITPP_EXPORT std::ostream &operator<<(std::ostream &os, const cmat &m);
475 template ITPP_EXPORT std::ostream &operator<<(std::ostream &os, const imat &m);
476 template ITPP_EXPORT std::ostream &operator<<(std::ostream &os, const smat &m);
477 template ITPP_EXPORT std::ostream &operator<<(std::ostream &os, const bmat &m);
478 
479 template ITPP_EXPORT std::istream &operator>>(std::istream &is, mat &m);
480 template ITPP_EXPORT std::istream &operator>>(std::istream &is, cmat &m);
481 template ITPP_EXPORT std::istream &operator>>(std::istream &is, imat &m);
482 template ITPP_EXPORT std::istream &operator>>(std::istream &is, smat &m);
483 template ITPP_EXPORT std::istream &operator>>(std::istream &is, bmat &m);
484 
485 } // namespace itpp
486 
SourceForge Logo

Generated on Sat Jul 6 2013 10:54:20 for IT++ by Doxygen 1.8.2