home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 January / usenetsourcesnewsgroupsinfomagicjanuary1994.iso / sources / misc / volume26 / newmat04 / part04 / newmat6.cxx < prev    next >
Encoding:
C/C++ Source or Header  |  1991-11-30  |  6.9 KB  |  268 lines
  1. //$$ newmat6.cxx            Operators, element access, submatrices
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  4.  
  5. #include "include.hxx"
  6.  
  7. #include "newmat.hxx"
  8. #include "newmatrc.hxx"
  9.  
  10.  
  11. //#define REPORT { static ExeCounter ExeCount(__LINE__,6); ExeCount++; }
  12.  
  13. #define REPORT {}
  14.  
  15. /*************************** general utilities *************************/
  16.  
  17. static int tristore(int n)                      // els in triangular matrix
  18. { return (n*(n+1))/2; }
  19.  
  20.  
  21. /****************************** operators *******************************/
  22.  
  23. real& Matrix::operator()(int m, int n)
  24. {
  25.    if (m<=0 || m>nrows || n<=0 || n>ncols) MatrixError("Index out of range");
  26.    return store[(m-1)*ncols+n-1];
  27. }
  28.  
  29. real& SymmetricMatrix::operator()(int m, int n)
  30. {
  31.    if (m<=0 || n<=0 || m>nrows || n>ncols) MatrixError("Index out of range");
  32.    if (m>=n) return store[tristore(m-1)+n-1];
  33.    else return store[tristore(n-1)+m-1];
  34. }
  35.  
  36. real& UpperTriangularMatrix::operator()(int m, int n)
  37. {
  38.    if (m<=0 || n<m || n>ncols) MatrixError("Index out of range");
  39.    return store[(m-1)*ncols+n-1-tristore(m-1)];
  40. }
  41.  
  42. real& LowerTriangularMatrix::operator()(int m, int n)
  43. {
  44.    if (n<=0 || m<n || m>nrows) MatrixError("Index out of range");
  45.    return store[tristore(m-1)+n-1];
  46. }
  47.  
  48. real& DiagonalMatrix::operator()(int m, int n)
  49. {
  50.    if (n<=0 || m!=n || m>nrows || n>ncols) MatrixError("Index out of range");
  51.    return store[n-1];
  52. }
  53.  
  54. real& DiagonalMatrix::operator()(int m)
  55. {
  56.    if (m<=0 || m>nrows) MatrixError("Index out of range");
  57.    return store[m-1];
  58. }
  59.  
  60. real& ColumnVector::operator()(int m)
  61. {
  62.    if (m<=0 || m> nrows) MatrixError("Index out of range");
  63.    return store[m-1];
  64. }
  65.  
  66. real& RowVector::operator()(int n)
  67. {
  68.    if (n<=0 || n> ncols) MatrixError("Index out of range");
  69.    return store[n-1];
  70. }
  71.  
  72. #ifndef __ZTC__
  73.  
  74. real Matrix::operator()(int m, int n) const
  75. {
  76.    if (m<=0 || m>nrows || n<=0 || n>ncols) MatrixError("Index out of range");
  77.    return store[(m-1)*ncols+n-1];
  78. }
  79.  
  80. real SymmetricMatrix::operator()(int m, int n) const
  81. {
  82.    if (m<=0 || n<=0 || m>nrows || n>ncols) MatrixError("Index out of range");
  83.    if (m>=n) return store[tristore(m-1)+n-1];
  84.    else return store[tristore(n-1)+m-1];
  85. }
  86.  
  87. real UpperTriangularMatrix::operator()(int m, int n) const
  88. {
  89.    if (m<=0 || n<m || n>ncols) MatrixError("Index out of range");
  90.    return store[(m-1)*ncols+n-1-tristore(m-1)];
  91. }
  92.  
  93. real LowerTriangularMatrix::operator()(int m, int n) const
  94. {
  95.    if (n<=0 || m<n || m>nrows) MatrixError("Index out of range");
  96.    return store[tristore(m-1)+n-1];
  97. }
  98.  
  99. real DiagonalMatrix::operator()(int m, int n) const
  100. {
  101.    if (n<=0 || m!=n || m>nrows || n>ncols) MatrixError("Index out of range");
  102.    return store[n-1];
  103. }
  104.  
  105. real DiagonalMatrix::operator()(int m) const
  106. {
  107.    if (m<=0 || m>nrows) MatrixError("Index out of range");
  108.    return store[m-1];
  109. }
  110.  
  111. real ColumnVector::operator()(int m) const
  112. {
  113.    if (m<=0 || m> nrows) MatrixError("Index out of range");
  114.    return store[m-1];
  115. }
  116.  
  117. real RowVector::operator()(int n) const
  118. {
  119.    if (n<=0 || n> ncols) MatrixError("Index out of range");
  120.    return store[n-1];
  121. }
  122.  
  123. #endif
  124.  
  125. BaseMatrix::operator real()
  126. {
  127.    REPORT
  128.    GeneralMatrix* gm = Evaluate();
  129.    if (gm->nrows!=1 || gm->ncols!=1)
  130.       MatrixError("Attempt to convert non 1x1 matrix to scalar");
  131.    real x = *(gm->store); gm->tDelete(); return x;
  132. }
  133.  
  134. AddedMatrix BaseMatrix::operator+(BaseMatrix& bm)
  135. { REPORT return AddedMatrix(this, &bm); }
  136.  
  137. MultipliedMatrix BaseMatrix::operator*(BaseMatrix& bm)
  138. { REPORT return MultipliedMatrix(this, &bm); }
  139.  
  140. SolvedMatrix InvertedMatrix::operator*(BaseMatrix& bmx)
  141. { REPORT return SolvedMatrix(bm, &bmx); }
  142.  
  143. SubtractedMatrix BaseMatrix::operator-(BaseMatrix& bm)
  144. { REPORT return SubtractedMatrix(this, &bm); }
  145.  
  146. ShiftedMatrix BaseMatrix::operator+(real f)
  147. { REPORT return ShiftedMatrix(this, f); }
  148.  
  149. ScaledMatrix BaseMatrix::operator*(real f)
  150. { REPORT return ScaledMatrix(this, f); }
  151.  
  152. ScaledMatrix BaseMatrix::operator/(real f)
  153. { REPORT return ScaledMatrix(this, 1.0/f); }
  154.  
  155. ShiftedMatrix BaseMatrix::operator-(real f)
  156. { REPORT return ShiftedMatrix(this, -f); }
  157.  
  158. TransposedMatrix BaseMatrix::t() { REPORT return TransposedMatrix(this); }
  159.  
  160. NegatedMatrix BaseMatrix::operator-() { REPORT return NegatedMatrix(this); }
  161.  
  162. InvertedMatrix BaseMatrix::i() { REPORT return InvertedMatrix(this); }
  163.  
  164. ConstMatrix GeneralMatrix::c() const
  165. {
  166.    if (tag != -1) MatrixError(".c() applied to temporary matrix");
  167.    REPORT return ConstMatrix(this);
  168. }
  169.  
  170. RowedMatrix BaseMatrix::CopyToRow() { REPORT return RowedMatrix(this); }
  171.  
  172. ColedMatrix BaseMatrix::CopyToColumn() { REPORT return ColedMatrix(this); }
  173.  
  174. DiagedMatrix BaseMatrix::CopyToDiagonal() { REPORT return DiagedMatrix(this); }
  175.  
  176. MatedMatrix BaseMatrix::CopyToMatrix(int nrx, int ncx)
  177. { REPORT return MatedMatrix(this,nrx,ncx); }
  178.  
  179. void GeneralMatrix::operator=(real f)
  180. { REPORT int i=storage; real* s=store; while (i--) { *s++ = f; } }
  181.  
  182. void Matrix::operator=(BaseMatrix& X)
  183. { REPORT CheckConversion(X); Eq(X,MatrixType::Rect); } 
  184.  
  185. void RowVector::operator=(BaseMatrix& X)
  186. {
  187.    REPORT CheckConversion(X); Eq(X,MatrixType::RowV);
  188.    if (nrows!=1) MatrixError("Illegal conversion to row vector");
  189. }
  190.  
  191. void ColumnVector::operator=(BaseMatrix& X)
  192. {
  193.    REPORT CheckConversion(X); Eq(X,MatrixType::ColV);
  194.    if (ncols!=1) MatrixError("Illegal conversion to column vector");
  195. }
  196.  
  197. void SymmetricMatrix::operator=(BaseMatrix& X)
  198. { REPORT CheckConversion(X); Eq(X,MatrixType::Sym); }
  199.  
  200. void UpperTriangularMatrix::operator=(BaseMatrix& X)
  201. { REPORT CheckConversion(X); Eq(X,MatrixType::UT); }
  202.  
  203. void LowerTriangularMatrix::operator=(BaseMatrix& X)
  204. { REPORT CheckConversion(X); Eq(X,MatrixType::LT); }
  205.  
  206. void DiagonalMatrix::operator=(BaseMatrix& X)
  207. { REPORT CheckConversion(X); Eq(X,MatrixType::Diag); }
  208.  
  209. void GeneralMatrix::operator<<(const real* r)
  210. {
  211.    REPORT
  212.    int i = storage; real* s=store;
  213.    while(i--) *s++ = *r++;
  214. }
  215.  
  216.  
  217. /************************* element access *********************************/
  218.  
  219. real& Matrix::element(int m, int n)
  220. {
  221.    if (m<0 || m>= nrows || n<0 || n>= ncols) MatrixError("Index out of range");
  222.    return store[m*ncols+n];
  223. }
  224.  
  225. real& SymmetricMatrix::element(int m, int n)
  226. {
  227.    if (m<0 || n<0 || m >= nrows || n>=ncols) MatrixError("Index out of range");
  228.    if (m>=n) return store[tristore(m)+n];
  229.    else return store[tristore(n)+m];
  230. }
  231.  
  232. real& UpperTriangularMatrix::element(int m, int n)
  233. {
  234.    if (m<0 || n<m || n>=ncols) MatrixError("Index out of range");
  235.    return store[m*ncols+n-tristore(m)];
  236. }
  237.  
  238. real& LowerTriangularMatrix::element(int m, int n)
  239. {
  240.    if (n<0 || m<n || m>=nrows) MatrixError("Index out of range");
  241.    return store[tristore(m)+n];
  242. }
  243.  
  244. real& DiagonalMatrix::element(int m, int n)
  245. {
  246.    if (n<0 || m!=n || m>=nrows || n>=ncols) MatrixError("Index out of range");
  247.    return store[n];
  248. }
  249.  
  250. real& DiagonalMatrix::element(int m)
  251. {
  252.    if (m<0 || m>=nrows) MatrixError("Index out of range");
  253.    return store[m];
  254. }
  255.  
  256. real& ColumnVector::element(int m)
  257. {
  258.    if (m<0 || m>= nrows) MatrixError("Index out of range");
  259.    return store[m];
  260. }
  261.  
  262. real& RowVector::element(int n)
  263. {
  264.    if (n<0 || n>= ncols) MatrixError("Index out of range");
  265.    return store[n];
  266. }
  267.  
  268.