/*---------------------------------------------------------------------------*
* IT++ *
*---------------------------------------------------------------------------*
* Copyright (c) 1995-2004 by Tony Ottosson, Thomas Eriksson, Pål Frenger, *
* Tobias Ringström, and Jonas Samuelsson. *
* *
* Permission to use, copy, modify, and distribute this software and its *
* documentation under the terms of the GNU General Public License is hereby *
* granted. No representations are made about the suitability of this *
* software for any purpose. It is provided "as is" without expressed or *
* implied warranty. See the GNU General Public License for more details. *
*---------------------------------------------------------------------------*/
/*!
\file
\brief Sparse Matrix Class Definitions
\author Tony Ottosson and Tobias Ringström
1.11
2003/06/06 10:16:12
*/
#ifndef __smat_h
#define __smat_h
#include "itconfig.h"
#include "base/svec.h"
#include "base/vec.h"
#include "base/mat.h"
namespace itpp {
// Declaration of class Vec
template <class T> class Vec;
// Declaration of class Mat
template <class T> class Mat;
// Declaration of class Sparse_Vec
template <class T> class Sparse_Vec;
// Declaration of class Sparse_Mat
template <class T> class Sparse_Mat;
// ------------------------ Sparse_Mat Friends -------------------------------------
//! m1+m2 where m1 and m2 are sparse matrices
template <class T>
Sparse_Mat<T> operator+(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
//! m1*m2 where m1 and m2 are sparse matrices
template <class T>
Sparse_Mat<T> operator*(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
//! m*v where m is a sparse matrix and v is a full column vector
template <class T>
Vec<T> operator*(const Sparse_Mat<T> &m, const Vec<T> &v);
//! v'*m where m is a sparse matrix and v is a full column vector
template <class T>
Vec<T> operator*(const Vec<T> &v, const Sparse_Mat<T> &m);
//! m'*m where m is a sparse matrix
template <class T>
Mat<T> trans_mult(const Sparse_Mat<T> &m);
//! m'*m where m is a sparse matrix
template <class T>
Sparse_Mat<T> trans_mult_s(const Sparse_Mat<T> &m);
//! m1'*m2 where m1 and m2 are sparse matrices
template <class T>
Sparse_Mat<T> trans_mult(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
//! m'*v where m is a sparse matrix and v is a full column vector
template <class T>
Vec<T> trans_mult(const Sparse_Mat<T> &m, const Vec<T> &v);
//! m1*m2' where m1 and m2 are sparse matrices
template <class T>
Sparse_Mat<T> mult_trans(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
/*!
\brief Templated Sparse Matrix Class
\author Tony Ottosson and Tobias Ringstrom
*/
template <class T>
class Sparse_Mat {
public:
//! Default constructor
Sparse_Mat();
/*!
\brief Initiate an empty sparse vector
A Sparse_Mat consists of colums that have the type Sparse_Vec. The maximum number of non-zero elements is each column
is denoted \c row_data_init.
\param rows Number of rows in the matrix
\param cols Number of columns in the matrix
\param row_data_init The maximum number of non-zero elements in each column (default value is 200)
*/
Sparse_Mat(int rows, int cols, int row_data_init=200);
//! Initiate a new sparse matrix. The elements of \c m are copied into the new sparse matrix
Sparse_Mat(const Sparse_Mat<T> &m);
//! Initiate a new sparse matrix from a dense matrix. The elements of \c m are copied into the new sparse matrix
Sparse_Mat(const Mat<T> &m);
/*!
\brief Initiate a new sparse matrix from a dense matrix. Elements of \c m larger than \c epsilon are copied into the new sparse matrix.
\note If the type T is double complex, then the elements of \c m larger than \c abs(epsilon) are copied into the new sparse matrix.
*/
Sparse_Mat(const Mat<T> &m, T epsilon);
//! Destructor
~Sparse_Mat();
/*!
\brief Set the size of the sparse matrix
A Sparse_Mat consists of colums that have the type Sparse_Vec. The maximum number of non-zero elements is each column
is denoted \c row_data_init, with default value =-1 indicating that the number of data elements is not changed.
\param rows Number of rows in the matrix
\param cols Number of columns in the matrix
\param row_data_init The maximum number of non-zero elements in each column (default value -1 \c => allocated size for the data is not changed)
*/
void set_size(int rows, int cols, int row_data_init=-1);
//! Returns the number of rows of the sparse matrix
int rows() const { return n_rows; }
//! Returns the number of columns of the sparse matrix
int cols() const { return n_cols; }
//! The number of non-zero elements in the sparse matrix
int nnz();
//! Returns the density of the sparse matrix: (number of non-zero elements)/(total number of elements)
double density();
//! Set the maximum number of non-zero elements in each column equal to the actual number of non-zero elements in each column
void compact();
//! Returns a full, dense matrix in \c m
void full(Mat<T> &m) const;
//! Returns a full, dense matrix
Mat<T> full() const;
//! Returns element of row \c r and column \c c
T operator()(int r, int c) const;
//! Set element (\c r, \c c ) equal to \c v
void set(int r, int c, T v);
//! Set a new element with index (\c r, \c c ) equal to \c v
void set_new(int r, int c, T v);
//! Add the element in row \c r and column \c c with \c v
void add_elem(const int r, const int c, const T v);
//! Set the sparse matrix to the all zero matrix (removes all non-zero elements)
void zeros();
//! Set the element in row \c r and column \c c to zero (i.e. clear that element if it contains a non-zero value)
void zero_elem(const int r, const int c);
//! Clear all non-zero elements of the sparse matrix
void clear();
//! Clear the element in row \c r and column \c c (if it contains a non-zero value)
void clear_elem(const int r, const int c);
//! Set submatrix defined by rows r1,r2 and columns c1,c2 to matrix m
void set_submatrix(int r1, int r2, int c1, int c2, const Mat<T> &m);
//! Set submatrix defined by upper-left element (\c r,\c c) and the size of matrix \c m to \c m
void set_submatrix(int r, int c, const Mat<T>& m);
//! Returns the sub-matrix from rows \c r1 to \c r2 and columns \c c1 to \c c2
Sparse_Mat<T> get_submatrix(int r1, int r2, int c1, int c2) const;
//! Returns the sub-matrix from columns \c c1 to \c c2 (all rows)
Sparse_Mat<T> get_submatrix_cols(int c1, int c2) const;
//! Returns column \c c of the Sparse_Mat in the Sparse_Vec \c v
void get_col(int c, Sparse_Vec<T> &v) const;
//! Returns column \c c of the Sparse_Mat
Sparse_Vec<T> get_col(int c) const;
//! Set column \c c of the Sparse_Mat
void set_col(int c, const Sparse_Vec<T> &v);
//! Transpose the sparse matrix, return the result in \c m
void transpose(Sparse_Mat<T> &m) const;
//! Returns the transpose of the sparse matrix
Sparse_Mat<T> transpose() const;
//! Returns the transpose of the sparse matrix
// Sparse_Mat<T> T() const { return this->transpose(); };
//! Assign sparse matrix the value and dimensions of the sparse matrix \c m
void operator=(const Sparse_Mat<T> &m);
//! Assign sparse matrix the value and dimensions of the dense matrix \c m
void operator=(const Mat<T> &m);
//! Returns the sign inverse of all elements in the sparse matrix
Sparse_Mat<T> operator-() const;
//! Compare two sparse matricies. False if wrong sizes or different values
bool operator==(const Sparse_Mat<T> &m) const;
//! Add sparse matrix \c v to all non-zero elements of the sparse matrix
void operator+=(const Sparse_Mat<T> &v);
//! Add matrix \c v to all non-zero elements of the sparse matrix
void operator+=(const Mat<T> &v);
//! Subtract sparse matrix \c v from all non-zero elements of the sparse matrix
void operator-=(const Sparse_Mat<T> &v);
//! Subtract matrix \c v from all non-zero elements of the sparse matrix
void operator-=(const Mat<T> &v);
//! Multiply all non-zero elements of the sparse matrix with the scalar \c v
void operator*=(const T &v);
//! Divide all non-zero elements of the sparse matrix with the scalar \c v
void operator/=(const T &v);
//! Addition m1+m2 where m1 and m2 are sparse matrices
friend Sparse_Mat<T> operator+<>(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
//! Multiplication m1*m2 where m1 and m2 are sparse matrices
friend Sparse_Mat<T> operator*<>(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
//! Multiplication m*v where m is a sparse matrix and v is a full column vector
friend Vec<T> operator*<>(const Sparse_Mat<T> &m, const Vec<T> &v);
//! Multiplication v'*m where m is a sparse matrix and v is a full column vector
friend Vec<T> operator*<>(const Vec<T> &v, const Sparse_Mat<T> &m);
//! Multiplication m'*m where m is a sparse matrix. Returns a full, dense matrix
friend Mat<T> trans_mult <>(const Sparse_Mat<T> &m);
//! Multiplication m'*m where m is a sparse matrix, Returns a sparse matrix
friend Sparse_Mat<T> trans_mult_s <>(const Sparse_Mat<T> &m);
//! Multiplication m1'*m2 where m1 and m2 are sparse matrices
friend Sparse_Mat<T> trans_mult <>(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
//! Multiplication m'*v where m is a sparse matrix and v is a full column vector
friend Vec<T> trans_mult <>(const Sparse_Mat<T> &m, const Vec<T> &v);
//! Multiplication m1*m2' where m1 and m2 are sparse matrices
friend Sparse_Mat<T> mult_trans <>(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2);
private:
void init();
void alloc_empty();
void alloc(int row_data_size=200);
void free();
int n_rows, n_cols;
Sparse_Vec<T> *col;
};
/*!
\relates Sparse_Mat
\brief Sparse integer matrix
*/
typedef Sparse_Mat<int> sparse_imat;
/*!
\relates Sparse_Mat
\brief Sparse double matrix
*/
typedef Sparse_Mat<double> sparse_mat;
/*!
\relates Sparse_Mat
\brief Sparse complex<double> matrix
*/
typedef Sparse_Mat<complex<double> > sparse_cmat;
//---------------------- Implementation starts here --------------------------------
template <class T>
void Sparse_Mat<T>::init()
{
n_rows = 0;
n_cols = 0;
col = 0;
}
template <class T>
void Sparse_Mat<T>::alloc_empty()
{
if (n_cols == 0)
col = 0;
else
col = new Sparse_Vec<T>[n_cols];
}
template <class T>
void Sparse_Mat<T>::alloc(int row_data_init)
{
if (n_cols == 0)
col = 0;
else
col = new Sparse_Vec<T>[n_cols];
for (int c=0; c<n_cols; c++)
col[c].set_size(n_rows, row_data_init);
}
template <class T>
void Sparse_Mat<T>::free()
{
delete [] col;
col = 0;
}
template <class T>
Sparse_Mat<T>::Sparse_Mat()
{
init();
}
template <class T>
Sparse_Mat<T>::Sparse_Mat(int rows, int cols, int row_data_init)
{
init();
n_rows = rows;
n_cols = cols;
alloc(row_data_init);
}
template <class T>
Sparse_Mat<T>::Sparse_Mat(const Sparse_Mat<T> &m)
{
init();
n_rows = m.n_rows;
n_cols = m.n_cols;
alloc_empty();
for (int c=0; c<n_cols; c++)
col[c] = m.col[c];
}
template <class T>
Sparse_Mat<T>::Sparse_Mat(const Mat<T> &m)
{
init();
n_rows = m.rows();
n_cols = m.cols();
alloc();
for (int c=0; c<n_cols; c++) {
for (int r=0; r<n_rows; r++) {
//if (abs(m(r,c)) != T(0))
if (m(r,c) != T(0))
col[c].set_new(r, m(r,c));
}
col[c].compact();
}
}
template <class T>
Sparse_Mat<T>::Sparse_Mat(const Mat<T> &m, T epsilon)
{
init();
n_rows = m.rows();
n_cols = m.cols();
alloc();
for (int c=0; c<n_cols; c++) {
for (int r=0; r<n_rows; r++) {
if (std::abs(m(r,c)) > std::abs(epsilon))
col[c].set_new(r, m(r,c));
}
col[c].compact();
}
}
template <class T>
Sparse_Mat<T>::~Sparse_Mat()
{
free();
}
template <class T>
void Sparse_Mat<T>::set_size(int rows, int cols, int row_data_init)
{
n_rows = rows;
//Allocate new memory for data if the number of columns has changed or if row_data_init != -1
if (cols!=n_cols||row_data_init!=-1) {
n_cols = cols;
free();
alloc(row_data_init);
}
}
template <class T>
int Sparse_Mat<T>::nnz()
{
int n=0;
for (int c=0; c<n_cols; c++)
n += col[c].nnz();
return n;
}
template <class T>
double Sparse_Mat<T>::density()
{
//return static_cast<double>(nnz())/(n_rows*n_cols);
return double(nnz())/(n_rows*n_cols);
}
template <class T>
void Sparse_Mat<T>::compact()
{
for (int c=0; c<n_cols; c++)
col[c].compact();
}
template <class T>
void Sparse_Mat<T>::full(Mat<T> &m) const
{
m.set_size(n_rows, n_cols);
m = T(0);
for (int c=0; c<n_cols; c++) {
for (int p=0; p<col[c].nnz(); p++)
m(col[c].get_nz_index(p),c) = col[c].get_nz_data(p);
}
}
template <class T>
Mat<T> Sparse_Mat<T>::full() const
{
Mat<T> r(n_rows, n_cols);
full(r);
return r;
}
template <class T>
T Sparse_Mat<T>::operator()(int r, int c) const
{
it_assert0(r>=0&&r<n_rows&&c>=0&&c<n_cols, "Incorrect input indexes given");
return col[c](r);
}
template <class T>
void Sparse_Mat<T>::set(int r, int c, T v)
{
it_assert0(r>=0&&r<n_rows&&c>=0&&c<n_cols, "Incorrect input indexes given");
col[c].set(r, v);
}
template <class T>
void Sparse_Mat<T>::set_new(int r, int c, T v)
{
it_assert0(r>=0&&r<n_rows&&c>=0&&c<n_cols, "Incorrect input indexes given");
col[c].set_new(r, v);
}
template <class T>
void Sparse_Mat<T>::add_elem(int r, int c, T v)
{
it_assert0(r>=0&&r<n_rows&&c>=0&&c<n_cols, "Incorrect input indexes given");
col[c].add_elem(r, v);
}
template <class T>
void Sparse_Mat<T>::zeros()
{
for (int c=0; c<n_cols; c++)
col[c].zeros();
}
template <class T>
void Sparse_Mat<T>::zero_elem(const int r, const int c)
{
it_assert0(r>=0&&r<n_rows&&c>=0&&c<n_cols, "Incorrect input indexes given");
col[c].zero_elem(r);
}
template <class T>
void Sparse_Mat<T>::clear()
{
for (int c=0; c<n_cols; c++)
col[c].clear();
}
template <class T>
void Sparse_Mat<T>::clear_elem(const int r, const int c)
{
it_assert0(r>=0&&r<n_rows&&c>=0&&c<n_cols, "Incorrect input indexes given");
col[c].clear_elem(r);
}
template <class T>
void Sparse_Mat<T>::set_submatrix(int r1, int r2, int c1, int c2, const Mat<T>& m)
{
if (r1 == -1) r1 = n_rows-1;
if (r2 == -1) r2 = n_rows-1;
if (c1 == -1) c1 = n_cols-1;
if (c2 == -1) c2 = n_cols-1;
it_assert1(r1>=0 && r2>=0 && r1<n_rows && r2<n_rows &&
c1>=0 && c2>=0 && c1<n_cols && c2<n_cols, "Sparse_Mat<Num_T>::set_submatrix(): index out of range");
it_assert1(r2>=r1 && c2>=c1, "Sparse_Mat<Num_T>::set_submatrix: r2<r1 or c2<c1");
it_assert1(m.rows() == r2-r1+1 && m.cols() == c2-c1+1, "Mat<Num_T>::set_submatrix(): sizes don't match");
for (int i=0 ; i<m.rows() ; i++) {
for (int j=0 ; j<m.cols() ; j++) {
set(r1+i, c1+j, m(i,j));
}
}
}
template <class T>
void Sparse_Mat<T>::set_submatrix(int r, int c, const Mat<T>& m)
{
it_assert1(r>=0 && r+m.rows()<=n_rows &&
c>=0 && c+m.cols()<=n_cols, "Sparse_Mat<Num_T>::set_submatrix(): index out of range");
for (int i=0 ; i<m.rows() ; i++) {
for (int j=0 ; j<m.cols() ; j++) {
set(r+i, c+j, m(i,j));
}
}
}
template <class T>
Sparse_Mat<T> Sparse_Mat<T>::get_submatrix(int r1, int r2, int c1, int c2) const
{
it_assert0(r1<=r2 && r1>=0 && r1<n_rows && c1<=c2 && c1>=0 && c1<n_cols,
"Sparse_Mat<T>::get_submatrix(): illegal input variables");
Sparse_Mat<T> r(r2-r1+1, c2-c1+1);
for (int c=c1; c<=c2; c++)
r.col[c-c1] = col[c].get_subvector(r1, r2);
r.compact();
return r;
}
template <class T>
Sparse_Mat<T> Sparse_Mat<T>::get_submatrix_cols(int c1, int c2) const
{
it_assert0(c1<=c2 && c1>=0 && c1<n_cols, "Sparse_Mat<T>::get_submatrix_cols()");
Sparse_Mat<T> r(n_rows, c2-c1+1, 0);
for (int c=c1; c<=c2; c++)
r.col[c-c1] = col[c];
r.compact();
return r;
}
template <class T>
void Sparse_Mat<T>::get_col(int c, Sparse_Vec<T> &v) const
{
it_assert(c>=0 && c<n_cols, "Sparse_Mat<T>::get_col()");
v = col[c];
}
template <class T>
Sparse_Vec<T> Sparse_Mat<T>::get_col(int c) const
{
it_assert(c>=0 && c<n_cols, "Sparse_Mat<T>::get_col()");
return col[c];
}
template <class T>
void Sparse_Mat<T>::set_col(int c, const Sparse_Vec<T> &v)
{
it_assert(c>=0 && c<n_cols, "Sparse_Mat<T>::set_col()");
col[c] = v;
}
template <class T>
void Sparse_Mat<T>::transpose(Sparse_Mat<T> &m) const
{
m.set_size(n_cols, n_rows);
for (int c=0; c<n_cols; c++) {
for (int p=0; p<col[c].nnz(); p++)
m.col[col[c].get_nz_index(p)].set_new(c, col[c].get_nz_data(p));
}
}
template <class T>
Sparse_Mat<T> Sparse_Mat<T>::transpose() const
{
Sparse_Mat<T> m;
transpose(m);
return m;
}
template <class T>
void Sparse_Mat<T>::operator=(const Sparse_Mat<T> &m)
{
free();
n_rows = m.n_rows;
n_cols = m.n_cols;
alloc_empty();
for (int c=0; c<n_cols; c++)
col[c] = m.col[c];
}
template <class T>
void Sparse_Mat<T>::operator=(const Mat<T> &m)
{
free();
n_rows = m.rows();
n_cols = m.cols();
alloc();
for (int c=0; c<n_cols; c++) {
for (int r=0; r<n_rows; r++) {
if (m(r,c) != T(0))
col[c].set_new(r, m(r,c));
}
col[c].compact();
}
}
template <class T>
Sparse_Mat<T> Sparse_Mat<T>::operator-() const
{
Sparse_Mat r(n_rows, n_cols, 0);
for (int c=0; c<n_cols; c++) {
r.col[c].resize_data(col[c].nnz());
for (int p=0; p<col[c].nnz(); p++)
r.col[c].set_new(col[c].get_nz_index(p), -col[c].get_nz_data(p));
}
return r;
}
template <class T>
bool Sparse_Mat<T>::operator==(const Sparse_Mat<T> &m) const
{
if (n_rows!=m.n_rows || n_cols!=m.n_cols)
return false;
for (int c=0; c<n_cols; c++) {
if (!(col[c] == m.col[c]))
return false;
}
// If they passed all tests, they must be equal
return true;
}
template <class T>
void Sparse_Mat<T>::operator+=(const Sparse_Mat<T> &m)
{
it_assert0(m.rows()==n_rows&&m.cols()==n_cols, "Addition of unequal sized matrices is not allowed");
Sparse_Vec<T> v;
for (int c=0; c<n_cols; c++) {
m.get_col(c,v);
col[c]+=v;
}
}
template <class T>
void Sparse_Mat<T>::operator+=(const Mat<T> &m)
{
it_assert0(m.rows()==n_rows&&m.cols()==n_cols, "Addition of unequal sized matrices is not allowed");
for (int c=0; c<n_cols; c++)
col[c]+=(m.get_col(c));
}
template <class T>
void Sparse_Mat<T>::operator-=(const Sparse_Mat<T> &m)
{
it_assert0(m.rows()==n_rows&&m.cols()==n_cols, "Subtraction of unequal sized matrices is not allowed");
Sparse_Vec<T> v;
for (int c=0; c<n_cols; c++) {
m.get_col(c,v);
col[c]-=v;
}
}
template <class T>
void Sparse_Mat<T>::operator-=(const Mat<T> &m)
{
it_assert0(m.rows()==n_rows&&m.cols()==n_cols, "Subtraction of unequal sized matrices is not allowed");
for (int c=0; c<n_cols; c++)
col[c]-=(m.get_col(c));
}
template <class T>
void Sparse_Mat<T>::operator*=(const T &m)
{
for (int c=0; c<n_cols; c++)
col[c]*=m;
}
template <class T>
void Sparse_Mat<T>::operator/=(const T &m)
{
for (int c=0; c<n_cols; c++)
col[c]/=m;
}
template <class T>
Sparse_Mat<T> operator+(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2)
{
it_assert0(m1.n_cols == m2.n_cols && m1.n_rows == m2.n_rows , "Sparse_Mat<T> + Sparse_Mat<T>");
Sparse_Mat<T> m(m1.n_rows, m1.n_cols, 0);
for (int c=0; c<m.n_cols; c++)
m.col[c] = m1.col[c] + m2.col[c];
return m;
}
template <class T>
Sparse_Mat<T> operator*(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2)
{
it_assert0(m1.n_cols == m2.n_rows, "Sparse_Mat<T> * Sparse_Mat<T>");
Sparse_Mat<T> ret(m1.n_rows, m2.n_cols);
ivec occupied_by(ret.n_rows), pos(ret.n_rows);
for (int rp=0; rp<m1.n_rows; rp++)
occupied_by[rp] = -1;
for (int c=0; c<ret.n_cols; c++) {
Sparse_Vec<T> &m2col = m2.col[c];
for (int p2=0; p2<m2col.nnz(); p2++) {
Sparse_Vec<T> &m1col = m1.col[m2col.get_nz_index(p2)];
for (int p1=0; p1<m1col.nnz(); p1++) {
int r = m1col.get_nz_index(p1);
T inc = m1col.get_nz_data(p1) * m2col.get_nz_data(p2);
if (occupied_by[r] == c) {
int index=ret.col[c].get_nz_index(pos[r]);
ret.col[c].add_elem(index,inc);
}
else {
occupied_by[r] = c;
pos[r] = ret.col[c].nnz();
ret.col[c].set_new(r, inc);
}
}
}
}
ret.compact();
return ret;
}
template <class T>
Vec<T> operator*(const Sparse_Mat<T> &m, const Vec<T> &v)
{
it_assert0(m.n_cols == v.size(), "Sparse_Mat<T> * Vec<T>");
Vec<T> r(m.n_rows);
r.clear();
//r = T(0);
for (int c=0; c<m.n_cols; c++) {
for (int p=0; p<m.col[c].nnz(); p++)
r(c) += m.col[c].get_nz_data(p) * v(c);
}
return r;
}
template <class T>
Vec<T> operator*(const Vec<T> &v, const Sparse_Mat<T> &m)
{
it_assert0(v.size() == m.n_rows, "Vec<T> * Sparse_Mat<T>");
Vec<T> r(m.n_cols);
r.clear();
//r = T(0);
for (int c=0; c<m.n_cols; c++)
r[c] = v * m.col[c];
return r;
}
template <class T>
Mat<T> trans_mult(const Sparse_Mat<T> &m)
{
Mat<T> ret(m.n_cols, m.n_cols);
Vec<T> col;
for (int c=0; c<ret.cols(); c++) {
m.col[c].full(col);
for (int r=0; r<c; r++) {
T tmp = m.col[r] * col;
ret(r,c) = tmp;
ret(c,r) = tmp;
}
ret(c,c) = m.col[c].sqr();
}
return ret;
}
template <class T>
Sparse_Mat<T> trans_mult_s(const Sparse_Mat<T> &m)
{
Sparse_Mat<T> ret(m.n_cols, m.n_cols);
Vec<T> col;
T tmp;
for (int c=0; c<ret.n_cols; c++) {
m.col[c].full(col);
for (int r=0; r<c; r++) {
tmp = m.col[r] * col;
if (tmp != T(0)) {
ret.col[c].set_new(r, tmp);
ret.col[r].set_new(c, tmp);
}
}
tmp = m.col[c].sqr();
if (tmp != T(0))
ret.col[c].set_new(c, tmp);
}
return ret;
}
template <class T>
Sparse_Mat<T> trans_mult(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2)
{
it_assert0(m1.n_rows == m2.n_rows, "trans_mult()");
Sparse_Mat<T> ret(m1.n_cols, m2.n_cols);
Vec<T> col;
for (int c=0; c<ret.n_cols; c++) {
m2.col[c].full(col);
for (int r=0; r<ret.n_rows; r++)
ret.col[c].set_new(r, m1.col[r] * col);
}
return ret;
}
template <class T>
Vec<T> trans_mult(const Sparse_Mat<T> &m, const Vec<T> &v)
{
Vec<T> r(m.n_cols);
for (int c=0; c<m.n_cols; c++)
r(c) = m.col[c] * v;
return r;
}
template <class T>
Sparse_Mat<T> mult_trans(const Sparse_Mat<T> &m1, const Sparse_Mat<T> &m2)
{
return trans_mult(m1.transpose(),m2.transpose());
}
template <class T>
inline Sparse_Mat<T> sparse(const Mat<T> &m, T epsilon)
{
Sparse_Mat<T> s(m, epsilon);
return s;
}
template <class T>
inline Mat<T> full(const Sparse_Mat<T> &s)
{
Mat<T> m;
s.full(m);
return m;
}
template <class T>
inline Sparse_Mat<T> transpose(const Sparse_Mat<T> &s)
{
Sparse_Mat<T> m;
s.transpose(m);
return m;
}
} //namespace itpp
#endif //__smat_h
syntax highlighted by Code2HTML, v. 0.9.1