/*---------------------------------------------------------------------------* * 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 Matrix Class Definitions \author Tony Ottosson and Tobias Ringstrom 1.41 2004/09/05 03:34:22 */ #ifndef __mat_h #define __mat_h #include #include #include #include #include #include #include #include "itconfig.h" //#include "base/vec.h" is done below Mat class definition (needed for g++ 3.4) #include "base/itassert.h" #include "base/binary.h" #include "base/scalfunc.h" #include "base/factory.h" //#ifndef NO_CBLAS //#include "base/cblas.h" //#endif using std::cout; using std::endl; using std::string; using std::ostream; using std::istream; using std::istringstream; using std::getline; using std::swap; using std::complex; namespace itpp { //! Define of minimum of two numbers x and y #define minimum(x,y) (((x)<(y)) ? (x) : (y)) // Declaration of Vec template class Vec; // Declaration of Mat template class Mat; // Declaration of bin class bin; // ------------------------------------------------------------------------------------- // Declaration of Mat Friends // ------------------------------------------------------------------------------------- //! Horizontal concatenation of two matrices template Mat concat_horizontal(const Mat &m1, const Mat &m2); //! Vertical concatenation of two matrices template Mat concat_vertical(const Mat &m1, const Mat &m2); //! Addition of two matricies template Mat operator+(const Mat &m1, const Mat &m2); //! Addition of a matrix and a scalar template Mat operator+(const Mat &m, Num_T t); //! Addition of a scalar and a matrix template Mat operator+(Num_T t, const Mat &m); //! Subtraction of two matrices template Mat operator-(const Mat &m1, const Mat &m2); //! Subtraction of matrix and scalar template Mat operator-(const Mat &m, Num_T t); //! Subtraction of scalar and matrix template Mat operator-(Num_T t, const Mat &m); //! Negation of matrix template Mat operator-(const Mat &m); //! Multiplication of two matricies template Mat operator*(const Mat &m1, const Mat &m2); //! Multiplication of matrix and vector template Vec operator*(const Mat &m, const Vec &v); //! Multiplication of vector and matrix (only works for a matrix that is a row vector) template Mat operator*(const Vec &v, const Mat &m); //! Multiplication of matrix and scalar template Mat operator*(const Mat &m, Num_T t); //! Multiplication of scalar and matrix template Mat operator*(Num_T t, const Mat &m); //! Element wise multiplication of two matricies template Mat elem_mult(const Mat &m1, const Mat &m2); //! Division of matrix and scalar template Mat operator/(const Mat &m, Num_T t); //! Element wise division of two matricies template Mat elem_div(const Mat &m1, const Mat &m2); // ------------------------------------------------------------------------------------- // Declaration of Mat // ------------------------------------------------------------------------------------- /*! \brief Templated Matrix Class \author Tony Ottosson and Tobias Ringstrom Matrices can be of arbitrarily types, but conversions and functions are prepared for \c bin, \c short, \c int, \c double, and \c complex vectors and these are predefined as: \c bmat, \c smat, \c imat, \c mat, and \c cmat. \c double and \c complex are usually \c double and \c complex respectively. However, this can be changed when \ compiling the it++ (see installation notes for more details). Examples: Matrix Constructors: When constructing a matrix without a dimensions (memory) use \code mat temp; \endcode For construction of a matrix of a given size use \code mat temp(rows, cols); \endcode It is also possible to assign the constructed matrix the value and dimension of another matrix by \code vec temp(inmatrix); \endcode If you have explicit values you would like to assign to the matrix it is possible to do this using strings as: \code mat a("0 0.7;5 9.3"); // that is a = [0, 0.7; 5, 9.3] mat a="0 0.7;5 9.3"; // the constructor are called implicitly \endcode It is also possible to change dimension by \code temp.set_size(new_rows, new_cols, false); \endcode where \c false is used to indicate that the old values in \c temp is not copied. If you like to preserve the values use \c true. There are a number of methods to access parts of a matrix. Examples are \code a(5,3); // Element number (5,3) a(5,9,3,5); // Sub-matrix from rows 5, 6, 7, 8, 9 the columns 3, 4, and 5 a.get_row(10); // Row 10 a.get_col(10); // Column 10 \endcode It is also possible to modify parts of a vector as e.g. in \code a.set_row(5, invector); // Set row 5 to \c invector a.set_col(3, invector); // Set column 3 to \c invector a.copy_col(1, 5); // Copy column 5 to column 1 a.swap_cols(1, 5); // Swap the contents of columns 1 and 5 \endcode It is of course also possible to perform the common linear algebra methods such as addition, subtraction, and matrix multiplication. Observe though, that vectors are assumed to be column-vectors in operations with matrices. Most elementary functions such as sin(), cosh(), log(), abs(), ..., are available as operations on the individual elements of the matrices. Please see the individual functions for more details. By default, the Mat elements are created using the default constructor for the element type. This can be changed by specifying a suitable Factory in the Mat constructor call; see Detailed Description for Factory. */ template class Mat { public: //! Default constructor. An element factory \c f can be specified explicit Mat(const Factory &f = DEFAULT_FACTORY) : factory(f) { init(); } //! Create a matrix of size (inrow, incol). An element factory \c f can be specified Mat(int inrow, int incol, const Factory &f = DEFAULT_FACTORY) : factory(f) { init(); it_assert1(inrow>=0 && incol>=0, "The rows and columns must be >= 0"); alloc(inrow, incol); } //! Copy constructor Mat(const Mat &m); //! Constructor, similar to the copy constructor, but also takes an element factory \c f as argument Mat(const Mat &m, const Factory &f); //! Create a copy of the vector \c invector treated as a column vector. An element factory \c f can be specified Mat(const Vec &invector, const Factory &f = DEFAULT_FACTORY); //! Set matrix equal to values in string. An element factory \c f can be specified Mat(const std::string &str, const Factory &f = DEFAULT_FACTORY) : factory(f) { init(); set(str); } //! Set matrix equal to values in string. An element factory \c f can be specified Mat(const char *str, const Factory &f = DEFAULT_FACTORY) : factory(f) { init(); set(str); } /*! \brief Constructor taking a C-array as input. Copies all data. An element factory \c f can be specified By default the matrix is stored as a RowMajor matrix (i.e. listing elements in sequence beginning with the first column). */ Mat(Num_T *c_array, int rows, int cols, bool RowMajor, const Factory &f = DEFAULT_FACTORY); //! Destructor ~Mat() { free(); } //! The number of columns int cols() const { return no_cols; } //! The number of rows int rows() const { return no_rows; } //! The number of elements int size() const { return datasize; } //! Set size of matrix. If copy = true then keep the data before resizing. void set_size(int inrow, int incol, bool copy=false); //! Set matrix equal to the all zero matrix void zeros(); //! Set matrix equal to the all zero matrix void clear() { zeros(); } //! Set matrix equal to the all one matrix void ones(); //! Set matrix equal to values in \c values bool set(const char *values); //! Set matrix equal to values in the string \c str bool set(const std::string &str); //! Get element (R,C) from matrix const Num_T &operator()(int R,int C) const { it_assert0(R>=0 && R=0 && C::operator(): index out of range"); return data[R+C*no_rows]; } //! Get element (R,C) from matrix Num_T &operator()(int R,int C) { it_assert0(R>=0 && R=0 && C::operator(): index out of range"); return data[R+C*no_rows]; } //! Get element \c index using linear addressing (by rows) Num_T &operator()(int index) { it_assert0(index=0,"Mat::operator(): index out of range"); return data[index]; } //! Get element \c index using linear addressing (by rows) const Num_T &operator()(int index) const { it_assert0(index=0,"Mat::operator(): index out of range"); return data[index]; } //! Get element (R,C) from matrix const Num_T &get(int R,int C) const { it_assert0(R>=0 && R=0 && C::get(): index out of range"); return data[R+C*no_rows]; } //! Set element (R,C) of matrix void set(int R,int C, const Num_T &v) { it_assert0(R>=0 && R=0 && C::set(): index out of range"); data[R+C*no_rows] = v; } /*! \brief Sub-matrix from row \c r1 to row \c r2 and columns \c c1 to \c c2. Value -1 indicates the last row and column, respectively. */ const Mat operator()(int r1, int r2, int c1, int c2) const; /*! \brief Sub-matrix from row \c r1 to row \c r2 and columns \c c1 to \c c2. Value -1 indicates the last row and column, respectively. */ const Mat get(int r1, int r2, int c1, int c2) const; //! Get row \c Index Vec get_row(int Index) const ; //! Get rows \c r1 through \c r2 Mat get_rows(int r1, int r2) const; //! Get the rows specified by \c indexlist Mat get_rows(const Vec &indexlist) const; //! Get column \c Index Vec get_col(int Index) const ; //! Get columns \c c1 through \c c2 Mat get_cols(int c1, int c2) const; //! Get the columns specified by \c indexlist Mat get_cols(const Vec &indexlist) const; //! Set row \c Index to \c invector void set_row(int Index, const Vec &invector); //! Set column \c Index to \c invector void set_col(int Index, const Vec &invector); //! Copy row \c from onto row \c to void copy_row(int to, int from); //! Copy column \c from onto column \c to void copy_col(int to, int from); //! Swap the rows \c r1 and \c r2 void swap_rows(int r1, int r2); //! Swap the columns \c c1 and \c c2 void swap_cols(int c1, int c2); //! 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 &m); //! Set submatrix defined by upper-left element (r,c) and the size of matrix m to m void set_submatrix(int r, int c, const Mat &m); //! Set all elements of submatrix defined by rows r1,r2 and columns c1,c2 to value t void set_submatrix(int r1, int r2, int c1, int c2, const Num_T t); //! Delete row number \c r void del_row(int r); //! Delete rows from \c r1 to \c r2 void del_rows(int r1, int r2); //! Delete column number \c c void del_col(int c); //! Delete columns from \c c1 to \c c2 void del_cols(int c1, int c2); //! Insert vector \c in at row number \c r, the matrix can be empty. void ins_row(int r, const Vec &in); //! Insert vector \c in at column number \c c, the matrix can be empty. void ins_col(int c, const Vec &in); // ! Append vector \c to the bottom of the matrix, the matrix can be empty. void append_row(const Vec &in); // ! Append vector \c to the right of the matrix, the matrix can be empty. void append_col(const Vec &in); //! Matrix transpose Mat transpose() const; //! Matrix transpose Mat T() const { return this->transpose(); } //! Hermitian matrix transpose (conjugate transpose) Mat hermitian_transpose() const; //! Hermitian matrix transpose (conjugate transpose) Mat H() const { return this->hermitian_transpose(); } //! Concatenate the matrices \c m1 and \c m2 horizontally friend Mat concat_horizontal <>(const Mat &m1, const Mat &m2); //! Concatenate the matrices \c m1 and \c m2 vertically friend Mat concat_vertical <>(const Mat &m1, const Mat &m2); //! Set all elements of the matrix equal to \c t void operator=(Num_T t); //! Set matrix equal to \c m void operator=(const Mat &m); //! Set matrix equal to the vector \c v, assuming column vector void operator=(const Vec &v); //! Set matrix equal to values in the string void operator=(const char *values) { set(values); } //! Addition of matrices void operator+=(const Mat &m); //! Addition of scalar to matrix void operator+=(Num_T t); //! Addition of two matrices friend Mat operator+<>(const Mat &m1, const Mat &m2); //! Addition of matrix and scalar friend Mat operator+<>(const Mat &m, Num_T t); //! Addition of scalar and matrix friend Mat operator+<>(Num_T t, const Mat &m); //! Subtraction of matrix void operator-=(const Mat &m); //! Subtraction of scalar from matrix void operator-=(Num_T t); //! Subtraction of \c m2 from \c m1 friend Mat operator-<>(const Mat &m1, const Mat &m2); //! Subraction of scalar from matrix friend Mat operator-<>(const Mat &m, Num_T t); //! Subtract matrix from scalar friend Mat operator-<>(Num_T t, const Mat &m); //! Subraction of matrix friend Mat operator-<>(const Mat &m); //! Matrix multiplication void operator*=(const Mat &m); //! Multiplication by a scalar void operator*=(Num_T t); //! Multiplication of two matrices friend Mat operator*<>(const Mat &m1, const Mat &m2); //! Multiplication of matrix \c m and vector \c v (column vector) friend Vec operator*<>(const Mat &m, const Vec &v); //! Multiplication of transposed vector \c v and matrix \c m friend Mat operator*<>(const Vec &v, const Mat &m); //! Multiplication of matrix and scalar friend Mat operator*<>(const Mat &m, Num_T t); //! Multiplication of scalar and matrix friend Mat operator*<>(Num_T t, const Mat &m); //! Elementwise multiplication of two matrices friend Mat elem_mult <>(const Mat &m1, const Mat &m2); //! Division by a scalar void operator/=(Num_T t); //! Division of matrix with scalar friend Mat operator/<>(const Mat &m, Num_T t); //! Elementwise division with the current matrix void operator/=(const Mat &m); //! Elementwise division of matrix \c m1 with matrix \c m2 friend Mat elem_div <>(const Mat &m1, const Mat &m2); //! Compare two matrices. False if wrong sizes or different values bool operator==(const Mat &m) const; //! Compare two matrices. True if different bool operator!=(const Mat &m) const; //! Get element (R,C) from matrix without boundary check (Not recommended to use). Num_T &_elem(int R,int C) { return data[R+C*no_rows]; } //! Get element (R,C) from matrix without boundary check (Not recommended to use). const Num_T &_elem(int R,int C) const { return data[R+C*no_rows]; } //! Get element \c index using linear addressing (by rows) without boundary check (Not recommended to use). Num_T &_elem(int index) { return data[index]; } //! Get element \c index using linear addressing (by rows) without boundary check (Not recommended to use). const Num_T &_elem(int index) const { return data[index]; } //! Access of the internal data structure. Don't use. May be changed! Num_T *_data() { return data; } //! Access to the internal data structure. Don't use. May be changed! const Num_T *_data() const { return data; } //! Access to the internal data structure. Don't use. May be changed! int _datasize() const { return datasize; } protected: //! Allocate memory for the matrix void alloc(int rows, int cols) { if ( datasize == rows * cols ) { // Reuse the memory no_rows = rows; no_cols = cols; return; } free(); // Free memory (if any allocated) if (rows == 0 || cols == 0) return; datasize = rows * cols; create_elements(data, datasize, factory); no_rows = rows; no_cols = cols; it_assert1(data, "Mat::alloc: Out of memory"); } //! Free the memory space of the matrix void free() { delete [] data; init(); } //! Protected integer variables int datasize, no_rows, no_cols; //! Protected data pointer Num_T *data; //! Element factory (set to DEFAULT_FACTORY to use Num_T default constructors only) const Factory &factory; private: void init() { data = 0; datasize = no_rows = no_cols = 0; } }; } //namespace itpp #include "base/vec.h" namespace itpp { // ------------------------------------------------------------------------------------- // Declaration of input and output streams for Mat // ------------------------------------------------------------------------------------- /*! \relates Mat \brief Output stream for matrices */ template std::ostream &operator<<(std::ostream &os, const Mat &m); /*! \relates Mat \brief Input stream for matrices The input can be on the form "1 2 3; 4 5 6" or "[[1 2 3][4 5 6]]", i.e. with brackets or semicolons as row delimiters. The first form is compatible with the set method, while the second form is compatible with the ostream operator. The elements on a row can be separated by blank space or commas. Rows that are shorter than the longest row are padded with zero elements. "[]" means an empty matrix. */ template std::istream &operator>>(std::istream &is, Mat &m); // ------------------------------------------------------------------------------------- // Type definitions of mat, cmat, imat, smat, and bmat // ------------------------------------------------------------------------------------- /*! \relates Mat \brief Default Matrix Type */ typedef Mat mat; /*! \relates Mat \brief Default Complex Matrix Type */ typedef Mat > cmat; /*! \relates Mat \brief Integer matrix */ typedef Mat imat; /*! \relates Mat \brief short int matrix */ typedef Mat smat; /*! \relates Mat \brief bin matrix */ typedef Mat bmat; // ------------------------------------------------------------------------------------- // Implementation of templated Mat members and friends // ------------------------------------------------------------------------------------- template inline Mat::Mat(const Mat &m) : factory(m.factory) { init(); alloc(m.no_rows, m.no_cols); copy_vector(m.datasize, m.data, data); } template inline Mat::Mat(const Mat &m, const Factory &f) : factory(f) { init(); alloc(m.no_rows, m.no_cols); copy_vector(m.datasize, m.data, data); } template inline Mat::Mat(const Vec &invector, const Factory &f) : factory(f) { init(); set_size(invector.length(), 1, false); set_col(0,invector); } template inline Mat::Mat(Num_T *c_array, int rows, int cols, bool RowMajor, const Factory &f) : factory(f) { init(); alloc(rows, cols); if (!RowMajor) copy_vector(datasize, c_array, data); else { // Row Major Num_T *ptr = c_array; for (int i=0; i inline void Mat::set_size(int inrow, int incol, bool copy) { it_assert1(inrow>=0 && incol>=0, "Mat::set_size: The rows and columns must be >= 0"); if (no_rows!=inrow || no_cols!=incol) { if (copy) { Mat temp(*this); int i, j; alloc(inrow, incol); for (i=0;i inline void Mat::zeros() { for(int i=0; i inline void Mat::ones() { for(int i=0; i bool cmat::set(const char *values); template<> bool bmat::set(const char *values); template bool Mat::set(const char *values) { istringstream buffer(values); int rows=0, maxrows=10, cols=0, nocols=0, maxcols=10; alloc(maxrows, maxcols); while (buffer.peek()!=EOF) { rows++; if (rows > maxrows) { maxrows=maxrows*2; set_size(maxrows, maxcols, true); } cols=0; while ( (buffer.peek() != ';') && (buffer.peek() != EOF) ) { if (buffer.peek()==',') { buffer.get(); } else { cols++; if (cols > nocols) { nocols=cols; if (cols > maxcols) { maxcols=maxcols*2; set_size(maxrows, maxcols, true); } } buffer >> this->operator()(rows-1,cols-1); while (buffer.peek()==' ') { buffer.get(); } } } if (!buffer.eof()) buffer.get(); } set_size(rows, nocols, true); return true; } template bool Mat::set(const std::string &str) { return set(str.c_str()); } template inline const Mat Mat::operator()(int r1, int r2, int c1, int c2) const { if (r1 == -1) r1 = no_rows-1; if (r2 == -1) r2 = no_rows-1; if (c1 == -1) c1 = no_cols-1; if (c2 == -1) c2 = no_cols-1; it_assert1(r1>=0 && r2>=0 && r1=0 && c2>=0 && c1=r1 && c2>=c1, "Mat::op(): r2>=r1 or c2>=c1 not fulfilled"); Mat s(r2-r1+1, c2-c1+1); for (int i=0;i inline const Mat Mat::get(int r1, int r2, int c1, int c2) const { return (*this)(r1, r2, c1, c2); } template inline Vec Mat::get_row(int Index) const { it_assert1(Index>=0 && Index::get_row: index out of range"); Vec a(no_cols); copy_vector(no_cols, data+Index, no_rows, a._data(), 1); return a; } template Mat Mat::get_rows(int r1, int r2) const { it_assert1(r1>=0 && r2::get_rows: index out of range"); Mat m(r2-r1+1, no_cols); for (int i=0; i Mat Mat::get_rows(const Vec &indexlist) const { Mat m(indexlist.size(),no_cols); for (int i=0;i=0 && indexlist(i)::get_rows: index out of range"); copy_vector(no_cols, data+indexlist(i), no_rows, m.data+i, m.no_rows); } return m; } template inline Vec Mat::get_col(int Index) const { it_assert1(Index>=0 && Index::get_col: index out of range"); Vec a(no_rows); copy_vector(no_rows, data+Index*no_rows, a._data()); return a; } template inline Mat Mat::get_cols(int c1, int c2) const { it_assert1(c1>=0 && c2::get_cols: index out of range"); Mat m(no_rows, c2-c1+1); for (int i=0; i inline Mat Mat::get_cols(const Vec &indexlist) const { Mat m(no_rows,indexlist.size()); for (int i=0; i=0 && indexlist(i)::get_cols: index out of range"); copy_vector(no_rows, data+indexlist(i)*no_rows, m.data+i*m.no_rows); } return m; } template inline void Mat::set_row(int Index, const Vec &v) { it_assert1(Index>=0 && Index::set_row: index out of range"); it_assert1(v.length() == no_cols, "Mat::set_row: lengths doesn't match"); copy_vector(v.size(), v._data(), 1, data+Index, no_rows); } template inline void Mat::set_col(int Index, const Vec &v) { it_assert1(Index>=0 && Index::set_col: index out of range"); it_assert1(v.length() == no_rows, "Mat::set_col: lengths doesn't match"); copy_vector(v.size(), v._data(), data+Index*no_rows); } template inline void Mat::copy_row(int to, int from) { it_assert1(to>=0 && from>=0 && to::copy_row: index out of range"); if (from == to) return; copy_vector(no_cols, data+from, no_rows, data+to, no_rows); } template inline void Mat::copy_col(int to, int from) { it_assert1(to>=0 && from>=0 && to::copy_col: index out of range"); if (from == to) return; copy_vector(no_rows, data+from*no_rows, data+to*no_rows); } template inline void Mat::swap_rows(int r1, int r2) { it_assert1(r1>=0 && r2>=0 && r1::swap_rows: index out of range"); if (r1 == r2) return; swap_vector(no_cols, data+r1, no_rows, data+r2, no_rows); } template inline void Mat::swap_cols(int c1, int c2) { it_assert1(c1>=0 && c2>=0 && c1::swap_cols: index out of range"); if (c1 == c2) return; swap_vector(no_rows, data+c1*no_rows, data+c2*no_rows); } template inline void Mat::set_submatrix(int r1, int r2, int c1, int c2, const Mat &m) { if (r1 == -1) r1 = no_rows-1; if (r2 == -1) r2 = no_rows-1; if (c1 == -1) c1 = no_cols-1; if (c2 == -1) c2 = no_cols-1; it_assert1(r1>=0 && r2>=0 && r1=0 && c2>=0 && c1::set_submatrix(): index out of range"); it_assert1(r2>=r1 && c2>=c1, "Mat::set_submatrix: r2::set_submatrix(): sizes don't match"); for (int i=0; i inline void Mat::set_submatrix(int r, int c, const Mat &m) { it_assert1(r>=0 && r+m.no_rows<=no_rows && c>=0 && c+m.no_cols<=no_cols, "Mat::set_submatrix(): index out of range"); for (int i=0; i inline void Mat::set_submatrix(int r1, int r2, int c1, int c2, const Num_T t) { if (r1 == -1) r1 = no_rows-1; if (r2 == -1) r2 = no_rows-1; if (c1 == -1) c1 = no_cols-1; if (c2 == -1) c2 = no_cols-1; it_assert1(r1>=0 && r2>=0 && r1=0 && c2>=0 && c1::set_submatrix(): i\ ndex out of range"); it_assert1(r2>=r1 && c2>=c1, "Mat::set_submatrix: r2 inline void Mat::del_row(int r) { it_assert1(r>=0 && r::del_row(): index out of range"); Mat Temp(*this); set_size(no_rows-1, no_cols, false); for (int i=0 ; i < r ; i++) { copy_vector(no_cols, &Temp.data[i], no_rows+1, &data[i], no_rows); } for (int i=r ; i < no_rows ; i++) { copy_vector(no_cols, &Temp.data[i+1], no_rows+1, &data[i], no_rows); } } template inline void Mat::del_rows(int r1, int r2) { it_assert1(r1>=0 && r2::del_rows(): index out of range"); for (int i=r1 ; i<=r2 ; i++) { del_row(r1); } // this should work, I don't understand why it does not. /* Mat Temp(*this); */ /* int n_deleted_rows = r2-r1+1; */ /* set_size(no_rows-n_deleted_rows, no_cols, false); */ /* for (int i=0 ; i < r1 ; i++) { */ /* copy_vector(no_cols, &Temp.data[i], Temp.no_rows, &data[i], no_rows); */ /* } */ /* for (int i=r2 ; i < no_rows ; i++) { */ /* copy_vector(no_cols, &Temp.data[i+1], Temp.no_rows, &data[i-n_deleted_rows+1], no_rows); */ /* } */ } template inline void Mat::del_col(int c) { it_assert1(c>=0 && c::del_col(): index out of range"); Mat Temp(*this); set_size(no_rows, no_cols-1, false); copy_vector(c*no_rows, Temp.data, data); copy_vector((no_cols - c)*no_rows, &Temp.data[(c+1)*no_rows], &data[c*no_rows]); } template inline void Mat::del_cols(int c1, int c2) { it_assert1(c1>=0 && c2::del_cols(): index out of range"); Mat Temp(*this); int n_deleted_cols = c2-c1+1; set_size(no_rows, no_cols-n_deleted_cols, false); copy_vector(c1*no_rows, Temp.data, data); copy_vector((no_cols-c1)*no_rows, &Temp.data[(c2+1)*no_rows], &data[c1*no_rows]); } template inline void Mat::ins_row(int r, const Vec &in) { it_assert1(r>=0 && r<=no_rows, "Mat::ins_row(): index out of range"); it_assert1((in.size() == no_cols) || no_rows==0, "Mat::ins_row(): vector size does not match"); if (no_cols==0) { no_cols = in.size(); } Mat Temp(*this); set_size(no_rows+1, no_cols, false); for (int i=0 ; i < r ; i++) { copy_vector(no_cols, &Temp.data[i], no_rows-1, &data[i], no_rows); } copy_vector(no_cols, in._data(), 1, &data[r], no_rows); for (int i=r+1 ; i < no_rows ; i++) { copy_vector(no_cols, &Temp.data[i-1], no_rows-1, &data[i], no_rows); } } template inline void Mat::ins_col(int c, const Vec &in) { it_assert1(c>=0 && c<=no_cols, "Mat::ins_col(): index out of range"); it_assert1(in.size() == no_rows || no_cols==0, "Mat::ins_col(): vector size does not match"); if (no_rows==0) { no_rows = in.size(); } Mat Temp(*this); set_size(no_rows, no_cols+1, false); copy_vector(c*no_rows, Temp.data, data); copy_vector(no_rows, in._data(), &data[c*no_rows]); copy_vector((no_cols-c-1)*no_rows, &Temp.data[c*no_rows], &data[(c+1)*no_rows]); } template inline void Mat::append_row(const Vec &in) { ins_row(no_rows, in); } template inline void Mat::append_col(const Vec &in) { ins_col(no_cols, in); } template Mat Mat::transpose() const { Mat temp(no_cols, no_rows); for (int i=0; i cmat cmat::hermitian_transpose() const; template Mat Mat::hermitian_transpose() const { Mat temp(no_cols, no_rows); for (int i=0; i Mat concat_horizontal(const Mat &m1, const Mat &m2) { it_assert1(m1.no_rows == m2.no_rows, "Mat::concat_horizontal; wrong sizes"); Mat temp(m1.no_rows, m1.no_cols+m2.no_cols); int i; for (i=0; i Mat concat_vertical(const Mat &m1, const Mat &m2) { it_assert1(m1.no_cols == m2.no_cols, "Mat::concat_vertical; wrong sizes"); Mat temp(m1.no_rows+m2.no_rows, m1.no_cols); int i; for (i=0; i inline void Mat::operator=(Num_T t) { for (int i=0; i inline void Mat::operator=(const Mat &m) { set_size(m.no_rows,m.no_cols, false); if (m.datasize == 0) return; copy_vector(m.datasize, m.data, data); } template inline void Mat::operator=(const Vec &v) { it_assert1( (no_rows == 1 && no_cols == v.size()) || (no_cols == 1 && no_rows == v.size()), "Mat::op=(vec); wrong size"); // construct a 1-d column of the vector set_size(v.size(), 1, false); set_col(0, v); } //-------------------- Templated friend functions -------------------------- template inline void Mat::operator+=(const Mat &m) { if (datasize == 0) operator=(m); else { int i, j, m_pos=0, pos=0; it_assert1(m.no_rows==no_rows && m.no_cols==no_cols,"Mat::operator+=: wrong sizes"); for (i=0; i inline Mat operator+(const Mat &m1, const Mat &m2) { Mat r(m1.no_rows, m1.no_cols); int i, j, m1_pos=0, m2_pos=0, r_pos=0; it_assert1(m1.no_rows==m2.no_rows && m1.no_cols == m2.no_cols, "Mat::operator+: wrong sizes"); for (i=0; i inline void Mat::operator+=(Num_T t) { for (int i=0; i inline Mat operator+(const Mat &m, Num_T t) { Mat r(m.no_rows, m.no_cols); for (int i=0; i inline Mat operator+(Num_T t, const Mat &m) { Mat r(m.no_rows, m.no_cols); for (int i=0; i inline void Mat::operator-=(const Mat &m) { int i,j, m_pos=0, pos=0; if (datasize == 0) { set_size(m.no_rows, m.no_cols, false); for (i=0; i::operator-=: wrong sizes"); for (i=0; i inline Mat operator-(const Mat &m1, const Mat &m2) { Mat r(m1.no_rows, m1.no_cols); int i, j, m1_pos=0, m2_pos=0, r_pos=0; it_assert1(m1.no_rows==m2.no_rows && m1.no_cols == m2.no_cols, "Mat::operator-: wrong sizes"); for (i=0; i inline void Mat::operator-=(Num_T t) { for (int i=0; i inline Mat operator-(const Mat &m, Num_T t) { Mat r(m.no_rows, m.no_cols); int i, j, m_pos=0, r_pos=0; for (i=0; i inline Mat operator-(Num_T t, const Mat &m) { Mat r(m.no_rows, m.no_cols); int i, j, m_pos=0, r_pos=0; for (i=0; i inline Mat operator-(const Mat &m) { Mat r(m.no_rows, m.no_cols); int i, j, m_pos=0, r_pos=0; for (i=0; i void mat::operator*=(const mat &m); template<> void cmat::operator*=(const cmat &m); #endif //NO_CBLAS template inline void Mat::operator*=(const Mat &m) { it_assert1(no_cols == m.no_rows,"Mat::operator*=: wrong sizes"); Mat r(no_rows, m.no_cols); // this consumes unnecessary memory Num_T tmp; int i,j,k, r_pos=0, pos=0, m_pos=0; for (i=0; i inline void Mat::operator*=(Num_T t) { for (int i=0; i mat operator*(const mat &m1, const mat &m2); template<> cmat operator*(const cmat &m1, const cmat &m2); #endif //NO_CBLAS template Mat operator*(const Mat &m1, const Mat &m2) { it_assert1(m1.no_cols == m2.no_rows,"Mat::operator*: wrong sizes"); Mat r(m1.no_rows, m2.no_cols); Num_T tmp; int i, j, k; Num_T *tr=r.data, *t1, *t2=m2.data; for (i=0; i