PageRenderTime 66ms CodeModel.GetById 30ms RepoModel.GetById 0ms app.codeStats 0ms

/_3rdParty/MathLib/Matrix.h

https://bitbucket.org/byzhang/mldemos
C Header | 1634 lines | 1338 code | 115 blank | 181 comment | 108 complexity | 7f28a7fc79938007d9ca2e45eee47331 MD5 | raw file
Possible License(s): BSD-3-Clause
  1. /*
  2. * Copyright (C) 2010 Learning Algorithms and Systems Laboratory, EPFL, Switzerland
  3. * Author: Eric Sauser
  4. * email: eric.sauser@a3.epf.ch
  5. * website: lasa.epfl.ch
  6. *
  7. * Permission is granted to copy, distribute, and/or modify this program
  8. * under the terms of the GNU General Public License, version 2 or any
  9. * later version published by the Free Software Foundation.
  10. *
  11. * This program is distributed in the hope that it will be useful, but
  12. * WITHOUT ANY WARRANTY; without even the implied warranty of
  13. * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
  14. * Public License for more details
  15. */
  16. #ifndef MATRIX_H
  17. #define MATRIX_H
  18. #include <assert.h>
  19. #include <vector>
  20. using namespace std;
  21. #include "MathLibCommon.h"
  22. #include "Macros.h"
  23. #include "Vector.h"
  24. #ifdef USE_MATHLIB_NAMESPACE
  25. namespace MathLib {
  26. #endif
  27. #ifdef USE_T_EXTENSIONS
  28. template<unsigned int ROW> class TMatrix;
  29. #endif
  30. /**
  31. * \class Matrix
  32. *
  33. * \ingroup MathLib
  34. *
  35. * \brief The basic matrix class
  36. *
  37. * This matrix class can be used for doing various matrix manipulation.
  38. * This shFould be combined with the Vector class for doing almost anything
  39. * you ever dreamt of. Please understand almost as almost...
  40. */
  41. class Matrix
  42. {
  43. friend class Vector;
  44. #ifdef USE_T_EXTENSIONS
  45. template<unsigned int ROW> friend class TMatrix;
  46. #endif
  47. protected:
  48. static int bInverseOk; ///< Tell if last inverse operation was sucessfull
  49. unsigned int row; ///< Number of rows
  50. unsigned int column; ///< Number of columns
  51. REALTYPE *_; ///< The data array
  52. public:
  53. /// Empty constructor
  54. inline Matrix() {
  55. row = 0;
  56. column = 0;
  57. _ = NULL;
  58. }
  59. /// Destructor
  60. inline virtual ~Matrix(){
  61. Release();
  62. }
  63. /// Copy constructor
  64. inline Matrix(const Matrix &matrix)
  65. {
  66. row = 0;
  67. column = 0;
  68. _ = NULL;
  69. Set(matrix);
  70. }
  71. /**
  72. * \brief Create a sized matrix
  73. * \param rowSize The row size
  74. * \param colSize The column size
  75. * \param clear Tell if the matrix sould be set to 0
  76. */
  77. inline Matrix(unsigned int rowSize, unsigned int colSize, bool clear = true)
  78. {
  79. row = 0;
  80. column = 0;
  81. _ = NULL;
  82. Resize(rowSize,colSize,false);
  83. if(clear) Zero();
  84. }
  85. /**
  86. * \brief Create a matrix by copying a memory area
  87. * \param _ The array of REALTYPE to be copied
  88. * \param rowSize The row size
  89. * \param colSize The column size
  90. */
  91. inline Matrix(const REALTYPE *array, unsigned int rowSize, unsigned int colSize)
  92. {
  93. row = 0;
  94. column = 0;
  95. _ = NULL;
  96. Set(array,rowSize,colSize);
  97. }
  98. #ifdef USE_T_EXTENSIONS
  99. /// Copy Contructor of a template Matrix (see TMatrix)
  100. template<unsigned int ROW> inline Matrix(const TMatrix<ROW> &matrix)
  101. {
  102. row = 0;
  103. column = 0;
  104. _ = NULL;
  105. Set(matrix._,ROW,ROW);
  106. }
  107. #endif
  108. /**
  109. * \brief Set a matrix by copying a memory area
  110. * \param _ The array of REALTYPE to be copied
  111. * \param rowSize The row size
  112. * \param colSize The column size
  113. */
  114. inline Matrix& Set(const REALTYPE *array, unsigned int rowSize, unsigned int colSize){
  115. Resize(rowSize,colSize,false);
  116. if((row)&&(column)&&(array))
  117. memcpy(_,array,rowSize*colSize*sizeof(REALTYPE));
  118. return *this;
  119. }
  120. /**
  121. * \brief Set a matrix by copying a memory area
  122. * \param _ The array of REALTYPE to be copied
  123. * \param rowSize The row size
  124. * \param colSize The column size
  125. */
  126. inline Matrix& Set(const Matrix & matrix){
  127. Resize(matrix.row,matrix.column,false);
  128. if((row)&&(column))
  129. memcpy(_,matrix._,row*column*sizeof(REALTYPE));
  130. return *this;
  131. }
  132. /// Get the number of rows
  133. inline unsigned int RowSize() const{
  134. return row;
  135. }
  136. /// Get the number of columns
  137. inline unsigned int ColumnSize() const{
  138. return column;
  139. }
  140. /// Get the data array
  141. inline REALTYPE *Array() const{
  142. return _;
  143. }
  144. /// Clear the matrix
  145. inline Matrix& Zero()
  146. {
  147. if((row)&&(column))
  148. memset(_,0,row*column*sizeof(REALTYPE));
  149. return *this;
  150. }
  151. /// Identity matrix
  152. inline Matrix& Identity()
  153. {
  154. unsigned int k = (row>column?column:row);
  155. Zero();
  156. REALTYPE *ptr = _;
  157. while(k--){
  158. *ptr = R_ONE;
  159. ptr += column+1;
  160. }
  161. return *this;
  162. }
  163. /// Set all elements to 1
  164. inline Matrix& One()
  165. {
  166. REALTYPE *src = _;
  167. unsigned int len = row*column;
  168. while(len--)
  169. *(src++) = R_ONE;
  170. return *this;
  171. }
  172. /// Access to a reference of any element of the matrix
  173. inline REALTYPE& operator() (const unsigned int row, const unsigned int col){
  174. if((row<this->row)&&(col<this->column))
  175. return _[row*column+col];
  176. return Vector::undef;
  177. }
  178. /// Access to a reference of any element of the matrix
  179. inline REALTYPE& Ref(const unsigned int row, const unsigned int col){
  180. if((row<this->row)&&(col<this->column))
  181. return _[row*column+col];
  182. return Vector::undef;
  183. }
  184. /// Access to an element of the matrix
  185. inline REALTYPE At(const unsigned int row, const unsigned int col) const{
  186. if((row<this->row)&&(col<this->column))
  187. return _[row*column+col];
  188. return Vector::undef;
  189. }
  190. /// Access to a reference of any element of the matrix without checks
  191. inline REALTYPE& RefNoCheck(const unsigned int row, const unsigned int col){
  192. return _[row*column+col];
  193. }
  194. /// Access to an element of the matrix without checks
  195. inline REALTYPE AtNoCheck(const unsigned int row, const unsigned int col) const{
  196. return _[row*column+col];
  197. }
  198. /// Get a copy of a row into a vector
  199. inline Vector GetRow(const unsigned int r) const
  200. {
  201. Vector result;
  202. return GetRow(r,result);
  203. }
  204. /// Get a copy of a row into a vector
  205. inline Vector& GetRow(const unsigned int r, Vector& result) const
  206. {
  207. result.Resize(column,false);
  208. if(r<this->row){
  209. REALTYPE *src = _ + r*column;
  210. REALTYPE *dst = result._;
  211. unsigned int len = column;
  212. while(len--)
  213. *(dst++) = *(src++);
  214. }else{
  215. result.Zero();
  216. }
  217. return result;
  218. }
  219. /// Get a copy of a column into a vector
  220. inline Vector GetColumn(const unsigned int col) const
  221. {
  222. Vector result;
  223. return GetColumn(col,result);
  224. }
  225. /// Get a copy of a column into a vector
  226. inline Vector& GetColumn(const unsigned int col, Vector& result) const
  227. {
  228. result.Resize(row,false);
  229. if(col<column){
  230. REALTYPE *src = _ + col;
  231. REALTYPE *dst = result._;
  232. unsigned int len = row;
  233. while(len--){
  234. *(dst++) = *(src);
  235. src += column;
  236. }
  237. }else{
  238. result.Zero();
  239. }
  240. return result;
  241. }
  242. inline Matrix GetRows(unsigned int rowStart, unsigned int rowEnd) const{
  243. Matrix result;
  244. return GetRows(rowStart,rowEnd,result);
  245. }
  246. inline Matrix& GetRows(unsigned int rowStart, unsigned int rowEnd, Matrix& result) const{
  247. if(rowStart>rowEnd){
  248. return result.Resize(0,0,false);
  249. }
  250. unsigned int rowLength = rowEnd-rowStart+1;
  251. result.Resize(rowLength,column,false);
  252. if(rowStart>=row){
  253. return result.Zero();
  254. }
  255. REALTYPE *src = _ + rowStart*column;
  256. REALTYPE *dst = result._;
  257. unsigned int len;
  258. if(rowEnd<row) len = rowLength*column;
  259. else len = (row-rowStart)*column;
  260. while(len--)
  261. *(dst++) = *(src++);
  262. if(rowEnd>=row){
  263. len = (rowEnd-row+1)*column;
  264. while(len--)
  265. *(dst++) = R_ZERO;
  266. }
  267. return result;
  268. }
  269. inline Matrix GetColumns(unsigned int colStart, unsigned int colEnd) const{
  270. Matrix result;
  271. return GetColumns(colStart,colEnd,result);
  272. }
  273. inline Matrix& GetColumns(unsigned int colStart, unsigned int colEnd, Matrix& result) const{
  274. if(colStart>colEnd){
  275. return result.Resize(0,0,false);
  276. }
  277. unsigned int colLength = colEnd-colStart+1;
  278. result.Resize(row,colLength,false);
  279. if(colStart>=column){
  280. return result.Zero();
  281. }
  282. REALTYPE *src = _ + colStart;
  283. REALTYPE *dst = result._;
  284. unsigned int dstOffset = 0;
  285. if(colEnd>=column){
  286. result.Zero();
  287. dstOffset = (colEnd-column+1);
  288. colEnd = column-1;
  289. colLength = colEnd-colStart+1;
  290. }
  291. unsigned int rowLen = row;
  292. while(rowLen--){
  293. unsigned int colLen = colLength;
  294. REALTYPE *cSrc = src;
  295. while(colLen--)
  296. *(dst++) = *(cSrc++);
  297. dst += dstOffset;
  298. src += column;
  299. }
  300. return result;
  301. }
  302. inline Matrix GetMatrix(unsigned int rowStart, unsigned int rowEnd, unsigned int colStart, unsigned int colEnd) const{
  303. Matrix result;
  304. return GetMatrix(rowStart,rowEnd,colStart,colEnd,result);
  305. }
  306. inline Matrix& GetMatrix(unsigned int rowStart, unsigned int rowEnd, unsigned int colStart, unsigned int colEnd, Matrix &result) const{
  307. if((rowStart>rowEnd)||(colStart>colEnd)){
  308. return result.Resize(0,0,false);
  309. }
  310. unsigned int rowLength = rowEnd-rowStart+1;
  311. unsigned int colLength = colEnd-colStart+1;
  312. result.Resize(rowLength,colLength,false);
  313. if((rowStart>=row)||(colStart>=column)){
  314. return result.Zero();
  315. }
  316. unsigned int dstOffset = 0;
  317. if((rowEnd>=row)||(colEnd>=column)){
  318. result.Zero();
  319. if(rowEnd>=row){
  320. rowEnd = row-1;
  321. rowLength = rowEnd-rowStart+1;
  322. }
  323. if(colEnd>=column){
  324. dstOffset = (colEnd-column+1);
  325. colEnd = column-1;
  326. colLength = colEnd-colStart+1;
  327. }
  328. }
  329. REALTYPE *src = _ + colStart + rowStart*column;
  330. REALTYPE *dst = result._;
  331. unsigned int rowLen = rowLength;
  332. while(rowLen--){
  333. REALTYPE *cSrc = src;
  334. unsigned int colLen = colLength;
  335. while(colLen--){
  336. *(dst++) = *(cSrc++);
  337. }
  338. src += column;
  339. dst += dstOffset;
  340. }
  341. return result;
  342. }
  343. /**
  344. * \brief Get a matrix spanning several rows of the matrix
  345. * \param row The starting row
  346. * \param rowSize The number of rows
  347. * \return The resulting matrix
  348. */
  349. inline Matrix GetRowSpace(const unsigned int row, const unsigned int len) const
  350. {
  351. Matrix result;
  352. return GetRowSpace(row,len,result);
  353. }
  354. /**
  355. * \brief Get a matrix spanning several rows of the matrix
  356. * \param row The starting row
  357. * \param len The number of rows
  358. * \param result The target matrix
  359. * \return The resulting matrix
  360. */
  361. inline Matrix& GetRowSpace(const unsigned int row, const unsigned int len, Matrix &result) const
  362. {
  363. if(len>0){
  364. return GetRows(row,row+len-1,result);
  365. }
  366. return result.Resize(0,0,false);
  367. }
  368. /**
  369. * \brief Get a matrix spanning several columns of the matrix
  370. * \param row The starting column
  371. * \param len The number of columns
  372. * \return The resulting matrix
  373. */
  374. inline Matrix GetColumnSpace(const unsigned int col, const unsigned int len) const
  375. {
  376. Matrix result;
  377. return GetColumnSpace(col,len,result);
  378. }
  379. /**
  380. * \brief Get a matrix spanning several columns of the matrix
  381. * \param row The starting column
  382. * \param len The number of columns
  383. * \param result The target matrix
  384. * \return The resulting matrix
  385. */
  386. inline Matrix& GetColumnSpace(const unsigned int col, const unsigned int len, Matrix &result) const
  387. {
  388. if(len>0){
  389. return GetColumns(col,col+len-1,result);
  390. }
  391. return result.Resize(0,0,false);
  392. }
  393. inline Matrix GetMatrixSpace(const unsigned int row, const unsigned int rowLen,
  394. const unsigned int col, const unsigned int colLen) const
  395. {
  396. Matrix result;
  397. return GetMatrixSpace(row,rowLen,col,colLen,result);
  398. }
  399. inline Matrix GetMatrixSpace(const unsigned int row, const unsigned int rowLen,
  400. const unsigned int col, const unsigned int colLen,
  401. Matrix& result) const
  402. {
  403. if((rowLen>0)&&(colLen>0)){
  404. return GetMatrix(row,row+rowLen-1,col,col+colLen-1,result);
  405. }
  406. return result.Resize(0,0,false);
  407. }
  408. /**
  409. * \brief Assign a value to a matrix row
  410. * \param value The value
  411. * \param row The row
  412. */
  413. inline Matrix& SetRow(const REALTYPE value, const unsigned int row)
  414. {
  415. if(row<this->row){
  416. REALTYPE *src = _+ row*column;
  417. unsigned int len = column;
  418. while(len--){
  419. *(src++) = value;
  420. }
  421. }
  422. return *this;
  423. }
  424. /**
  425. * \brief Assign a vector to a matrix row
  426. * \param vector The input vector
  427. * \param row The row
  428. */
  429. inline Matrix& SetRow(const Vector &vector, const unsigned int row)
  430. {
  431. if(row<this->row){
  432. REALTYPE *src = vector._;
  433. REALTYPE *dst = _+ row*column;
  434. unsigned int len = (column<=vector.row?column:vector.row);
  435. while(len--){
  436. *(dst++) = *(src++);
  437. }
  438. }
  439. return *this;
  440. }
  441. /**
  442. * \brief Assign a value to a matrix column
  443. * \param value The value
  444. * \param row The row
  445. */
  446. inline Matrix& SetColumn(const REALTYPE value, const unsigned int col)
  447. {
  448. if(col<this->column){
  449. REALTYPE *src = _ + col;
  450. unsigned int len = row;
  451. while(len--){
  452. *src = value;
  453. src += column;
  454. }
  455. }
  456. return *this;
  457. }
  458. /**
  459. * \brief Assign a vector to a matrix column
  460. * \param vector The input vector
  461. * \param col The column
  462. */
  463. inline Matrix& SetColumn(const Vector &vector, const unsigned int col)
  464. {
  465. if(col<this->column){
  466. REALTYPE *src = vector._;
  467. REALTYPE *dst = _+ col;
  468. unsigned int len = (row<=vector.row?row:vector.row);
  469. while(len--){
  470. *dst = *(src++);
  471. dst += column;
  472. }
  473. }
  474. return *this;
  475. }
  476. /**
  477. * \brief Assign a matrix to the current matrix rows
  478. * \param vector The input matrix
  479. * \param row The starting row
  480. */
  481. inline Matrix& SetRowSpace(const Matrix &matrix, const unsigned int row)
  482. {
  483. if(row<this->row){
  484. const unsigned int ki = (column<=matrix.column?column:matrix.column);
  485. const unsigned int kj = (row+matrix.row<=this->row?row+matrix.row:this->row);
  486. for (unsigned int j = row; j < kj; j++)
  487. for (unsigned int i = 0; i < ki; i++)
  488. _[j*column+i] = matrix._[(j-row)*matrix.column+i];
  489. }
  490. return *this;
  491. }
  492. /**
  493. * \brief Assign a matrix to the current matrix columns
  494. * \param vector The input matrix
  495. * \param col The starting column
  496. */
  497. inline Matrix& SetColumnSpace(const Matrix &matrix, const unsigned int col)
  498. {
  499. if(col<this->column){
  500. const unsigned int kj = (row<=matrix.row?row:matrix.row);
  501. const unsigned int ki = (col+matrix.column<=this->column?col+matrix.column:this->column);
  502. for (unsigned int j = 0; j < kj; j++)
  503. for (unsigned int i = col; i < ki; i++)
  504. _[j*column+i] = matrix._[j*matrix.column+(i-col)];
  505. }
  506. return *this;
  507. }
  508. inline Matrix GetRowSpace(const Vector& ids, Matrix &result) const
  509. {
  510. int s = ids.Size();
  511. IndicesVector id;
  512. for(int i=0;i<s;i++) id.push_back(int(ROUND(ids.At(i))));
  513. return GetRowSpace(id,result);
  514. }
  515. /**
  516. * \brief Get a matrix spanning several rows of the matrix
  517. * \param ids The indices of the desired rows
  518. * \return The resulting matrix
  519. */
  520. inline Matrix GetRowSpace(const IndicesVector& ids) const
  521. {
  522. Matrix result(ids.size(),column);
  523. return GetRowSpace(ids,result);
  524. }
  525. inline Matrix& GetRowSpace(const IndicesVector& ids, Matrix &result) const
  526. {
  527. const unsigned int k=ids.size();
  528. result.Resize(k,column,false);
  529. for(unsigned int i=0;i<k;i++){
  530. const unsigned int g = ids[i];
  531. const unsigned int offset = i*column;
  532. if(g<row){
  533. for(unsigned int j=0;j<column;j++)
  534. result._[offset+j] = _[g*column+j];
  535. }else{
  536. for(unsigned int j=0;j<column;j++)
  537. result._[offset+j] = R_ZERO;
  538. }
  539. }
  540. return result;
  541. }
  542. inline Matrix GetColumnSpace(const Vector& ids, Matrix &result) const
  543. {
  544. int s = ids.Size();
  545. IndicesVector id;
  546. for(int i=0;i<s;i++) id.push_back(int(ROUND(ids.At(i))));
  547. return GetColumnSpace(id,result);
  548. }
  549. /**
  550. * \brief Get a matrix spanning several columns of the matrix
  551. * \param ids The indices of the desired columns
  552. * \return The resulting matrix
  553. */
  554. inline Matrix GetColumnSpace(const IndicesVector& ids) const
  555. {
  556. Matrix result(row,ids.size());
  557. return GetColumnSpace(ids,result);
  558. }
  559. inline Matrix& GetColumnSpace(const IndicesVector& ids, Matrix &result) const
  560. {
  561. const unsigned int k=ids.size();
  562. result.Resize(row,k);
  563. for(unsigned int i=0;i<k;i++){
  564. const unsigned int g = ids[i];
  565. if(g<column){
  566. for(unsigned int j=0;j<row;j++)
  567. result._[j*k+i] = _[j*column+g];
  568. }else{
  569. for(unsigned int j=0;j<row;j++)
  570. result._[j*k+i] = R_ZERO;
  571. }
  572. }
  573. return result;
  574. }
  575. inline Matrix GetMatrixSpace(const Vector& rowIds, const Vector& colIds, Matrix &result) const
  576. {
  577. int r = rowIds.Size();
  578. int c = colIds.Size();
  579. IndicesVector idr,idc;
  580. for(int i=0;i<r;i++) idr.push_back(int(ROUND(rowIds.At(i))));
  581. for(int i=0;i<c;i++) idc.push_back(int(ROUND(colIds.At(i))));
  582. return GetMatrixSpace(idr,idc,result);
  583. }
  584. /**
  585. * \brief Get a matrix spanning several rows and columns of the matrix
  586. * \param rowIds The indices of the desired rows
  587. * \param colIds The indices of the desired columns
  588. * \return The resulting matrix
  589. */
  590. inline Matrix GetMatrixSpace(const IndicesVector& rowIds,const IndicesVector& colIds) const
  591. {
  592. Matrix result(rowIds.size(),colIds.size());
  593. return GetMatrixSpace(rowIds,colIds,result);
  594. }
  595. inline Matrix& GetMatrixSpace(const IndicesVector& rowIds,const IndicesVector& colIds, Matrix &result) const
  596. {
  597. const unsigned int k1=rowIds.size();
  598. const unsigned int k2=colIds.size();
  599. result.Resize(k1,k2);
  600. for(unsigned int i=0;i<k1;i++){
  601. const unsigned int g1 = rowIds[i];
  602. if(g1<row){
  603. for(unsigned int j=0;j<k2;j++){
  604. const unsigned int g2 = colIds[j];
  605. if(g2<column){
  606. result._[i*k2+j] = _[g1*column+g2];
  607. }else{
  608. result._[i*k2+j] = R_ZERO;
  609. }
  610. }
  611. }else{
  612. for(unsigned int j=0;j<k2;j++)
  613. result._[i*k2+j] = R_ZERO;
  614. }
  615. }
  616. return result;
  617. }
  618. /**
  619. * \brief Set a specified column set of the matrix with columns of the source matrix
  620. * \param ids The indices of the desired columns
  621. * \param source The source matrix
  622. * \return The resulting matrix
  623. */
  624. inline Matrix& SetColumnSpace(const IndicesVector& ids, const Matrix &source)
  625. {
  626. const unsigned int k = MIN(ids.size(),source.column);
  627. const unsigned int r = MIN(row,source.row);
  628. for(unsigned int i=0;i<k;i++){
  629. const unsigned int g = ids[i];
  630. if(g<column){
  631. for(unsigned int j=0;j<r;j++)
  632. _[j*column+g] = source._[j*source.column+i];
  633. }
  634. }
  635. return *this;
  636. }
  637. inline Matrix& InsertSubRow(unsigned int startRow, unsigned int startColumn,
  638. const Matrix& matrix,
  639. unsigned int matrixRow,
  640. unsigned int matrixStartColumn, unsigned int matrixColumnLength){
  641. // Check submatrix boundaries
  642. if((matrixRow >= matrix.row)||
  643. (matrixStartColumn >= matrix.column)) return *this;
  644. // Check matrix boundaries
  645. if((startRow >= row)||
  646. (startColumn >= column)) return *this;
  647. if(matrixStartColumn+matrixColumnLength > matrix.column) matrixColumnLength = matrix.column-matrixStartColumn;
  648. if(startColumn+matrixColumnLength > column) matrixColumnLength = column-startColumn;
  649. unsigned int rowOffset = startRow*column;
  650. unsigned int matrixRowOffset = matrixRow*matrix.column;
  651. unsigned int colOffset = startColumn;
  652. unsigned int matrixColOffset = matrixStartColumn;
  653. for(unsigned int j=0;j<matrixColumnLength;j++){
  654. _[rowOffset+colOffset] = matrix._[matrixRowOffset+matrixColOffset];
  655. colOffset++;
  656. matrixColOffset++;
  657. }
  658. return *this;
  659. }
  660. inline Matrix& InsertSubColumn(unsigned int startRow, unsigned int startColumn,
  661. const Matrix& matrix,
  662. unsigned int matrixStartRow, unsigned int matrixRowLength,
  663. unsigned int matrixColumn){
  664. if((matrixStartRow >= matrix.row)||
  665. (matrixColumn >= matrix.column)) return *this;
  666. if((startRow >= row)||
  667. (startColumn >= column)) return *this;
  668. if(matrixStartRow+ matrixRowLength > matrix.row) matrixRowLength = matrix.row -matrixStartRow;
  669. if(startRow+ matrixRowLength > row) matrixRowLength = row -startRow;
  670. unsigned int rowOffset = startRow*column;
  671. unsigned int matrixRowOffset = matrixStartRow*matrix.column;
  672. for(unsigned int i=0;i<matrixRowLength;i++){
  673. _[rowOffset+startColumn] = matrix._[matrixRowOffset+matrixColumn];
  674. rowOffset +=column;
  675. matrixRowOffset +=matrix.column;
  676. }
  677. return *this;
  678. }
  679. inline Matrix& InsertSubMatrix(unsigned int startRow, unsigned int startColumn,
  680. const Matrix& matrix,
  681. unsigned int matrixStartRow, unsigned int matrixRowLength,
  682. unsigned int matrixStartColumn, unsigned int matrixColumnLength){
  683. // Check submatrix boundaries
  684. if((matrixStartRow >= matrix.row)||
  685. (matrixStartColumn >= matrix.column)) return *this;
  686. // Check matrix boundaries
  687. if((startRow >= row)||
  688. (startColumn >= column)) return *this;
  689. if(matrixStartRow+ matrixRowLength > matrix.row) matrixRowLength = matrix.row -matrixStartRow;
  690. if(matrixStartColumn+matrixColumnLength > matrix.column) matrixColumnLength = matrix.column-matrixStartColumn;
  691. if(startRow+ matrixRowLength > row) matrixRowLength = row -startRow;
  692. if(startColumn+matrixColumnLength > column) matrixColumnLength = column-startColumn;
  693. unsigned int rowOffset = startRow*column;
  694. unsigned int matrixRowOffset = matrixStartRow*matrix.column;
  695. for(unsigned int i=0;i<matrixRowLength;i++){
  696. unsigned int colOffset = startColumn;
  697. unsigned int matrixColOffset = matrixStartColumn;
  698. for(unsigned int j=0;j<matrixColumnLength;j++){
  699. _[rowOffset+colOffset] = matrix._[matrixRowOffset+matrixColOffset];
  700. colOffset++;
  701. matrixColOffset++;
  702. }
  703. rowOffset +=column;
  704. matrixRowOffset +=matrix.column;
  705. }
  706. return *this;
  707. }
  708. inline Matrix operator - () const
  709. {
  710. Matrix result(row,column,false);
  711. REALTYPE *src = _;
  712. REALTYPE *dst = result._;
  713. unsigned int len = row*column;
  714. while(len--){
  715. *(dst++) = -(*(src++));
  716. }
  717. return result;
  718. }
  719. inline Matrix& SMinus()
  720. {
  721. REALTYPE *src = _;
  722. unsigned int len = row*column;
  723. while(len--){
  724. *(src) = -(*(src));
  725. src++;
  726. }
  727. return *this;
  728. }
  729. inline virtual Matrix& operator = (const Matrix &matrix)
  730. {
  731. return Set(matrix);
  732. }
  733. inline Matrix& operator += (const Matrix &matrix)
  734. {
  735. const unsigned int kj = (row<=matrix.row?row:matrix.row);
  736. REALTYPE *dst = _;
  737. REALTYPE *src = matrix._;
  738. if(column==matrix.column){
  739. unsigned int len = column*kj;
  740. while(len--){
  741. *(dst++) += *(src++);
  742. }
  743. }else if(column < matrix.column){
  744. unsigned int colLen = column;
  745. unsigned int colDif = matrix.column-column;
  746. while(colLen--){
  747. unsigned int rowLen = kj;
  748. while(rowLen--){
  749. *(dst++) += *(src++);
  750. }
  751. src += colDif;
  752. }
  753. }else{
  754. unsigned int colLen = matrix.column;
  755. unsigned int colDif = column - matrix.column;
  756. while(colLen--){
  757. unsigned int rowLen = kj;
  758. while(rowLen--){
  759. *(dst++) += *(src++);
  760. }
  761. dst += colDif;
  762. }
  763. }
  764. return *this;
  765. }
  766. inline Matrix& operator -= (const Matrix &matrix)
  767. {
  768. const unsigned int kj = (row<=matrix.row?row:matrix.row);
  769. REALTYPE *dst = _;
  770. REALTYPE *src = matrix._;
  771. if(column==matrix.column){
  772. unsigned int len = column*kj;
  773. while(len--){
  774. *(dst++) -= *(src++);
  775. }
  776. }else if(column < matrix.column){
  777. unsigned int colLen = column;
  778. unsigned int colDif = matrix.column-column;
  779. while(colLen--){
  780. unsigned int rowLen = kj;
  781. while(rowLen--){
  782. *(dst++) -= *(src++);
  783. }
  784. src += colDif;
  785. }
  786. }else{
  787. unsigned int colLen = matrix.column;
  788. unsigned int colDif = column - matrix.column;
  789. while(colLen--){
  790. unsigned int rowLen = kj;
  791. while(rowLen--){
  792. *(dst++) -= *(src++);
  793. }
  794. dst += colDif;
  795. }
  796. }
  797. return *this;
  798. }
  799. /// matrix multiplication
  800. inline Matrix& operator *= (const Matrix &matrix){
  801. Matrix tmp;
  802. Mult(matrix,tmp);
  803. return Swap(tmp);
  804. }
  805. /// Element-wise multiplication
  806. inline Matrix& operator ^= (const Matrix &matrix)
  807. {
  808. const unsigned int kj = (row<=matrix.row?row:matrix.row);
  809. REALTYPE *dst = _;
  810. REALTYPE *src = matrix._;
  811. if(column==matrix.column){
  812. unsigned int len = column*kj;
  813. while(len--){
  814. *(dst++) *= *(src++);
  815. }
  816. }else if(column < matrix.column){
  817. unsigned int colLen = column;
  818. unsigned int colDif = matrix.column-column;
  819. while(colLen--){
  820. unsigned int rowLen = kj;
  821. while(rowLen--){
  822. *(dst++) *= *(src++);
  823. }
  824. src += colDif;
  825. }
  826. }else{
  827. unsigned int colLen = matrix.column;
  828. unsigned int colDif = column - matrix.column;
  829. while(colLen--){
  830. unsigned int rowLen = kj;
  831. while(rowLen--){
  832. *(dst++) *= *(src++);
  833. }
  834. dst += colDif;
  835. }
  836. }
  837. return *this;
  838. }
  839. /* Row wise multiplication by a vector */
  840. inline Matrix& operator ^= (Vector vec)
  841. {
  842. for (unsigned int j = 0; j < row; j++)
  843. {
  844. REALTYPE scalar = vec(j);
  845. for (unsigned int i = 0; i < column; i++)
  846. _[j*column+i] *= scalar;
  847. }
  848. return *this;
  849. }
  850. /* Row wise division by a vector */
  851. inline Matrix& operator /= (Vector vec)
  852. {
  853. for (unsigned int j = 0; j < row; j++)
  854. {
  855. REALTYPE scalar = R_ONE/vec(j);
  856. for (unsigned int i = 0; i < column; i++)
  857. _[j*column+i] *= scalar;
  858. }
  859. return *this;
  860. }
  861. /// Element-wise division
  862. inline Matrix& operator /= (const Matrix &matrix)
  863. {
  864. const unsigned int kj = (row<=matrix.row?row:matrix.row);
  865. REALTYPE *dst = _;
  866. REALTYPE *src = matrix._;
  867. if(column==matrix.column){
  868. unsigned int len = column*kj;
  869. while(len--){
  870. *(dst++) /= *(src++);
  871. }
  872. }else if(column < matrix.column){
  873. unsigned int colLen = column;
  874. unsigned int colDif = matrix.column-column;
  875. while(colLen--){
  876. unsigned int rowLen = kj;
  877. while(rowLen--){
  878. *(dst++) /= *(src++);
  879. }
  880. src += colDif;
  881. }
  882. }else{
  883. unsigned int colLen = matrix.column;
  884. unsigned int colDif = column - matrix.column;
  885. while(colLen--){
  886. unsigned int rowLen = kj;
  887. while(rowLen--){
  888. *(dst++) /= *(src++);
  889. }
  890. dst += colDif;
  891. }
  892. }
  893. return *this;
  894. }
  895. inline Matrix& operator += (REALTYPE scalar)
  896. {
  897. REALTYPE *src = _;
  898. unsigned int len = row*column;
  899. while(len--){
  900. *(src++) += scalar;
  901. }
  902. return *this;
  903. }
  904. inline Matrix& operator -= (REALTYPE scalar)
  905. {
  906. REALTYPE *src = _;
  907. unsigned int len = row*column;
  908. while(len--){
  909. *(src++) -= scalar;
  910. }
  911. return *this;
  912. }
  913. inline Matrix& operator *= (REALTYPE scalar)
  914. {
  915. REALTYPE *src = _;
  916. unsigned int len = row*column;
  917. while(len--){
  918. *(src++) *= scalar;
  919. }
  920. return *this;
  921. }
  922. inline Matrix& operator /= (REALTYPE scalar)
  923. {
  924. scalar = R_ONE/scalar;
  925. return *(this) *= scalar;
  926. }
  927. inline Matrix operator + (const Matrix &matrix) const
  928. {
  929. Matrix result(row,column,false);
  930. return Add(matrix,result);
  931. }
  932. inline Matrix operator - (const Matrix &matrix) const
  933. {
  934. Matrix result(row,column,false);
  935. return Sub(matrix,result);
  936. }
  937. /// Element-wise multiplication
  938. inline Matrix operator ^ (const Matrix &matrix) const
  939. {
  940. Matrix result(row,column,false);
  941. return PMult(matrix,result);
  942. }
  943. inline Matrix operator / (const Matrix &matrix) const
  944. {
  945. Matrix result(row,column,false);
  946. return PDiv(matrix,result);
  947. }
  948. inline Matrix& ScaleAddTo(REALTYPE scale, Matrix &result) const
  949. {
  950. const unsigned int kj = (row<=result.row?row:result.row);
  951. REALTYPE *dst = result._;
  952. REALTYPE *src = _;
  953. if(column==result.column){
  954. unsigned int len = column*kj;
  955. while(len--){
  956. *(dst++) += *(src++) * scale;
  957. }
  958. }else if(result.column < column){
  959. unsigned int colLen = result.column;
  960. unsigned int colDif = column-result.column;
  961. while(colLen--){
  962. unsigned int rowLen = kj;
  963. while(rowLen--){
  964. *(dst++) += *(src++) * scale;
  965. }
  966. src += colDif;
  967. }
  968. }else{
  969. unsigned int colLen = column;
  970. unsigned int colDif = result.column - column;
  971. while(colLen--){
  972. unsigned int rowLen = kj;
  973. while(rowLen--){
  974. *(dst++) += *(src++) * scale;
  975. }
  976. dst += colDif;
  977. }
  978. }
  979. return result;
  980. }
  981. inline Matrix& Add(const Matrix &matrix, Matrix &result) const
  982. {
  983. result.Resize(row,column,false);
  984. const unsigned int kj = (row<=matrix.row?row:matrix.row);
  985. REALTYPE *dst = result._;
  986. REALTYPE *src0 = _;
  987. REALTYPE *src = matrix._;
  988. if(column==matrix.column){
  989. unsigned int len = column*kj;
  990. while(len--){
  991. *(dst++) = *(src0++) + *(src++);
  992. }
  993. }else if(column < matrix.column){
  994. unsigned int colLen = column;
  995. unsigned int colDif = matrix.column-column;
  996. while(colLen--){
  997. unsigned int rowLen = kj;
  998. while(rowLen--){
  999. *(dst++) = *(src0++) + *(src++);
  1000. }
  1001. src += colDif;
  1002. }
  1003. }else{
  1004. unsigned int colLen = matrix.column;
  1005. unsigned int colDif = column - matrix.column;
  1006. while(colLen--){
  1007. unsigned int rowLen = kj;
  1008. while(rowLen--){
  1009. *(dst++) = *(src0++) + *(src++);
  1010. }
  1011. rowLen = colDif;
  1012. while(rowLen--){
  1013. *(dst++) = *(src0++);
  1014. }
  1015. }
  1016. }
  1017. if(kj!=row){
  1018. unsigned int len = column*(kj-matrix.row);
  1019. while(len--){
  1020. *(dst++) = *(src0++);
  1021. }
  1022. }
  1023. return result;
  1024. }
  1025. inline Matrix& Sub(const Matrix &matrix, Matrix &result) const
  1026. {
  1027. result.Resize(row,column,false);
  1028. const unsigned int kj = (row<=matrix.row?row:matrix.row);
  1029. REALTYPE *dst = result._;
  1030. REALTYPE *src0 = _;
  1031. REALTYPE *src = matrix._;
  1032. if(column==matrix.column){
  1033. unsigned int len = column*kj;
  1034. while(len--){
  1035. *(dst++) = *(src0++) - *(src++);
  1036. }
  1037. }else if(column < matrix.column){
  1038. unsigned int colLen = column;
  1039. unsigned int colDif = matrix.column-column;
  1040. while(colLen--){
  1041. unsigned int rowLen = kj;
  1042. while(rowLen--){
  1043. *(dst++) = *(src0++) - *(src++);
  1044. }
  1045. src += colDif;
  1046. }
  1047. }else{
  1048. unsigned int colLen = matrix.column;
  1049. unsigned int colDif = column - matrix.column;
  1050. while(colLen--){
  1051. unsigned int rowLen = kj;
  1052. while(rowLen--){
  1053. *(dst++) = *(src0++) - *(src++);
  1054. }
  1055. rowLen = colDif;
  1056. while(rowLen--){
  1057. *(dst++) = *(src0++);
  1058. }
  1059. }
  1060. }
  1061. if(kj!=row){
  1062. unsigned int len = column*(kj-matrix.row);
  1063. while(len--){
  1064. *(dst++) = *(src0++);
  1065. }
  1066. }
  1067. return result;
  1068. }
  1069. /// Element-wise multiplication
  1070. inline Matrix& PMult(const Matrix &matrix, Matrix &result) const
  1071. {
  1072. result.Resize(row,column,false);
  1073. const unsigned int kj = (row<=matrix.row?row:matrix.row);
  1074. REALTYPE *dst = result._;
  1075. REALTYPE *src0 = _;
  1076. REALTYPE *src = matrix._;
  1077. if(column==matrix.column){
  1078. unsigned int len = column*kj;
  1079. while(len--){
  1080. *(dst++) = *(src0++) * *(src++);
  1081. }
  1082. }else if(column < matrix.column){
  1083. unsigned int colLen = column;
  1084. unsigned int colDif = matrix.column-column;
  1085. while(colLen--){
  1086. unsigned int rowLen = kj;
  1087. while(rowLen--){
  1088. *(dst++) = *(src0++) * *(src++);
  1089. }
  1090. src += colDif;
  1091. }
  1092. }else{
  1093. unsigned int colLen = matrix.column;
  1094. unsigned int colDif = column - matrix.column;
  1095. while(colLen--){
  1096. unsigned int rowLen = kj;
  1097. while(rowLen--){
  1098. *(dst++) = *(src0++) * *(src++);
  1099. }
  1100. rowLen = colDif;
  1101. while(rowLen--){
  1102. *(dst++) = *(src0++);
  1103. }
  1104. }
  1105. }
  1106. if(kj!=row){
  1107. unsigned int len = column*(kj-matrix.row);
  1108. while(len--){
  1109. *(dst++) = *(src0++);
  1110. }
  1111. }
  1112. return result;
  1113. }
  1114. /// Element-wise division
  1115. inline Matrix& PDiv(const Matrix &matrix, Matrix &result) const
  1116. {
  1117. result.Resize(row,column,false);
  1118. const unsigned int kj = (row<=matrix.row?row:matrix.row);
  1119. REALTYPE *dst = result._;
  1120. REALTYPE *src0 = _;
  1121. REALTYPE *src = matrix._;
  1122. if(column==matrix.column){
  1123. unsigned int len = column*kj;
  1124. while(len--){
  1125. *(dst++) = *(src0++) / *(src++);
  1126. }
  1127. }else if(column < matrix.column){
  1128. unsigned int colLen = column;
  1129. unsigned int colDif = matrix.column-column;
  1130. while(colLen--){
  1131. unsigned int rowLen = kj;
  1132. while(rowLen--){
  1133. *(dst++) = *(src0++) / *(src++);
  1134. }
  1135. src += colDif;
  1136. }
  1137. }else{
  1138. unsigned int colLen = matrix.column;
  1139. unsigned int colDif = column - matrix.column;
  1140. while(colLen--){
  1141. unsigned int rowLen = kj;
  1142. while(rowLen--){
  1143. *(dst++) = *(src0++) / *(src++);
  1144. }
  1145. rowLen = colDif;
  1146. while(rowLen--){
  1147. *(dst++) = *(src0++);
  1148. }
  1149. }
  1150. }
  1151. if(kj!=row){
  1152. unsigned int len = column*(kj-matrix.row);
  1153. while(len--){
  1154. *(dst++) = *(src0++);
  1155. }
  1156. }
  1157. return result;
  1158. }
  1159. inline Matrix operator + (REALTYPE scalar) const
  1160. {
  1161. Matrix result(row,column,false);
  1162. return Add(scalar,result);
  1163. }
  1164. inline Matrix operator - (REALTYPE scalar) const
  1165. {
  1166. Matrix result(row,column,false);
  1167. return Sub(scalar,result);
  1168. }
  1169. inline Matrix operator * (REALTYPE scalar) const
  1170. {
  1171. Matrix result(row,column,false);
  1172. return Mult(scalar,result);
  1173. }
  1174. inline Matrix operator / (REALTYPE scalar) const
  1175. {
  1176. Matrix result(row,column,false);
  1177. return Div(scalar,result);
  1178. }
  1179. inline Matrix& Add(REALTYPE scalar, Matrix& result) const
  1180. {
  1181. result.Resize(row,column,false);
  1182. REALTYPE *src = _;
  1183. REALTYPE *dst = result._;
  1184. unsigned int len = row*column;
  1185. while(len--){
  1186. *(dst++) = *(src++) + scalar;
  1187. }
  1188. return result;
  1189. }
  1190. inline Matrix& Sub(REALTYPE scalar, Matrix& result) const
  1191. {
  1192. return Add(-scalar,result);
  1193. }
  1194. inline Matrix& Mult(REALTYPE scalar, Matrix& result) const
  1195. {
  1196. result.Resize(row,column,false);
  1197. REALTYPE *src = _;
  1198. REALTYPE *dst = result._;
  1199. unsigned int len = row*column;
  1200. while(len--){
  1201. *(dst++) = *(src++) * scalar;
  1202. }
  1203. return result;
  1204. }
  1205. inline Matrix& Div(REALTYPE scalar, Matrix& result) const
  1206. {
  1207. return Mult(R_ONE/scalar,result);
  1208. }
  1209. inline bool operator == (const Matrix& matrix) const
  1210. {
  1211. if((row!=matrix.row)||(column!=matrix.column)) return false;
  1212. REALTYPE *src = _;
  1213. REALTYPE *dst = matrix._;
  1214. unsigned int len = row*column;
  1215. while(len--){
  1216. if( *(dst++) != *(src++) ) return false;
  1217. }
  1218. return true;
  1219. }
  1220. inline bool operator != (const Matrix& matrix) const
  1221. {
  1222. return !(*this == matrix);
  1223. }
  1224. /// Matrix multiplication against a vector
  1225. inline Vector operator * (const Vector &vector) const
  1226. {
  1227. Vector result(row,false);
  1228. return Mult(vector,result);
  1229. }
  1230. /// Matrix multiplication against a vector
  1231. inline Vector Mult(const Vector &vector) const
  1232. {
  1233. Vector result(row,false);
  1234. return Mult(vector,result);
  1235. }
  1236. /// Matrix multiplication against a vector
  1237. inline Vector& Mult(const Vector &vector, Vector &result) const
  1238. {
  1239. result.Resize(row,false);
  1240. const unsigned int ki = (column<=vector.row?column:vector.row);
  1241. unsigned int rowLen = row;
  1242. REALTYPE *dst = result._;
  1243. REALTYPE *srcM = _;
  1244. while(rowLen--){
  1245. REALTYPE *srcV = vector._;
  1246. REALTYPE sum = R_ZERO;
  1247. unsigned int colLen = ki;
  1248. while(colLen--){
  1249. sum += (*(srcM++)) * (*(srcV++));
  1250. }
  1251. srcM += column-ki;
  1252. *(dst++) = sum;
  1253. }
  1254. return result;
  1255. }
  1256. /// Transpose Matrix multiplication against a vector
  1257. inline Vector TransposeMult(const Vector &vector) const
  1258. {
  1259. Vector result(column,false);
  1260. return TransposeMult(vector,result);
  1261. }
  1262. /// Transpose Matrix multiplication against a vector
  1263. inline Vector& TransposeMult(const Vector &vector, Vector &result) const
  1264. {
  1265. result.Resize(column,false);
  1266. result.Zero();
  1267. const unsigned int ki = (row<vector.row?row:vector.row);
  1268. REALTYPE *srcV = vector._;
  1269. REALTYPE *srcM = _;
  1270. unsigned int rowLen = ki;
  1271. while(rowLen--){
  1272. REALTYPE *dst = result._;
  1273. unsigned int colLen = column;
  1274. while(colLen--){
  1275. *(dst++) += (*(srcM++)) * (*(srcV));
  1276. }
  1277. srcV++;
  1278. }
  1279. return result;
  1280. }
  1281. inline Matrix& SAddToRow(const Vector &vector){
  1282. const unsigned int k = (column<=vector.row?column:vector.row);
  1283. REALTYPE *ptr = _;
  1284. unsigned int rowLen = row;
  1285. while(rowLen--){
  1286. REALTYPE *ptr2 = ptr;
  1287. REALTYPE *ptrv = vector._;
  1288. unsigned int colLen = k;
  1289. while(colLen--){
  1290. *(ptr2++) += *(ptrv++);
  1291. }
  1292. ptr += column;
  1293. }
  1294. return *this;
  1295. }
  1296. inline Matrix& SAddToColumn(const Vector &vector){
  1297. REALTYPE *ptr = _;
  1298. REALTYPE *ptrv = vector._;
  1299. unsigned int rowLen = (row<=vector.row?row:vector.row);
  1300. while(rowLen--){
  1301. unsigned int colLen = column;
  1302. while(colLen--)
  1303. *(ptr++) += *(ptrv);
  1304. ptrv++;
  1305. }
  1306. return *this;
  1307. }
  1308. /// Multiplication of each ith row by the ith element of a vector
  1309. inline Matrix MultRow(const Vector &vector) const {
  1310. Matrix result;
  1311. return MultRow(vector,result);
  1312. }
  1313. inline Matrix& MultRow(const Vector &vector, Matrix& result) const {
  1314. result.Set(*this);
  1315. return result.SMultRow(vector);
  1316. }
  1317. /// Self multiplication of each ith row by the ith element of a vector
  1318. inline Matrix& SMultRow(const Vector &vector)
  1319. {
  1320. REALTYPE *cDst = _;
  1321. REALTYPE *cSrc = vector._;
  1322. unsigned int ki = (row<=vector.row?row:vector.row);
  1323. while(ki--){
  1324. REALTYPE v = *(cSrc++);
  1325. unsigned int len = column;
  1326. while(len--) (*cDst++) *= v;
  1327. }
  1328. if(row>vector.row){
  1329. unsigned int len = column * (row-vector.row);
  1330. while(len--) *(cDst++) = R_ZERO;
  1331. }
  1332. return *this;
  1333. }
  1334. /// Multiplication of each ith column by the ith element of a vector
  1335. inline Matrix MultColumn(const Vector &vector) const {
  1336. Matrix result;
  1337. return MultColumn(vector,result);
  1338. }
  1339. inline Matrix& MultColumn(const Vector &vector, Matrix& result) const {
  1340. result.Set(*this);
  1341. return result.SMultColumn(vector);
  1342. }
  1343. /// Self multiplication of each ith row by the ith element of a vector
  1344. inline Matrix& SMultColumn(const Vector &vector)
  1345. {
  1346. REALTYPE *cDst = _;
  1347. if(column<=vector.row){
  1348. unsigned int len = row;
  1349. while(len--){
  1350. REALTYPE *cSrc = vector._;
  1351. unsigned int colLen = column;
  1352. while(colLen--) *(cDst++) *= *(cSrc++);
  1353. }
  1354. }else{
  1355. unsigned int len = row;
  1356. while(len--){
  1357. REALTYPE *cSrc = vector._;
  1358. unsigned int colLen = vector.row;
  1359. while(colLen--) *(cDst++) *= *(cSrc++);
  1360. colLen = row-vector.row;
  1361. while(colLen--) *(cDst++) = R_ZERO;
  1362. }
  1363. }
  1364. return *this;
  1365. }
  1366. /// Matrix multiplication
  1367. inline Matrix operator * (const Matrix &matrix) const {
  1368. Matrix result;
  1369. return Mult(matrix,result);
  1370. }
  1371. /// Matrix multiplication
  1372. inline Matrix Mult(const Matrix &matrix) const {
  1373. Matrix result;
  1374. return Mult(matrix,result);
  1375. }
  1376. /// Matrix multiplication
  1377. inline Matrix& Mult(const Matrix &matrix, Matrix &result) const
  1378. {
  1379. result.Resize(row,matrix.column,false);
  1380. result.Zero();
  1381. const unsigned int rcol = result.column;
  1382. const unsigned int kk = (column<=matrix.row?column:matrix.row);
  1383. REALTYPE *cP1 = _;
  1384. REALTYPE *eP1 = cP1 + row*column;
  1385. REALTYPE *cD = result._;
  1386. while(cP1!=eP1){
  1387. REALTYPE *currP1 = cP1;
  1388. REALTYPE *endP1 = currP1 + kk;
  1389. REALTYPE *currP2 = matrix._;
  1390. while(currP1!=endP1){
  1391. REALTYPE *currPD = cD;
  1392. REALTYPE curr1 = *currP1;
  1393. REALTYPE *endP2 = currP2 + rcol;
  1394. while(currP2!=endP2){
  1395. (*currPD++) += curr1 * (*(currP2++));
  1396. }
  1397. currP1++;
  1398. }
  1399. cD += rcol;
  1400. cP1 += column;
  1401. }
  1402. return result;
  1403. }
  1404. inline Matrix& TransposeMult(const Matrix &matrix) const{
  1405. Matrix result;
  1406. return TransposeMult(matrix,result);
  1407. }
  1408. /// Transpose then Multiply by the matrix :
  1409. // result = this^T * matrix (avoiding a useless call to .Transpose() )
  1410. inline Matrix& TransposeMult(const Matrix &matrix, Matrix &result) const
  1411. {
  1412. //[a0 a1 .. an] [b0 b1 .. bn]
  1413. result.Resize(column,matrix.column,false);
  1414. result.Zero();
  1415. unsigned int kk = (row<=matrix.row?row:matrix.row);
  1416. REALTYPE *cP1 = _;
  1417. REALTYPE *cP2 = matrix._;
  1418. while(kk--){
  1419. REALTYPE *currD = result._;
  1420. REALTYPE *currP1 = cP1;
  1421. unsigned int len1 = column;
  1422. while(len1--){
  1423. REALTYPE *currP2 = cP2;
  1424. unsigned int len2 = matrix.column;
  1425. while(len2--){
  1426. *(currD++) += *(currP1) *(*(currP2++));
  1427. }
  1428. currP1++;
  1429. }
  1430. cP2 += matrix.column;
  1431. cP1 += column;
  1432. }
  1433. return result;
  1434. }
  1435. inline Matrix& MultTranspose2(const Matrix &matrix) const{
  1436. Matrix result;
  1437. return MultTranspose2(matrix,result);
  1438. }
  1439. /// Transpose then Multiply by the matrix :
  1440. // result = this^T * matrix (avoiding a useless call to .Transpose() )
  1441. inline Matrix& MultTranspose2(const Matrix &matrix, Matrix &result) const
  1442. {
  1443. result.Resize(row,matrix.row,false);
  1444. //[a0 a1 .. an] [b0 b1 .. bn]
  1445. if(column==matrix.column){
  1446. REALTYPE *cP1 = _;
  1447. REALTYPE *currD = result._;
  1448. unsigned int len1 = row;
  1449. while(len1--){
  1450. REALTYPE *currP2 = matrix._;
  1451. unsigned int len2 = matrix.row;
  1452. while(len2--){
  1453. REALTYPE *currP1 = cP1;
  1454. REALTYPE sum=0.0;
  1455. unsigned int len = column;
  1456. while(len--)
  1457. sum += *(currP1++) * (*(currP2++));
  1458. *(currD++) = sum;
  1459. }
  1460. cP1 += column;
  1461. }
  1462. }else{
  1463. REALTYPE *cP1 = _;
  1464. REALTYPE *currD = result._;
  1465. unsigned int len1 = row;
  1466. const unsigned int kk = (column<=matrix.column?column:matrix.column);
  1467. unsigned int off2 = matrix.column - kk;
  1468. while(len1--){
  1469. REALTYPE *currP2 = matrix._;
  1470. unsigned int len2 = matrix.row;
  1471. while(len2--){
  1472. REALTYPE *currP1 = cP1;
  1473. REALTYPE sum=0.0;
  1474. unsigned int len = kk;
  1475. while(len--)
  1476. sum += *(currP1++) * (*(currP2++));
  1477. *(currD++) = sum;
  1478. currP2 += off2;
  1479. }
  1480. cP1 += column;
  1481. }
  1482. }
  1483. return result;
  1484. }
  1485. /// Set a diagonal matrix given a vector of diagonal elements
  1486. inline Matrix& Diag(const Vector &vector)
  1487. {
  1488. Resize(vector.row,vector.row,false);
  1489. Zero();
  1490. REALTYPE *src = vector._;
  1491. REALTYPE *dst = _;
  1492. unsigned int len = row;
  1493. while(len--){
  1494. *dst = *(src++);
  1495. dst += column+1;
  1496. }
  1497. return *this;
  1498. }
  1499. /// Set a random matrix with value uniformly distributed between 0 and 1
  1500. inline Matrix& Random(){
  1501. REALTYPE *src = _;
  1502. unsigned int len = row*column;
  1503. while(len--){
  1504. *(src++) = RND(R_ONE);
  1505. }
  1506. return *this;
  1507. }
  1508. /// Return the transpose of a matrix
  1509. inline Matrix Transpose() const
  1510. {
  1511. Matrix result;
  1512. return Transpose(result);
  1513. }
  1514. /// Compute the transpose of a matrix
  1515. inline Matrix& Transpose(Matrix &result) const
  1516. {
  1517. result.Resize(column,row,false);
  1518. REALTYPE *src = _;
  1519. REALTYPE *dst = re