A lite library for operations in linear algebra
Embed:
(wiki syntax)
Show/hide line numbers
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
Generated on Fri Jul 15 2022 15:12:08 by
1.7.2