A lite library for operations in linear algebra

Embed: (wiki syntax)

« Back to documentation index

Show/hide line numbers MATRIX_LIB.cpp Source File

MATRIX_LIB.cpp

00001 #include "MATRIX_LIB.h"
00002 
00003 // No pre-assignments, empty matrix
00004 Mat::Mat(void):
00005         _nRow(0), _nCol(0),
00006         data(0, vector<float>())
00007 {
00008 
00009 }
00010 // Pre-assign sizes
00011 Mat::Mat(size_t nRow_in, size_t nCol_in):
00012         _nRow(nRow_in), _nCol(nCol_in),
00013         data(nRow_in, vector<float>(nCol_in))
00014 {
00015 
00016 }
00017 // Pre-assign sizes and elements
00018 Mat::Mat(size_t nRow_in, size_t nCol_in, float element_in):
00019         _nRow(nRow_in), _nCol(nCol_in),
00020         data(nRow_in, vector<float>(nCol_in, element_in))
00021 {
00022 
00023 }
00024 
00025 // Public methods
00026 //----------------------------------//
00027 // Size ---
00028 // "Get" the size of the matrix
00029 size_t Mat::nRow() {
00030     return this->_nRow;
00031 }
00032 size_t Mat::nCol() {
00033     return this->_nCol;
00034 }
00035 size_t Mat::nElement() { // Number of elements
00036     return (_nRow*_nCol);
00037 }
00038 // const functions
00039 size_t Mat::nRow() const {
00040     return this->_nRow;
00041 }
00042 size_t Mat::nCol() const {
00043     return this->_nCol;
00044 }
00045 size_t Mat::nElement() const { // Number of elements
00046     return (_nRow*_nCol);
00047 }
00048 // "Set" and "get" the new size of the matrix
00049 size_t Mat::nRow(size_t nRow_in){
00050     _nRow = nRow_in;
00051     // Resize the matrix
00052     data.resize(_nRow, vector<float>(_nCol));
00053     return _nRow;
00054 }
00055 size_t Mat::nCol(size_t nCol_in){
00056     _nCol = nCol_in;
00057     // Resize the matrix
00058     for (size_t i = 0; i < _nRow; ++i){
00059         data[i].resize(_nCol);
00060     }
00061     return _nCol;
00062 }
00063 // "Set" and "get" the new size of the matrix with assignment to new elements
00064 size_t Mat::nRow(size_t nRow_in, float element_in){
00065     _nRow = nRow_in;
00066     // Resize the matrix
00067     data.resize(_nRow, vector<float>(_nCol, element_in));
00068     return _nRow;
00069 }
00070 size_t Mat::nCol(size_t nCol_in, float element_in){
00071     _nCol = nCol_in;
00072     // Resize the matrix
00073     for (size_t i = 0; i < _nRow; ++i){
00074         data[i].resize(_nCol, element_in);
00075     }
00076     return _nCol;
00077 }
00078 // Resize on both row and column
00079 void Mat::resize(size_t nRow_in, size_t nCol_in){
00080     // Check the size
00081     if (_nRow == nRow_in && _nCol == nCol_in){
00082         return; // Nothing to do
00083     }
00084 
00085     // Change the nCol first to resize the old row vector
00086     _nCol = nCol_in;
00087     // Resize the matrix
00088     for (size_t i = 0; i < _nRow; ++i){
00089         data[i].resize(_nCol);
00090     }
00091     // Then extend the row
00092     _nRow = nRow_in;
00093     data.resize(_nRow, vector<float>(_nCol));
00094 }
00095 void Mat::resize(size_t nRow_in, size_t nCol_in, float element_in){ // Assign the new element_in to new elements
00096     // Check the size
00097     if (_nRow == nRow_in && _nCol == nCol_in){
00098         return; // Nothing to do
00099     }
00100 
00101     // Change the nCol first to resize the old row vector
00102     _nCol = nCol_in;
00103     // Resize the matrix
00104     for (size_t i = 0; i < _nRow; ++i){
00105         data[i].resize(_nCol, element_in);
00106     }
00107     // Then extend the row
00108     _nRow = nRow_in;
00109     data.resize(_nRow, vector<float>(_nCol, element_in));
00110 }
00111 // Reshape the matrix, keep all the elements but change the shpae of the matrix (eg. 1 by 6 -> 2 by 3 or 3 by 2)
00112 void Mat::reshape(size_t nRow_in, size_t nCol_in, bool is_Row_first){ // If is_Row_first, keep the elements' order the same in the row direction; otherwise, keep the order the same in column direction.
00113     // Note: nRow_in * nCol_in should be better be the same as this->nElement(), or some elemnets will be lost or be added (add zeros, actually)
00114     //
00115     Mat M_new(nRow_in, nCol_in);
00116     // The index for the old matrix
00117     size_t io = 0;
00118     size_t jo = 0;
00119     //
00120     if (is_Row_first){ // Row go first
00121         //
00122         for (size_t i = 0; i < nRow_in; ++i){
00123             for (size_t j = 0; j < nCol_in; ++j){ // Row go first
00124                 //-------------------//
00125                 // In-range check:
00126                 if (io < this->_nRow && jo < this->_nCol){ // In range
00127                     M_new.data[i][j] = this->data[io][jo];
00128                     // Row go first
00129                     jo++;
00130                     // Change row
00131                     if (jo >= this->_nCol){
00132                         jo = 0;
00133                         io++;
00134                     }
00135                     //
00136                 }else{ // Out of range
00137                     M_new.data[i][j] = 0.0;
00138                 }
00139                 //-------------------//
00140             }
00141         }
00142     }else{ // Column go first
00143         //
00144         for (size_t j = 0; j < nCol_in; ++j){
00145             for (size_t i = 0; i < nRow_in; ++i){ // Column go first
00146                 //-------------------//
00147                 // In-range check:
00148                 if (io < this->_nRow && jo < this->_nCol){ // In range
00149                     M_new.data[i][j] = this->data[io][jo];
00150                     // Column go fist
00151                     io++;
00152                     // Change row
00153                     if (io >= this->_nRow){
00154                         io = 0;
00155                         jo++;
00156                     }
00157                     //
00158                 }else{ // Out of range
00159                     M_new.data[i][j] = 0.0;
00160                 }
00161                 //-------------------//
00162             }
00163         }
00164     }
00165     //
00166     // Update this matrix
00167     *this = M_new;
00168 }
00169 // Synchronize the _nRow and _nCol with the true size of data
00170 void Mat::syncSizeInData(){
00171     //
00172     this->_nRow = this->data.size();
00173     //
00174     if (this->_nRow == 0){ // Empty
00175         this->_nCol = 0;
00176     }else{
00177         this->_nCol = this->data[0].size();
00178     }
00179 }
00180 
00181 // Assignment ---
00182 void Mat::assign(std::stringstream &ss_in, size_t nRow_in, size_t nCol_in){ // The stringstream ss_in contains a nRow_in by nCol_in matrix.
00183     // Resize
00184     this->resize(nRow_in, nCol_in);
00185     //
00186     for (size_t i = 0; i < this->_nRow; ++i){
00187         for (size_t j = 0; j < this->_nCol; ++j){
00188             ss_in >> this->data[i][j];
00189         }
00190     }
00191 }
00192 // From primitive array in c++
00193 void Mat::assign(float* Matrix_in, size_t nRow_in, size_t nCol_in){ // From primitive 1-D array in c++, Matrix_in is a nRow_in by nCol_in matrix.
00194     // Resize
00195     this->resize(nRow_in, nCol_in);
00196     //
00197     for (size_t i = 0; i < this->_nRow; ++i){
00198         for (size_t j = 0; j < this->_nCol; ++j){
00199             // this->data[i][j] = Matrix_in[i][j];
00200             this->data[i][j] = *Matrix_in;
00201             Matrix_in++;
00202         }
00203     }
00204 }
00205 void Mat::assign(const vector<float> &vec_in, bool is_RowVec){ // From 1-D vector
00206     size_t n = vec_in.size();
00207     //
00208     if (is_RowVec){ // Row vector
00209         //
00210         this->resize(1,n);
00211         //
00212         for (size_t j = 0; j < _nCol; ++j){
00213             this->data[0][j] = vec_in[j];
00214         }
00215 
00216         //
00217     }else{ // Column vector
00218         //
00219         this->resize(n,1);
00220         //
00221         for (size_t i = 0; i < _nRow; ++i){
00222             this->data[i][0] = vec_in[i];
00223         }
00224         //
00225     }
00226 }
00227 void Mat::assign(const vector<vector<float> > &MatData_in){ // A = MatData_in, assign a primitive matrix directly
00228     this->data = MatData_in;
00229     this->syncSizeInData();
00230 }
00231 // Partial asignment
00232 void Mat::setPart(const Mat &M_part, size_t m_from, size_t n_from){ // The block starts from (m_from, n_from)
00233     size_t m_to, n_to;
00234     // The index of the last element of the inserted block
00235     m_to = m_from + M_part.nRow() - 1;
00236     n_to = n_from + M_part.nCol() - 1;
00237 
00238     //-------------------------------//
00239     // Note: if the region of the insertion block exceed the region of this matrix,
00240     //       resize to a biger one and insert zero at the blank sides.
00241 
00242     // The size of the new matrix
00243     size_t M,N;
00244     //
00245     M = this->_nRow;
00246     N = this->_nCol;
00247     //
00248     // Check if block exceeds the region
00249     if (m_to >= M ){
00250         M = m_to + 1;
00251     }
00252     if(n_to >= N){
00253         N = n_to + 1;
00254     }
00255     // Deside if the resize action is needed
00256     this->resize(M, N, 0.0);
00257 
00258     // Start the insertion
00259     //
00260     // Initialization
00261     size_t ip = 0;
00262     size_t jp = 0;
00263     for (size_t i = m_from; i <= m_to; ++i ){
00264         jp = 0;
00265         for (size_t j = n_from; j <= n_to; ++j){
00266             // M_part.data[ip][jp] = data[i][j]; // <- get
00267             this->data[i][j] = M_part.data[ip][jp]; // <- insert
00268             //
00269             jp++;
00270         }
00271         ip++;
00272     }
00273 
00274 }
00275 // Spetial assignments
00276 void Mat::zeros(size_t nRow_in, size_t nCol_in){ // zeros(m,n)
00277     resize(0,0);
00278     resize(nRow_in,nCol_in, 0.0);
00279 }
00280 void Mat::ones(size_t nRow_in, size_t nCol_in){ // ones(m,n)
00281     resize(0,0);
00282     resize(nRow_in,nCol_in, 1.0);
00283 }
00284 void Mat::eye(size_t nRow_in, size_t nCol_in){ // Identity matrix, eye(m,n)
00285     zeros(nRow_in, nCol_in);
00286     //
00287     size_t K;
00288     //
00289     // min()
00290     if (nRow_in < nCol_in){
00291         K = nRow_in;
00292     }else{
00293         K = nCol_in;
00294     }
00295     //
00296     for (size_t i = 0; i < K; ++i){
00297         data[i][i] = 1.0;
00298     }
00299     //
00300 }
00301 void Mat::diag(const Mat &M_diag){ // Transform the row/column vector to a diagonal matrix and assign into this matrix
00302     size_t M = M_diag.nRow();
00303     size_t N = M_diag.nCol();
00304     size_t nE = M_diag.nElement();
00305     //
00306     this->zeros(nE,nE); // nE by nE matrix
00307     //
00308     size_t id = 0;
00309     size_t jd = 0;
00310     for (size_t i = 0; i < nE; ++i){
00311         this->data[i][i] = M_diag.data[id][jd];
00312         jd++;
00313         if (jd >= N){
00314             jd = 0;
00315             id++;
00316         }
00317     }
00318 }
00319 
00320 
00321 // Get elements ---
00322 float Mat::at(size_t m, size_t n){ // Get the element at (m,n)
00323     return data[m][n];
00324 }
00325 Mat Mat::getPart(size_t m_from, size_t m_to, size_t n_from, size_t n_to){ // Get partial, M_part = M(m_from:m_to, n_from:n_to)
00326     int M = m_to - m_from + 1;
00327     int N = n_to - n_from + 1;
00328     //
00329     if (M < 0 || N < 0){
00330         Mat Empty;
00331         return Empty;
00332     }
00333     //
00334     Mat M_part(M,N);
00335 
00336     //
00337     // Initialization
00338     size_t ip = 0;
00339     size_t jp = 0;
00340     //
00341     for (size_t i = m_from; i <= m_to; ++i ){
00342         jp = 0;
00343         for (size_t j = n_from; j <= n_to; ++j){
00344             // In-range check
00345             if (i < this->_nRow && j < this->_nCol){
00346                 M_part.data[ip][jp] = this->data[i][j];
00347             }else{ // Out of bound, set elements to zeros
00348                 M_part.data[ip][jp] = 0.0;
00349             }
00350             //
00351             jp++;
00352         }
00353         ip++;
00354     }
00355     return M_part;
00356 }
00357 Mat Mat::getCol(size_t n){ // Get a specific column vector in 2-D matrix form
00358     return getPart(0, (this->_nRow-1), n, n);
00359 }
00360 Mat Mat::getRow(size_t m){ // Get a specific row vector in 2-D matrix form
00361     return getPart(m, m, 0, (this->_nCol-1));
00362 }
00363 vector<float> Mat::getColVec(size_t n){ // Return a c++ vector, M(:,n)
00364     vector<float> ColVec(this->_nRow);
00365     for (size_t i = 0; i < this->_nRow; ++i){
00366         ColVec[i] = this->data[i][n];
00367     }
00368     return ColVec;
00369 }
00370 vector<float> Mat::getRowVec(size_t m){ // Return a c++ vector, M(m,:)
00371     return this->data[m];
00372 }
00373 
00374 // Print out matrix ---
00375 std::string Mat::print(void){ // Print this matrix out as string
00376     /*
00377     std::string str_out;
00378     //
00379     std::string endl("\n"); // The end line character
00380 
00381     //
00382     str_out += endl;
00383     for (size_t i = 0; i < this->_nRow; ++i){
00384         for (size_t j = 0; j < this->_nCol; ++j){
00385             str_out += std::to_string( this->data[i][j] ) + "\t";
00386         }
00387         str_out += endl;
00388     }
00389     str_out += endl;
00390     //
00391     return str_out;
00392     */
00393 
00394     // Using stringstream
00395     std::stringstream ss_out;
00396     //
00397     ss_out << std::endl;
00398     //
00399     for (size_t i = 0; i < this->_nRow; ++i){
00400         for (size_t j = 0; j < this->_nCol; ++j){
00401             ss_out << this->data[i][j] << "\t";
00402         }
00403         ss_out << std::endl;
00404     }
00405     ss_out << std::endl;
00406     //
00407     return ss_out.str();
00408 }
00409 
00410 
00411 // Operations =====
00412 // Comparison
00413 bool Mat::is_equal(const Mat &M_in){
00414     //
00415     if (M_in.nRow() != this->_nRow || M_in.nCol() != this->_nCol){
00416         return false;
00417     }
00418 
00419     // They are of the same size.
00420     for (size_t i = 0; i < this->_nRow; ++i){
00421         for (size_t j = 0; j < this->_nCol; ++j){
00422             if ( this->data[i][j] != M_in.data[i][j] ){
00423                 return false;
00424             }
00425         }
00426     }
00427     return true;
00428 }
00429 
00430 // Self-operation (affecting on this matrix) ---
00431 void Mat::scaleUp(float scale){ // M *= scale
00432     //
00433     for (size_t i = 0; i < _nRow; ++i){
00434         MP::Get_VectorScaleUp(data[i], scale);
00435     }
00436 }
00437 void Mat::scaleUp_Mat(const Mat &M_right){ // M = M.times(M_right), element-wise multiplication
00438     //
00439     for (size_t i = 0; i < _nRow; ++i){
00440         for (size_t j = 0; j < _nCol; ++j){
00441             data[i][j] *= M_right.data[i][j];
00442         }
00443     }
00444 }
00445 void Mat::increase(float scale){ // M += scale, for each element in M
00446     //
00447     for (size_t i = 0; i < _nRow; ++i){
00448         for (size_t j = 0; j < _nCol; ++j){
00449             data[i][j] += scale;
00450         }
00451     }
00452 }
00453 void Mat::decrease(float scale){ // M -= scale, for each element in M
00454     //
00455     for (size_t i = 0; i < _nRow; ++i){
00456         for (size_t j = 0; j < _nCol; ++j){
00457             data[i][j] -= scale;
00458         }
00459     }
00460 }
00461 void Mat::increase(const Mat &M_right){ // M += M_right
00462     //
00463     for (size_t i = 0; i < _nRow; ++i){
00464         MP::Get_VectorIncrement(data[i], M_right.data[i], false); // +=
00465     }
00466 }
00467 void Mat::decrease(const Mat &M_right){ // M -= M_right
00468     //
00469     for (size_t i = 0; i < _nRow; ++i){
00470         MP::Get_VectorIncrement(data[i], M_right.data[i], true); // -=
00471     }
00472 }
00473 // Single-operation, output another matrix ---
00474 // Transpose
00475 Mat& Mat::T(void){ // Return the transposed version of this matrix
00476     //
00477     // Mat M_t(_nCol, _nRow);
00478     static Mat M_t;
00479     M_t.resize(_nCol, _nRow);
00480     //
00481     for (size_t i = 0; i < _nRow; ++i){
00482         for (size_t j = 0; j < _nCol; ++j){
00483             M_t.data[j][i] = data[i][j];
00484         }
00485     }
00486     return M_t;
00487 }
00488 // Inverse
00489 Mat& Mat::inverse(void){ // Return the inversion of this matrix
00490     // The output matrix
00491     static Mat M_inv;
00492 
00493     if (this->_nRow != this->_nCol){
00494         M_inv.resize(0,0); // Empty matrix
00495         return M_inv;
00496     }
00497 
00498 
00499     // M_inv.resize( _nCol, _nRow); // Maybe we will impliment the psodo inverse next time
00500 
00501     //
00502     M_inv.data = MP::MatrixInversion(this->data);
00503 
00504     // Update the _nRow and _nCol
00505     M_inv.syncSizeInData();
00506 
00507     return M_inv;
00508 }
00509 Mat& Mat::intPower(int power){ // M^power, M^0 -> I, M^-1 -> M_inverse
00510     // The output matrix
00511     static Mat M_out;
00512 
00513     // Spetial case M^1 -> M
00514     if (power == 1){
00515         M_out = *this;
00516         return M_out;
00517     }
00518     // Spetial case M^0 -> I
00519     if (power == 0){
00520         M_out.eye(_nRow,_nCol);
00521         return M_out;
00522     }
00523 
00524     // Check for square matrix
00525     if (this->_nRow != this->_nCol){
00526         M_out.resize(0,0); // Error, return empty matrix
00527         return M_out;
00528     }
00529     //-----------------------------------------------//
00530 
00531     // Spetial case M^-1 -> this->inverse()
00532     if (power == -1){
00533         // M^-1, special case for speed up the frequent useage
00534         return this->inverse();
00535     }
00536     //
00537 
00538     // Other cases
00539     //
00540     if (power > 1){ // power > 1
00541         M_out = *this;
00542         for (size_t i = 0; i < (power-1); ++i){
00543             M_out = this->dot(M_out);
00544         }
00545     }else if (power < -1){ // power < -1
00546         Mat M_inv = this->inverse();
00547         M_out = M_inv; // this->inverse()
00548         for (size_t i = 0; i < (-power - 1); ++i){
00549             M_out = M_inv.dot(M_out);
00550         }
00551     }else{
00552         // Nothing, power \in {-1, 0, 1} has been processed
00553     }
00554     //
00555     return M_out;
00556 
00557 }
00558 
00559 // Duo-operation, output another matrix ---
00560 // Plus
00561 Mat& Mat::plus(float scale){ // (M + scale), for each element in M
00562     // Create the output matrix
00563     // Mat M_out(_nRow,_nCol);
00564     static Mat M_out;
00565     M_out.resize(_nRow,_nCol);
00566     //
00567     for (size_t i = 0; i < _nRow; ++i){
00568         for (size_t j = 0; j < _nCol; ++j){
00569             M_out.data[i][j] = this->data[i][j] + scale;
00570         }
00571     }
00572     return M_out;
00573 }
00574 Mat& Mat::minus(float scale){ // (M - scale), for each element in M
00575     // Create the output matrix
00576     // Mat M_out(_nRow,_nCol);
00577     static Mat M_out;
00578     M_out.resize(_nRow,_nCol);
00579     //
00580     for (size_t i = 0; i < _nRow; ++i){
00581         for (size_t j = 0; j < _nCol; ++j){
00582             M_out.data[i][j] = this->data[i][j] - scale;
00583         }
00584     }
00585     return M_out;
00586 }
00587 Mat& Mat::minus(float scale, bool is_reversed){ // is_reversed -> (scale - M), for each element in M
00588     // Create the output matrix
00589     // Mat M_out(_nRow,_nCol);
00590     static Mat M_out;
00591     M_out.resize(_nRow,_nCol);
00592     //
00593     if (is_reversed){ // Reversed
00594         for (size_t i = 0; i < _nRow; ++i){
00595             for (size_t j = 0; j < _nCol; ++j){
00596                 M_out.data[i][j] = scale - this->data[i][j]; // (scale - M)
00597             }
00598         }
00599     }else{
00600         for (size_t i = 0; i < _nRow; ++i){
00601             for (size_t j = 0; j < _nCol; ++j){
00602                 M_out.data[i][j] = this->data[i][j] - scale; // (M - scale)
00603             }
00604         }
00605     }
00606     return M_out;
00607 }
00608 Mat& Mat::plus(const Mat &M_right){
00609     // Create the output matrix
00610     // Mat M_out(_nRow,_nCol);
00611     static Mat M_out;
00612     M_out.resize(_nRow,_nCol);
00613     //
00614     for (size_t i = 0; i < _nRow; ++i){
00615         for (size_t j = 0; j < _nCol; ++j){
00616             M_out.data[i][j] = this->data[i][j] + M_right.data[i][j];
00617         }
00618     }
00619     return M_out;
00620 }
00621 Mat& Mat::minus(const Mat &M_right){
00622     // Create the output matrix
00623     // Mat M_out(_nRow,_nCol);
00624     static Mat M_out;
00625     M_out.resize(_nRow,_nCol);
00626     //
00627     for (size_t i = 0; i < _nRow; ++i){
00628         for (size_t j = 0; j < _nCol; ++j){
00629             M_out.data[i][j] = this->data[i][j] - M_right.data[i][j];
00630         }
00631     }
00632     return M_out;
00633 }
00634 // Concatenation
00635 Mat Mat::cat_below(const Mat &M_in){ // Below this matrix, [A; B]
00636     Mat M_cat;
00637     // Using setPart(), the size will be automatically changed.
00638     //------------------------------//
00639     // Put this matrix at (0, 0)
00640     M_cat.setPart(*this, 0, 0);
00641     // Put M_in at (_nRow, 0)
00642     M_cat.setPart(M_in, this->_nRow, 0);
00643     //------------------------------//
00644     return M_cat;
00645 }
00646 Mat Mat::cat_right(const Mat &M_in){ // Right-side of this matrix, [A, B]
00647     Mat M_cat;
00648     // Using setPart(), the size will be automatically changed.
00649     //------------------------------//
00650     // Put this matrix at (0, 0)
00651     M_cat.setPart(*this, 0, 0);
00652     // Put M_in at (0, _nCol)
00653     M_cat.setPart(M_in, 0, this->_nCol);
00654     //------------------------------//
00655     return M_cat;
00656 }
00657 Mat Mat::cat(const Mat &M_in, bool is_horizontal){ // is_horizontal --> cat_Right(); otherwise --> cat_Below()
00658     //
00659     if(is_horizontal){ // [A, B]
00660         return this->cat_right(M_in);
00661     }else{ // [A; B]
00662         return this->cat_below(M_in);
00663     }
00664 }
00665 
00666 // Scalar/Element-wise multiplication
00667 Mat& Mat::times(float scale){ // Scalar multiplication
00668     // Create the output matrix
00669     // Mat M_out(_nRow,_nCol);
00670     //
00671     static Mat M_out;
00672     M_out.resize(_nRow, _nCol);
00673 
00674 
00675     if(scale == -1.0){
00676         //
00677         for (size_t i = 0; i < _nRow; ++i){
00678             for (size_t j = 0; j < _nCol; ++j){
00679                 M_out.data[i][j] = -(this->data[i][j]);
00680             }
00681         }
00682     }else{
00683         //
00684         for (size_t i = 0; i < _nRow; ++i){
00685             for (size_t j = 0; j < _nCol; ++j){
00686                 M_out.data[i][j] = this->data[i][j] * scale;
00687             }
00688         }
00689     }
00690     return M_out;
00691 }
00692 Mat& Mat::times(const Mat &M_right){ // Element-wise multiplication
00693     // Create the output matrix
00694     // Mat M_out(_nRow,_nCol);
00695     static Mat M_out;
00696     M_out.resize(_nRow,_nCol);
00697     //
00698     for (size_t i = 0; i < _nRow; ++i){
00699         for (size_t j = 0; j < _nCol; ++j){
00700             M_out.data[i][j] = this->data[i][j] * M_right.data[i][j];
00701         }
00702     }
00703     return M_out;
00704 }
00705 
00706 
00707 // Matrix multiplication
00708 // Note: this matrix is the "left" one
00709 Mat& Mat::dot(const Mat &M_right){ // Similar to the nomenclature of numpy in Python
00710 
00711     // Size check
00712     // mxn = (mxk)*(kxn)
00713     size_t M, K, N;
00714     //
00715     M = _nRow;
00716     K = _nCol;
00717     N = M_right.nCol();
00718 
00719     // The output matrix
00720     // Mat M_out(M, N);
00721     static Mat M_out;
00722     M_out.resize(M, N);
00723 
00724 
00725     // cout << "here in\n";
00726     float* ptr_data = NULL;
00727     // Using indexing
00728     for (size_t i = 0; i < M; ++i){ // row in m_out
00729         //
00730         // cout << "#i: " << i << "\n";
00731         //
00732         for (size_t j = 0; j < N; ++j){ // column in m_out
00733             //
00734             // cout << "#j: " << j << "\n";
00735             //
00736             ptr_data = &(M_out.data[i][j]);
00737             //
00738             *ptr_data = 0.0;
00739             for (size_t k = 0; k < K; ++k){
00740                 if (data[i][k] != 0.0 && M_right.data[k][j] != 0.0)
00741                     *ptr_data += data[i][k]*M_right.data[k][j];
00742             }
00743         }
00744     }
00745 
00746     return M_out;
00747 }
00748 Mat& Mat::dot(bool Transpose_left, const Mat &M_right, bool Transpose_right){ // Extended version for conbining the ability of transpose of both mtrices
00749 
00750     // Size check
00751     // mxn = (mxk)*(kxn)
00752     size_t M,K,N;
00753     // Left transpose
00754     if (Transpose_left){
00755         M = _nCol;
00756         K = _nRow;
00757     }else{
00758         M = _nRow;
00759         K = _nCol;
00760     }
00761     // Right transpose
00762     if (Transpose_right){
00763         N = M_right.nRow();
00764     }else{
00765         N = M_right.nCol();
00766     }
00767 
00768     // The output matrix
00769     // Mat M_out(M, N);
00770     static Mat M_out;
00771     M_out.resize(M, N);
00772 
00773     //
00774     float* ptr_data = NULL;
00775     //
00776 
00777     // Check the conditions of transpotations
00778     if(Transpose_left){
00779         if(Transpose_right){ // Both transposed
00780             //
00781             // Using indexing
00782             for (size_t i = 0; i < M; ++i){ // row in m_out
00783                 for (size_t j = 0; j < N; ++j){ // column in m_out
00784                     // M_out.data[i][j] = 0.0;
00785                     ptr_data = &(M_out.data[i][j]);
00786                     *ptr_data = 0.0;
00787                     //
00788                     for (size_t k = 0; k < K; ++k){
00789                         if (data[k][i] != 0.0 && M_right.data[j][k] != 0.0) // (i,k) -> (k,i), (k,j) -> (j,k)
00790                             *ptr_data += data[k][i]*M_right.data[j][k]; // (i,k) -> (k,i), (k,j) -> (j,k)
00791                     }
00792                 }
00793             }
00794             //
00795         }else{ // Only left transpose
00796             //
00797             // Using indexing
00798             for (size_t i = 0; i < M; ++i){ // row in m_out
00799                 for (size_t j = 0; j < N; ++j){ // column in m_out
00800                     // M_out.data[i][j] = 0.0;
00801                     ptr_data = &(M_out.data[i][j]);
00802                     *ptr_data = 0.0;
00803                     //
00804                     for (size_t k = 0; k < K; ++k){
00805                         if (data[k][i] != 0.0 && M_right.data[k][j] != 0.0) // (i,k) -> (k,i)
00806                             *ptr_data += data[k][i]*M_right.data[k][j]; // (i,k) -> (k,i)
00807                     }
00808                 }
00809             }
00810             //
00811         }
00812     }else{
00813         if(Transpose_right){ // Only right transpose
00814             //
00815             // Using indexing
00816             for (size_t i = 0; i < M; ++i){ // row in m_out
00817                 for (size_t j = 0; j < N; ++j){ // column in m_out
00818                     // M_out.data[i][j] = 0.0;
00819                     ptr_data = &(M_out.data[i][j]);
00820                     *ptr_data = 0.0;
00821                     //
00822                     for (size_t k = 0; k < K; ++k){
00823                         if (data[i][k] != 0.0 && M_right.data[j][k] != 0.0) // (k,j) -> (j,k)
00824                             *ptr_data += data[i][k]*M_right.data[j][k]; // (k,j) -> (j,k)
00825                     }
00826                 }
00827             }
00828             //
00829         }else{ // Normal
00830             //
00831             // Using indexing
00832             for (size_t i = 0; i < M; ++i){ // row in m_out
00833                 for (size_t j = 0; j < N; ++j){ // column in m_out
00834                     // M_out.data[i][j] = 0.0;
00835                     ptr_data = &(M_out.data[i][j]);
00836                     *ptr_data = 0.0;
00837                     //
00838                     for (size_t k = 0; k < K; ++k){
00839                         if (data[i][k] != 0.0 && M_right.data[k][j] != 0.0)
00840                             *ptr_data += data[i][k]*M_right.data[k][j];
00841                     }
00842                 }
00843             }
00844             //
00845         }
00846     }
00847 
00848     return M_out;
00849 }
00850 // Operator overloading
00851 //---------------------------------------//
00852 // Member
00853 // A <- b, assignment
00854 Mat& Mat::operator = (float scale){ // Assign a real number as 1 by 1 matrix
00855     // this->resize(1,1, scale);
00856     this->resize(1,1);
00857     this->data[0][0] = scale;
00858     return *this;
00859 }
00860 Mat& Mat::operator = (const vector<float> &colVec_in){ // A = vec_in, assign as a column vector
00861     this->assign(colVec_in, false); // A column vector
00862     return *this;
00863 }
00864 Mat& Mat::operator = (const vector<vector<float> > &MatData_in){ // A = MatData_in, assign a primitive matrix directly
00865     this->assign(MatData_in);
00866     return *this;
00867 }
00868 // b <- A, get value, using implicit type conversion
00869 // Note: no implicit conversion to float, since it would be ambiguous in other operator overriding
00870 Mat::operator float (){ // Get the first element as float, good for a 1 by 1 matrix
00871     return this->data[0][0];
00872 }
00873 Mat::operator std::vector<float> (){ // Get the column vector as std::vector<float> in c++
00874     return this->getColVec(0);
00875 }
00876 // A[]
00877 vector<float>& Mat::operator [] (size_t m){ // Indexing the m-th row, equivalent to data[m]
00878     return this->data[m];
00879 }
00880 // A == B
00881 bool Mat::operator == (Mat const& B){ // is_equal()
00882     return this->is_equal(B);
00883 }
00884 // A^z, z \in Z
00885 Mat& Mat::operator ^ (int power){ // A^power, this->intPower()
00886     return this->intPower(power);
00887 }
00888 
00889 
00890 // Non-member
00891 //
00892 // -A
00893 Mat& operator - (const Mat &A){
00894     return Mat(A).times(-1.0);
00895 }
00896 // A + B
00897 Mat& operator + (float a, const Mat &B){
00898     return Mat(B).plus(a);
00899 }
00900 Mat& operator + (const Mat &A, float b){
00901     return Mat(A).plus(b);
00902 }
00903 Mat& operator + (const Mat &A, const Mat &B){
00904     return Mat(A).plus(B);
00905 }
00906 // A - B
00907 Mat& operator - (float a, const Mat &B){
00908     return Mat(B).minus(a, true); // (a - B)
00909 }
00910 Mat& operator - (const Mat &A, float b){
00911     return Mat(A).minus(b);
00912 }
00913 Mat& operator - (const Mat &A, const Mat &B){
00914     return Mat(A).minus(B);
00915 }
00916 // A * B
00917 Mat& operator * (float a, const Mat &B){
00918     return Mat(B).times(a);
00919 }
00920 Mat& operator * (const Mat &A, float b){
00921     return Mat(A).times(b);
00922 }
00923 Mat& operator * (const Mat &A, const Mat &B){
00924     // Matrix multiplication
00925     return Mat(A).dot(B);
00926 }
00927 // A/B, including matrix inversion (eg, (1.0/B) <--> B^-1)
00928 Mat& operator / (float a, const Mat &B){ // a*(B^-1), for that B to be inversed
00929     if (a == 1.0){
00930         return Mat(B).inverse();
00931     }
00932     return ( ( Mat(B).inverse()).times(a) );
00933 }
00934 Mat& operator / (const Mat &A, float b){ // A*(1/b), scalar multiplication of 1/b
00935     return Mat(A).times(1.0/b);
00936 }
00937 Mat& operator / (const Mat &A, const Mat &B){ // A*(B^-1), for that B to be inversed
00938     return ( Mat(A).dot(Mat(B).inverse()) );
00939 }
00940 // A += B
00941 Mat& operator += (Mat &A, float b){
00942     A.increase(b);
00943     return A;
00944 }
00945 Mat& operator += (Mat &A, const Mat &B){
00946     A.increase(B);
00947     return A;
00948 }
00949 // A -= B
00950 Mat& operator -= (Mat &A, float b){
00951     A.decrease(b);
00952     return A;
00953 }
00954 Mat& operator -= (Mat &A, const Mat &B){
00955     A.decrease(B);
00956     return A;
00957 }
00958 // A *= B
00959 Mat& operator *= (Mat &A, float b){
00960     A.scaleUp(b);
00961     return A;
00962 }
00963 Mat& operator *= (Mat &A, const Mat &B){ // Matrix multiplication
00964     A = A.dot(B);
00965     return A;
00966 }
00967 
00968 //---------------------------------------//
00969 // end Operator overloading