home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / misc / volume21 / newmat02 / part02 / newmat6.cxx < prev    next >
C/C++ Source or Header  |  1991-08-01  |  11KB  |  394 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 "boolean.hxx"
  8.  
  9. typedef double real;                 // all floating point double
  10.  
  11. #include "newmat.hxx"
  12.  
  13. #include "newmatrc.hxx"
  14.  
  15.  
  16. //#define REPORT { static ExeCounter ExeCount(__LINE__,6); ExeCount++; }
  17.  
  18. #define REPORT {}
  19.  
  20. /*************************** general utilities *************************/
  21.  
  22. static int tristore(int n)                      // els in triangular matrix
  23. { return (n*(n+1))/2; }
  24.  
  25.  
  26. /****************************** operators *******************************/
  27.  
  28. real& Matrix::operator()(int m, int n)
  29. {
  30.    if (m<=0 || m>nrows || n<=0 || n>ncols) MatrixError("Index out of range");
  31.    return store[(m-1)*ncols+n-1];
  32. }
  33.  
  34. real& SymmetricMatrix::operator()(int m, int n)
  35. {
  36.    if (m<=0 || n<=0 || m>nrows || n>ncols) MatrixError("Index out of range");
  37.    if (m>=n) return store[tristore(m-1)+n-1];
  38.    else return store[tristore(n-1)+m-1];
  39. }
  40.  
  41. real& UpperTriangularMatrix::operator()(int m, int n)
  42. {
  43.    if (m<=0 || n<m || n>ncols) MatrixError("Index out of range");
  44.    return store[(m-1)*ncols+n-1-tristore(m-1)];
  45. }
  46.  
  47. real& LowerTriangularMatrix::operator()(int m, int n)
  48. {
  49.    if (n<=0 || m<n || m>nrows) MatrixError("Index out of range");
  50.    return store[tristore(m-1)+n-1];
  51. }
  52.  
  53. real& DiagonalMatrix::operator()(int m, int n)
  54. {
  55.    if (n<=0 || m!=n || m>nrows || n>ncols) MatrixError("Index out of range");
  56.    return store[n-1];
  57. }
  58.  
  59. real& DiagonalMatrix::operator()(int m)
  60. {
  61.    if (m<=0 || m>nrows) MatrixError("Index out of range");
  62.    return store[m-1];
  63. }
  64.  
  65. real& ColumnVector::operator()(int m)
  66. {
  67.    if (m<=0 || m> nrows) MatrixError("Index out of range");
  68.    return store[m-1];
  69. }
  70.  
  71. real& RowVector::operator()(int n)
  72. {
  73.    if (n<=0 || n> ncols) MatrixError("Index out of range");
  74.    return store[n-1];
  75. }
  76.  
  77. #ifndef __ZTC__
  78.  
  79. real Matrix::operator()(int m, int n) const
  80. {
  81.    if (m<=0 || m>nrows || n<=0 || n>ncols) MatrixError("Index out of range");
  82.    return store[(m-1)*ncols+n-1];
  83. }
  84.  
  85. real SymmetricMatrix::operator()(int m, int n) const
  86. {
  87.    if (m<=0 || n<=0 || m>nrows || n>ncols) MatrixError("Index out of range");
  88.    if (m>=n) return store[tristore(m-1)+n-1];
  89.    else return store[tristore(n-1)+m-1];
  90. }
  91.  
  92. real UpperTriangularMatrix::operator()(int m, int n) const
  93. {
  94.    if (m<=0 || n<m || n>ncols) MatrixError("Index out of range");
  95.    return store[(m-1)*ncols+n-1-tristore(m-1)];
  96. }
  97.  
  98. real LowerTriangularMatrix::operator()(int m, int n) const
  99. {
  100.    if (n<=0 || m<n || m>nrows) MatrixError("Index out of range");
  101.    return store[tristore(m-1)+n-1];
  102. }
  103.  
  104. real DiagonalMatrix::operator()(int m, int n) const
  105. {
  106.    if (n<=0 || m!=n || m>nrows || n>ncols) MatrixError("Index out of range");
  107.    return store[n-1];
  108. }
  109.  
  110. real DiagonalMatrix::operator()(int m) const
  111. {
  112.    if (m<=0 || m>nrows) MatrixError("Index out of range");
  113.    return store[m-1];
  114. }
  115.  
  116. real ColumnVector::operator()(int m) const
  117. {
  118.    if (m<=0 || m> nrows) MatrixError("Index out of range");
  119.    return store[m-1];
  120. }
  121.  
  122. real RowVector::operator()(int n) const
  123. {
  124.    if (n<=0 || n> ncols) MatrixError("Index out of range");
  125.    return store[n-1];
  126. }
  127.  
  128. #endif
  129.  
  130.  
  131. AddedMatrix BaseMatrix::operator+(BaseMatrix& bm)
  132. { REPORT return AddedMatrix(this, &bm); }
  133.  
  134. MultipliedMatrix BaseMatrix::operator*(BaseMatrix& bm)
  135. { REPORT return MultipliedMatrix(this, &bm); }
  136.  
  137. SolvedMatrix InvertedMatrix::operator*(BaseMatrix& bmx)
  138. { REPORT return SolvedMatrix(bm, &bmx); }
  139.  
  140. SubtractedMatrix BaseMatrix::operator-(BaseMatrix& bm)
  141. { REPORT return SubtractedMatrix(this, &bm); }
  142.  
  143. ShiftedMatrix BaseMatrix::operator+(real f)
  144. { REPORT return ShiftedMatrix(this, f); }
  145.  
  146. ScaledMatrix BaseMatrix::operator*(real f)
  147. { REPORT return ScaledMatrix(this, f); }
  148.  
  149. ScaledMatrix BaseMatrix::operator/(real f)
  150. { REPORT return ScaledMatrix(this, 1.0/f); }
  151.  
  152. ShiftedMatrix BaseMatrix::operator-(real f)
  153. { REPORT return ShiftedMatrix(this, -f); }
  154.  
  155. TransposedMatrix BaseMatrix::t() { REPORT return TransposedMatrix(this); }
  156.  
  157. NegatedMatrix BaseMatrix::operator-() { REPORT return NegatedMatrix(this); }
  158.  
  159. InvertedMatrix BaseMatrix::i() { REPORT return InvertedMatrix(this); }
  160.  
  161. ConstMatrix GeneralMatrix::c() const { REPORT return ConstMatrix(this); }
  162.  
  163. RowedMatrix BaseMatrix::CopyToRow() { REPORT return RowedMatrix(this); }
  164.  
  165. ColedMatrix BaseMatrix::CopyToColumn() { REPORT return ColedMatrix(this); }
  166.  
  167. DiagedMatrix BaseMatrix::CopyToDiagonal() { REPORT return DiagedMatrix(this); }
  168.  
  169. MatedMatrix BaseMatrix::CopyToMatrix(int nrx, int ncx)
  170. { REPORT return MatedMatrix(this,nrx,ncx); }
  171.  
  172. GetSubMatrix BaseMatrix::SubMatrix(int first_row, int last_row, int first_col,
  173.    int last_col)
  174. {
  175.    REPORT
  176.    int a = first_row - 1; int b = last_row - first_row + 1;
  177.    int c = first_col - 1; int d = last_col - first_col + 1;
  178.    if (a<0 || b<=0 || c<0 || d<=0) MatrixError("SubMatrix dimension error");
  179.    return GetSubMatrix(this, a, b, c, d, Type().sub());
  180. }
  181.  
  182. GetSubMatrix BaseMatrix::SymSubMatrix(int first_row, int last_row)
  183. {
  184.    REPORT
  185.    int a = first_row - 1; int b = last_row - first_row + 1;
  186.    if (a<0 || b<=0) MatrixError("SubMatrix dimension error");
  187.    return GetSubMatrix( this, a, b, a, b, Type().ssub());
  188. }
  189.  
  190. GetSubMatrix BaseMatrix::Row(int first_row)
  191. {
  192.    REPORT
  193.    int a = first_row - 1;
  194.    if (a<0) MatrixError("SubMatrix dimension error");
  195.    return GetSubMatrix(this, a, 1, 0, NcolsV(), MatrixType::RowV);
  196. }
  197.  
  198. GetSubMatrix BaseMatrix::Rows(int first_row, int last_row)
  199. {
  200.    REPORT
  201.    int a = first_row - 1; int b = last_row - first_row + 1;
  202.    if (a<0 || b<=0) MatrixError("SubMatrix dimension error");
  203.    return GetSubMatrix(this, a, b, 0, NcolsV(), Type().sub());
  204. }
  205.  
  206. GetSubMatrix BaseMatrix::Col(int first_col)
  207. {
  208.    REPORT
  209.    int c = first_col - 1;
  210.    if (c<0) MatrixError("SubMatrix dimension error");
  211.    return GetSubMatrix(this, 0, NrowsV(), c, 1, MatrixType::ColV);
  212. }
  213.  
  214. GetSubMatrix BaseMatrix::Cols(int first_col, int last_col)
  215. {
  216.    REPORT
  217.    int c = first_col - 1; int d = last_col - first_col + 1;
  218.    if (c<0 || d<=0) MatrixError("SubMatrix dimension error");
  219.    return GetSubMatrix(this, 0, NrowsV(), c, d, Type().sub());
  220. }
  221.  
  222. void GeneralMatrix::operator=(real f)
  223. { REPORT int i=storage; real* s=store; while (i--) { *s++ = f; } }
  224.  
  225. void Matrix::operator=(BaseMatrix& X)
  226. { REPORT CheckConversion(X); Eq(X,MatrixType::Rect); } 
  227.  
  228. void RowVector::operator=(BaseMatrix& X)
  229. {
  230.    REPORT CheckConversion(X); Eq(X,MatrixType::RowV);
  231.    if (nrows!=1) MatrixError("Illegal conversion to row vector");
  232. }
  233.  
  234. void ColumnVector::operator=(BaseMatrix& X)
  235. {
  236.    REPORT CheckConversion(X); Eq(X,MatrixType::ColV);
  237.    if (ncols!=1) MatrixError("Illegal conversion to column vector");
  238. }
  239.  
  240. void SymmetricMatrix::operator=(BaseMatrix& X)
  241. { REPORT CheckConversion(X); Eq(X,MatrixType::Sym); }
  242.  
  243. void UpperTriangularMatrix::operator=(BaseMatrix& X)
  244. { REPORT CheckConversion(X); Eq(X,MatrixType::UT); }
  245.  
  246. void LowerTriangularMatrix::operator=(BaseMatrix& X)
  247. { REPORT CheckConversion(X); Eq(X,MatrixType::LT); }
  248.  
  249. void DiagonalMatrix::operator=(BaseMatrix& X)
  250. { REPORT CheckConversion(X); Eq(X,MatrixType::Diag); }
  251.  
  252. void GeneralMatrix::operator<<(const real* r)
  253. {
  254.    REPORT
  255.    int i = storage; real* s=store;
  256.    while(i--) *s++ = *r++;
  257. }
  258.  
  259.  
  260. /************************* element access *********************************/
  261.  
  262. real& Matrix::element(int m, int n)
  263. {
  264.    if (m<0 || m>= nrows || n<0 || n>= ncols) MatrixError("Index out of range");
  265.    return store[m*ncols+n];
  266. }
  267.  
  268. real& SymmetricMatrix::element(int m, int n)
  269. {
  270.    if (m<0 || n<0 || m >= nrows || n>=ncols) MatrixError("Index out of range");
  271.    if (m>=n) return store[tristore(m)+n];
  272.    else return store[tristore(n)+m];
  273. }
  274.  
  275. real& UpperTriangularMatrix::element(int m, int n)
  276. {
  277.    if (m<0 || n<m || n>=ncols) MatrixError("Index out of range");
  278.    return store[m*ncols+n-tristore(m)];
  279. }
  280.  
  281. real& LowerTriangularMatrix::element(int m, int n)
  282. {
  283.    if (n<0 || m<n || m>=nrows) MatrixError("Index out of range");
  284.    return store[tristore(m)+n];
  285. }
  286.  
  287. real& DiagonalMatrix::element(int m, int n)
  288. {
  289.    if (n<0 || m!=n || m>=nrows || n>=ncols) MatrixError("Index out of range");
  290.    return store[n];
  291. }
  292.  
  293. real& DiagonalMatrix::element(int m)
  294. {
  295.    if (m<0 || m>=nrows) MatrixError("Index out of range");
  296.    return store[m];
  297. }
  298.  
  299. real& ColumnVector::element(int m)
  300. {
  301.    if (m<0 || m>= nrows) MatrixError("Index out of range");
  302.    return store[m];
  303. }
  304.  
  305. real& RowVector::element(int n)
  306. {
  307.    if (n<0 || n>= ncols) MatrixError("Index out of range");
  308.    return store[n];
  309. }
  310.  
  311.  
  312. /****************************** submatrices *********************************/
  313.  
  314. void GetSubMatrix::operator<<(BaseMatrix& bmx)
  315. {
  316.    REPORT
  317.    GeneralMatrix* gmx = bmx.Evaluate();
  318.    GeneralMatrix* gm = bm->Evaluate();
  319.    if ((BaseMatrix*)gm!=bm) MatrixError("Illegal argument for SubMatrix");
  320.    if (row_number != gmx->Nrows() || col_number != gmx->Ncols())
  321.       MatrixError("Dimensions don't match");
  322.    if (row_skip+row_number > gm->Nrows() || col_skip+col_number > gm->Ncols())
  323.       MatrixError("SubMatrix dimension error");
  324.    int i = row_number;
  325.    MatrixRow mrx(gmx, LoadOnEntry); 
  326.    MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
  327.                                   // do need LoadOnEntry
  328.    MatrixRowCol sub;
  329.    while (i--)
  330.    {
  331.       mr.SubRowCol(sub, col_skip, col_number);   // put values in sub
  332.       sub.Copy(mrx); mr.Next(); mrx.Next();
  333.    }
  334.    gmx->tDelete();
  335. }   
  336.  
  337. void GetSubMatrix::operator<<(const real* r)
  338. {
  339.    REPORT
  340.    GeneralMatrix* gm = bm->Evaluate();
  341.    if ((BaseMatrix*)gm!=bm) MatrixError("Illegal argument for SubMatrix");
  342.    if (row_skip+row_number > gm->Nrows() || col_skip+col_number > gm->Ncols())
  343.       MatrixError("SubMatrix dimension error");
  344.    int i = row_number;
  345.    MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
  346.                                   // do need LoadOnEntry
  347.    MatrixRowCol sub;
  348.    while (i--)
  349.    {
  350.       mr.SubRowCol(sub, col_skip, col_number);   // put values in sub
  351.       sub.Copy(r); mr.Next();
  352.    }
  353. }   
  354.  
  355. void GetSubMatrix::operator<<(real r)
  356. {
  357.    REPORT
  358.    GeneralMatrix* gm = bm->Evaluate();
  359.    if ((BaseMatrix*)gm!=bm) MatrixError("Illegal argument for SubMatrix");
  360.    if (row_skip+row_number > gm->Nrows() || col_skip+col_number > gm->Ncols())
  361.       MatrixError("SubMatrix dimension error");
  362.    int i = row_number;
  363.    MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
  364.                                   // do need LoadOnEntry
  365.    MatrixRowCol sub;
  366.    while (i--)
  367.    {
  368.       mr.SubRowCol(sub, col_skip, col_number);   // put values in sub
  369.       sub.Copy(r); mr.Next();
  370.    }
  371. }   
  372.  
  373. void GetSubMatrix::Inject(const GeneralMatrix& gmx)
  374. {
  375.    REPORT
  376.    GeneralMatrix* gm = bm->Evaluate();
  377.    if ((BaseMatrix*)gm!=bm) MatrixError("Illegal argument for SubMatrix");
  378.    if (row_number != gmx.Nrows() || col_number != gmx.Ncols())
  379.       MatrixError("Dimensions don't match");
  380.    if (row_skip+row_number > gm->Nrows() || col_skip+col_number > gm->Ncols())
  381.       MatrixError("SubMatrix dimension error");
  382.    int i = row_number;
  383.    MatrixRow mrx((GeneralMatrix*)(&gmx), LoadOnEntry);
  384.    MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart, row_skip);
  385.                                   // do need LoadOnEntry
  386.    MatrixRowCol sub;
  387.    while (i--)
  388.    {
  389.       mr.SubRowCol(sub, col_skip, col_number);   // put values in sub
  390.       sub.Inject(mrx); mr.Next(); mrx.Next();
  391.    }
  392. }
  393.  
  394.