/*!
* \file
* \brief Matrix Class Implementation
* \author Tony Ottosson, Tobias Ringstrom, Adam Piatyszek and Conrad Sanderson
*
* -------------------------------------------------------------------------
*
* Copyright (C) 1995-2010 (see AUTHORS file for a list of contributors)
*
* This file is part of IT++ - a C++ library of mathematical, signal
* processing, speech processing, and communications classes and functions.
*
* IT++ is free software: you can redistribute it and/or modify it under the
* terms of the GNU General Public License as published by the Free Software
* Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* IT++ is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
* FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
* details.
*
* You should have received a copy of the GNU General Public License along
* with IT++. If not, see <http://www.gnu.org/licenses/>.
*
* -------------------------------------------------------------------------
*/
#include <itpp/base/mat.h>
#ifndef _MSC_VER
# include <itpp/config.h>
#else
# include <itpp/config_msvc.h>
#endif
#if defined (HAVE_BLAS)
# include <itpp/base/blas.h>
#endif
//! \cond
namespace itpp
{
template<>
cmat cmat::hermitian_transpose() const
{
cmat temp(no_cols, no_rows);
for (int i = 0; i < no_rows; i++)
for (int j = 0; j < no_cols; j++)
temp(j, i) = std::conj(operator()(i,j));
return temp;
}
// -------- Multiplication operator -------------
#if defined(HAVE_BLAS)
template<>
mat& mat::operator*=(const mat &m)
{
it_assert_debug(no_cols == m.no_rows, "mat::operator*=(): Wrong sizes");
mat r(no_rows, m.no_cols); // unnecessary memory??
double alpha = 1.0;
double beta = 0.0;
char trans = 'n';
blas::dgemm_(&trans, &trans, &no_rows, &m.no_cols, &no_cols, &alpha, data,
&no_rows, m.data, &m.no_rows, &beta, r.data, &r.no_rows);
operator=(r); // time consuming
return *this;
}
template<>
cmat& cmat::operator*=(const cmat &m)
{
it_assert_debug(no_cols == m.no_rows, "cmat::operator*=(): Wrong sizes");
cmat r(no_rows, m.no_cols); // unnecessary memory??
std::complex<double> alpha = std::complex<double>(1.0);
std::complex<double> beta = std::complex<double>(0.0);
char trans = 'n';
blas::zgemm_(&trans, &trans, &no_rows, &m.no_cols, &no_cols, &alpha, data,
&no_rows, m.data, &m.no_rows, &beta, r.data, &r.no_rows);
operator=(r); // time consuming
return *this;
}
#else
template<>
mat& mat::operator*=(const mat &m)
{
it_assert_debug(no_cols == m.no_rows, "Mat<>::operator*=(): Wrong sizes");
mat r(no_rows, m.no_cols);
int r_pos = 0, pos = 0, m_pos = 0;
for (int i = 0; i < r.no_cols; i++) {
for (int j = 0; j < r.no_rows; j++) {
double tmp = 0.0;
pos = 0;
for (int k = 0; k < no_cols; k++) {
tmp += data[pos+j] * m.data[m_pos+k];
pos += no_rows;
}
r.data[r_pos+j] = tmp;
}
r_pos += r.no_rows;
m_pos += m.no_rows;
}
operator=(r); // time consuming
return *this;
}
template<>
cmat& cmat::operator*=(const cmat &m)
{
it_assert_debug(no_cols == m.no_rows, "Mat<>::operator*=(): Wrong sizes");
cmat r(no_rows, m.no_cols);
int r_pos = 0, pos = 0, m_pos = 0;
for (int i = 0; i < r.no_cols; i++) {
for (int j = 0; j < r.no_rows; j++) {
std::complex<double> tmp(0.0);
pos = 0;
for (int k = 0; k < no_cols; k++) {
tmp += data[pos+j] * m.data[m_pos+k];
pos += no_rows;
}
r.data[r_pos+j] = tmp;
}
r_pos += r.no_rows;
m_pos += m.no_rows;
}
operator=(r); // time consuming
return *this;
}
#endif // HAVE_BLAS
#if defined(HAVE_BLAS)
template<>
mat operator*(const mat &m1, const mat &m2)
{
it_assert_debug(m1.cols() == m2.rows(), "mat::operator*(): Wrong sizes");
int m1_r = m1.rows(); int m1_c = m1.cols();
int m2_r = m2.rows(); int m2_c = m2.cols();
mat r(m1_r, m2_c);
double alpha = 1.0;
double beta = 0.0;
char trans = 'n';
blas::dgemm_(&trans, &trans, &m1_r, &m2_c, &m1_c, &alpha,
m1._data(), &m1_r, m2._data(), &m2_r, &beta, r._data(),
&m1_r);
return r;
}
template<>
cmat operator*(const cmat &m1, const cmat &m2)
{
it_assert_debug(m1.cols() == m2.rows(), "cmat::operator*(): Wrong sizes");
int m1_r = m1.rows(); int m1_c = m1.cols();
int m2_r = m2.rows(); int m2_c = m2.cols();
cmat r(m1_r, m2_c);
std::complex<double> alpha = std::complex<double>(1.0);
std::complex<double> beta = std::complex<double>(0.0);
char trans = 'n';
blas::zgemm_(&trans, &trans, &m1_r, &m2_c, &m1_c, &alpha,
m1._data(), &m1_r, m2._data(), &m2_r, &beta, r._data(),
&m1_r);
return r;
}
#else
template<>
mat operator*(const mat &m1, const mat &m2)
{
it_assert_debug(m1.rows() == m2.cols(),
"Mat<>::operator*(): Wrong sizes");
mat r(m1.rows(), m2.cols());
double *tr = r._data();
const double *t1;
const double *t2 = m2._data();
for (int i = 0; i < r.cols(); i++) {
for (int j = 0; j < r.rows(); j++) {
double tmp = 0.0;
t1 = m1._data() + j;
for (int k = m1.cols(); k > 0; k--) {
tmp += *(t1) * *(t2++);
t1 += m1.rows();
}
*(tr++) = tmp;
t2 -= m2.rows();
}
t2 += m2.rows();
}
return r;
}
template<>
cmat operator*(const cmat &m1, const cmat &m2)
{
it_assert_debug(m1.cols() == m2.rows(),
"Mat<>::operator*(): Wrong sizes");
cmat r(m1.rows(), m2.cols());
std::complex<double> *tr = r._data();
const std::complex<double> *t1;
const std::complex<double> *t2 = m2._data();
for (int i = 0; i < r.cols(); i++) {
for (int j = 0; j < r.rows(); j++) {
std::complex<double> tmp(0.0);
t1 = m1._data() + j;
for (int k = m1.cols(); k > 0; k--) {
tmp += *(t1) * *(t2++);
t1 += m1.rows();
}
*(tr++) = tmp;
t2 -= m2.rows();
}
t2 += m2.rows();
}
return r;
}
#endif // HAVE_BLAS
#if defined(HAVE_BLAS)
template<>
vec operator*(const mat &m, const vec &v)
{
it_assert_debug(m.cols() == v.size(), "mat::operator*(): Wrong sizes");
int m_r = m.rows(); int m_c = m.cols();
vec r(m_r);
double alpha = 1.0;
double beta = 0.0;
char trans = 'n';
int incr = 1;
blas::dgemv_(&trans, &m_r, &m_c, &alpha, m._data(), &m_r,
v._data(), &incr, &beta, r._data(), &incr);
return r;
}
template<>
cvec operator*(const cmat &m, const cvec &v)
{
it_assert_debug(m.cols() == v.size(), "cmat::operator*(): Wrong sizes");
int m_r = m.rows(); int m_c = m.cols();
cvec r(m_r);
std::complex<double> alpha = std::complex<double>(1.0);
std::complex<double> beta = std::complex<double>(0.0);
char trans = 'n';
int incr = 1;
blas::zgemv_(&trans, &m_r, &m_c, &alpha, m._data(), &m_r,
v._data(), &incr, &beta, r._data(), &incr);
return r;
}
#else
template<>
vec operator*(const mat &m, const vec &v)
{
it_assert_debug(m.cols() == v.size(),
"Mat<>::operator*(): Wrong sizes");
vec r(m.rows());
for (int i = 0; i < m.rows(); i++) {
r(i) = 0.0;
int m_pos = 0;
for (int k = 0; k < m.cols(); k++) {
r(i) += m._data()[m_pos+i] * v(k);
m_pos += m.rows();
}
}
return r;
}
template<>
cvec operator*(const cmat &m, const cvec &v)
{
it_assert_debug(m.cols() == v.size(),
"Mat<>::operator*(): Wrong sizes");
cvec r(m.rows());
for (int i = 0; i < m.rows(); i++) {
r(i) = std::complex<double>(0.0);
int m_pos = 0;
for (int k = 0; k < m.cols(); k++) {
r(i) += m._data()[m_pos+i] * v(k);
m_pos += m.rows();
}
}
return r;
}
#endif // HAVE_BLAS
//---------------------------------------------------------------------
// Instantiations
//---------------------------------------------------------------------
// class instantiations
template class ITPP_EXPORT Mat<double>;
template class ITPP_EXPORT Mat<std::complex<double> >;
template class ITPP_EXPORT Mat<int>;
template class ITPP_EXPORT Mat<short int>;
template class ITPP_EXPORT Mat<bin>;
// addition operators
template ITPP_EXPORT mat operator+(const mat &m1, const mat &m2);
template ITPP_EXPORT cmat operator+(const cmat &m1, const cmat &m2);
template ITPP_EXPORT imat operator+(const imat &m1, const imat &m2);
template ITPP_EXPORT smat operator+(const smat &m1, const smat &m2);
template ITPP_EXPORT bmat operator+(const bmat &m1, const bmat &m2);
template ITPP_EXPORT mat operator+(const mat &m, double t);
template ITPP_EXPORT cmat operator+(const cmat &m, std::complex<double> t);
template ITPP_EXPORT imat operator+(const imat &m, int t);
template ITPP_EXPORT smat operator+(const smat &m, short t);
template ITPP_EXPORT bmat operator+(const bmat &m, bin t);
template ITPP_EXPORT mat operator+(double t, const mat &m);
template ITPP_EXPORT cmat operator+(std::complex<double> t, const cmat &m);
template ITPP_EXPORT imat operator+(int t, const imat &m);
template ITPP_EXPORT smat operator+(short t, const smat &m);
template ITPP_EXPORT bmat operator+(bin t, const bmat &m);
// subraction operators
template ITPP_EXPORT mat operator-(const mat &m1, const mat &m2);
template ITPP_EXPORT cmat operator-(const cmat &m1, const cmat &m2);
template ITPP_EXPORT imat operator-(const imat &m1, const imat &m2);
template ITPP_EXPORT smat operator-(const smat &m1, const smat &m2);
template ITPP_EXPORT bmat operator-(const bmat &m1, const bmat &m2);
template ITPP_EXPORT mat operator-(const mat &m, double t);
template ITPP_EXPORT cmat operator-(const cmat &m, std::complex<double> t);
template ITPP_EXPORT imat operator-(const imat &m, int t);
template ITPP_EXPORT smat operator-(const smat &m, short t);
template ITPP_EXPORT bmat operator-(const bmat &m, bin t);
template ITPP_EXPORT mat operator-(double t, const mat &m);
template ITPP_EXPORT cmat operator-(std::complex<double> t, const cmat &m);
template ITPP_EXPORT imat operator-(int t, const imat &m);
template ITPP_EXPORT smat operator-(short t, const smat &m);
template ITPP_EXPORT bmat operator-(bin t, const bmat &m);
// unary minus
template ITPP_EXPORT mat operator-(const mat &m);
template ITPP_EXPORT cmat operator-(const cmat &m);
template ITPP_EXPORT imat operator-(const imat &m);
template ITPP_EXPORT smat operator-(const smat &m);
template ITPP_EXPORT bmat operator-(const bmat &m);
// multiplication operators
template ITPP_EXPORT imat operator*(const imat &m1, const imat &m2);
template ITPP_EXPORT smat operator*(const smat &m1, const smat &m2);
template ITPP_EXPORT bmat operator*(const bmat &m1, const bmat &m2);
template ITPP_EXPORT ivec operator*(const imat &m, const ivec &v);
template ITPP_EXPORT svec operator*(const smat &m, const svec &v);
template ITPP_EXPORT bvec operator*(const bmat &m, const bvec &v);
template ITPP_EXPORT mat operator*(const mat &m, double t);
template ITPP_EXPORT cmat operator*(const cmat &m, std::complex<double> t);
template ITPP_EXPORT imat operator*(const imat &m, int t);
template ITPP_EXPORT smat operator*(const smat &m, short t);
template ITPP_EXPORT bmat operator*(const bmat &m, bin t);
template ITPP_EXPORT mat operator*(double t, const mat &m);
template ITPP_EXPORT cmat operator*(std::complex<double> t, const cmat &m);
template ITPP_EXPORT imat operator*(int t, const imat &m);
template ITPP_EXPORT smat operator*(short t, const smat &m);
template ITPP_EXPORT bmat operator*(bin t, const bmat &m);
// elementwise multiplication
template ITPP_EXPORT mat elem_mult(const mat &m1, const mat &m2);
template ITPP_EXPORT cmat elem_mult(const cmat &m1, const cmat &m2);
template ITPP_EXPORT imat elem_mult(const imat &m1, const imat &m2);
template ITPP_EXPORT smat elem_mult(const smat &m1, const smat &m2);
template ITPP_EXPORT bmat elem_mult(const bmat &m1, const bmat &m2);
template ITPP_EXPORT void elem_mult_out(const mat &m1, const mat &m2, mat &out);
template ITPP_EXPORT void elem_mult_out(const cmat &m1, const cmat &m2, cmat &out);
template ITPP_EXPORT void elem_mult_out(const imat &m1, const imat &m2, imat &out);
template ITPP_EXPORT void elem_mult_out(const smat &m1, const smat &m2, smat &out);
template ITPP_EXPORT void elem_mult_out(const bmat &m1, const bmat &m2, bmat &out);
template ITPP_EXPORT void elem_mult_out(const mat &m1, const mat &m2,
const mat &m3, mat &out);
template ITPP_EXPORT void elem_mult_out(const cmat &m1, const cmat &m2,
const cmat &m3, cmat &out);
template ITPP_EXPORT void elem_mult_out(const imat &m1, const imat &m2,
const imat &m3, imat &out);
template ITPP_EXPORT void elem_mult_out(const smat &m1, const smat &m2,
const smat &m3, smat &out);
template ITPP_EXPORT void elem_mult_out(const bmat &m1, const bmat &m2,
const bmat &m3, bmat &out);
template ITPP_EXPORT void elem_mult_out(const mat &m1, const mat &m2, const mat &m3,
const mat &m4, mat &out);
template ITPP_EXPORT void elem_mult_out(const cmat &m1, const cmat &m2,
const cmat &m3, const cmat &m4, cmat &out);
template ITPP_EXPORT void elem_mult_out(const imat &m1, const imat &m2,
const imat &m3, const imat &m4, imat &out);
template ITPP_EXPORT void elem_mult_out(const smat &m1, const smat &m2,
const smat &m3, const smat &m4, smat &out);
template ITPP_EXPORT void elem_mult_out(const bmat &m1, const bmat &m2,
const bmat &m3, const bmat &m4, bmat &out);
template ITPP_EXPORT void elem_mult_inplace(const mat &m1, mat &m2);
template ITPP_EXPORT void elem_mult_inplace(const cmat &m1, cmat &m2);
template ITPP_EXPORT void elem_mult_inplace(const imat &m1, imat &m2);
template ITPP_EXPORT void elem_mult_inplace(const smat &m1, smat &m2);
template ITPP_EXPORT void elem_mult_inplace(const bmat &m1, bmat &m2);
template ITPP_EXPORT double elem_mult_sum(const mat &m1, const mat &m2);
template ITPP_EXPORT std::complex<double> elem_mult_sum(const cmat &m1, const cmat &m2);
template ITPP_EXPORT int elem_mult_sum(const imat &m1, const imat &m2);
template ITPP_EXPORT short elem_mult_sum(const smat &m1, const smat &m2);
template ITPP_EXPORT bin elem_mult_sum(const bmat &m1, const bmat &m2);
// division operator
template ITPP_EXPORT mat operator/(double t, const mat &m);
template ITPP_EXPORT cmat operator/(std::complex<double> t, const cmat &m);
template ITPP_EXPORT imat operator/(int t, const imat &m);
template ITPP_EXPORT smat operator/(short t, const smat &m);
template ITPP_EXPORT bmat operator/(bin t, const bmat &m);
template ITPP_EXPORT mat operator/(const mat &m, double t);
template ITPP_EXPORT cmat operator/(const cmat &m, std::complex<double> t);
template ITPP_EXPORT imat operator/(const imat &m, int t);
template ITPP_EXPORT smat operator/(const smat &m, short t);
template ITPP_EXPORT bmat operator/(const bmat &m, bin t);
// elementwise division
template ITPP_EXPORT mat elem_div(const mat &m1, const mat &m2);
template ITPP_EXPORT cmat elem_div(const cmat &m1, const cmat &m2);
template ITPP_EXPORT imat elem_div(const imat &m1, const imat &m2);
template ITPP_EXPORT smat elem_div(const smat &m1, const smat &m2);
template ITPP_EXPORT bmat elem_div(const bmat &m1, const bmat &m2);
template ITPP_EXPORT void elem_div_out(const mat &m1, const mat &m2, mat &out);
template ITPP_EXPORT void elem_div_out(const cmat &m1, const cmat &m2, cmat &out);
template ITPP_EXPORT void elem_div_out(const imat &m1, const imat &m2, imat &out);
template ITPP_EXPORT void elem_div_out(const smat &m1, const smat &m2, smat &out);
template ITPP_EXPORT void elem_div_out(const bmat &m1, const bmat &m2, bmat &out);
template ITPP_EXPORT double elem_div_sum(const mat &m1, const mat &m2);
template ITPP_EXPORT std::complex<double> elem_div_sum(const cmat &m1,
const cmat &m2);
template ITPP_EXPORT int elem_div_sum(const imat &m1, const imat &m2);
template ITPP_EXPORT short elem_div_sum(const smat &m1, const smat &m2);
template ITPP_EXPORT bin elem_div_sum(const bmat &m1, const bmat &m2);
// concatenation
template ITPP_EXPORT mat concat_horizontal(const mat &m1, const mat &m2);
template ITPP_EXPORT cmat concat_horizontal(const cmat &m1, const cmat &m2);
template ITPP_EXPORT imat concat_horizontal(const imat &m1, const imat &m2);
template ITPP_EXPORT smat concat_horizontal(const smat &m1, const smat &m2);
template ITPP_EXPORT bmat concat_horizontal(const bmat &m1, const bmat &m2);
template ITPP_EXPORT mat concat_vertical(const mat &m1, const mat &m2);
template ITPP_EXPORT cmat concat_vertical(const cmat &m1, const cmat &m2);
template ITPP_EXPORT imat concat_vertical(const imat &m1, const imat &m2);
template ITPP_EXPORT smat concat_vertical(const smat &m1, const smat &m2);
template ITPP_EXPORT bmat concat_vertical(const bmat &m1, const bmat &m2);
// I/O streams
template ITPP_EXPORT std::ostream &operator<<(std::ostream &os, const mat &m);
template ITPP_EXPORT std::ostream &operator<<(std::ostream &os, const cmat &m);
template ITPP_EXPORT std::ostream &operator<<(std::ostream &os, const imat &m);
template ITPP_EXPORT std::ostream &operator<<(std::ostream &os, const smat &m);
template ITPP_EXPORT std::ostream &operator<<(std::ostream &os, const bmat &m);
template ITPP_EXPORT std::istream &operator>>(std::istream &is, mat &m);
template ITPP_EXPORT std::istream &operator>>(std::istream &is, cmat &m);
template ITPP_EXPORT std::istream &operator>>(std::istream &is, imat &m);
template ITPP_EXPORT std::istream &operator>>(std::istream &is, smat &m);
template ITPP_EXPORT std::istream &operator>>(std::istream &is, bmat &m);
} // namespace itpp
//! \endcond