CYCLUS
l_matrix.cc
Go to the documentation of this file.
1 //-----------------------------------------------------------------------------
2 // A LMatrix object contains n rows and m columns. When a LMatrix object is
3 // created, it is initialized with zeros in all of the elements by default.
4 // To change the value of the element aij at row i and column j, use the
5 // SetElement(int i, int j, long double aij) function. The only difference
6 // between LMatrix and Matrix is that LMatrix elements are long doubles instead
7 // of doubles.
8 //
9 // The function NumRows() returns the number of rows n in the LMatrix object,
10 // and the function NumCols() returns the number of columns m in the LMatrix
11 // object. To access the element in the ith row and jth column of a LMatrix A,
12 // use the syntax A(i,j). To print the matrix, use the Print() function.
13 //
14 // Mathematical functions that can be performed on two LMatrix objects include:
15 //
16 // * matrix assignment: A = B
17 // * matrix addition: A + B
18 // * matrix subtraction: A - B
19 // * matrix multiplication: A * B
20 //
21 // Mathematical functions that can be performed on a LMatrix object A and a
22 // double k include:
23 //
24 // * powers of a matrix: A ^ k
25 // * scalar multiplication: k * A or A * k
26 //
27 // The non-member function identity(int n) creates and returns an nxn identity
28 // matrix.
29 //
30 // Note: when referring to the elements of a LMatrix object, the indices for the
31 // rows and columns start from 1.
32 //-----------------------------------------------------------------------------
33
34 #include "l_matrix.h"
35
36 #include <iostream>
37 #include <iomanip>
38
39 namespace cyclus {
40
41 // constructs a 1x1 matrix of zeroes
43  rows_ = 1; // sets number of rows
44  cols_ = 1; // sets number of columns
45  std::vector<long double> element(1); // creates a single element
46  M_.push_back(element); // adds element to the matrix
47 }
48
49 // constructs an nxm matrix of zeroes
50 LMatrix::LMatrix(int n, int m) {
51  rows_ = n; // sets number of rows
52  cols_ = m; // sets number of columns
53  std::vector<long double> row(m); // creates a row with m elements
54  M_.assign(n, row); // adds n rows_ to the matrix
55 }
56
57 // returns the number of rows n in the matrix
58 int LMatrix::NumRows() const {
59  return rows_;
60 }
61
62 // returns the number of columns m in the matrix
63 int LMatrix::NumCols() const {
64  return cols_;
65 }
66
67 // overloads the () operator so that A(i,j) will return a reference to
68 // the element aij
69 const long double& LMatrix::operator()(int i, int j) const {
70  return M_[i - 1][j - 1];
71 }
72
73 // sets the value for the element aij at row i and column j
74 void LMatrix::SetElement(int i, int j, long double aij) {
75  M_[i - 1][j - 1] = aij;
76 }
77
78 // overloads the () operator so that A(i,j) will write the element aij
79 long double& LMatrix::operator()(int i, int j) {
80  return M_[i - 1][j - 1];
81 }
82
83 // adds a row at the end of the Matrix if it contains the same number of
84 // elements as the number of columns in the Matrix
85 void LMatrix::AddRow(std::vector<long double> row) {
86  int size = row.size();
87  if (size == cols_) {
88  M_.push_back(row);
89  ++rows_; // increase the number of rows
90  }
91 }
92
93 // prints the matrix to standard output
94 void LMatrix::Print() const {
95  std::cout.setf(std::ios::showpoint);
96  std::cout.setf(std::ios::scientific);
97
98  // sets elements to display 6 decimal places
99  std::cout << std::setiosflags(std::ios::fixed) << std::setprecision(10);
100
101  // prints single element if M is a 1x1 matrix
102  if (M_.capacity() == 1 && M_[0].capacity() == 1) {
103  std::cout << "[ " << M_[0][0] << " ]" << std::endl;
104  } else {
105  // loops through the rows of the matrix
106  for (int i = 0; i < rows_; i++) {
107  // prints all of the elements in the ith row of the matrix
108  std::cout << "[";
109  for (int j = 0; j < cols_ - 1; j++) {
110  std::cout << " " << std::setw(9) << M_[i][j] << " ";
111  }
112  std::cout << std::setw(9) << M_[i][cols_ - 1] << " ]" << std::endl;
113  }
114  }
115 }
116
117 // overloads the assignment operator "A = B" for matrix objects
118 const LMatrix& LMatrix::operator=(const LMatrix& rhs) {
119  if (this == &rhs) {
120  return *this; // returns A if it is already the same matrix as B
121  } else {
122  rows_ = rhs.rows_; // resets the number of rows to match B
123  cols_ = rhs.cols_; // resets the number of columns to match B
124
125  // rebuilds A with the dimensions of B
126  std::vector<long double> row(cols_); // creates a row with m elements
127  M_.assign(rows_, row);
128
129  // copies all of the elements from B into A
130  for (int i = 0; i < rows_; i++) {
131  for (int j = 0; j < cols_; j++) {
132  M_[i][j] = rhs.M_[i][j];
133  }
134  }
135
136  return *this; // returns the new matrix "A = B"
137  }
138 }
139
140 // overloads the assignment operator "A = A + B" for matrix objects
141 // Note: if the matrix dimensions do not match, then A is returned unchanged
142 const LMatrix& LMatrix::operator+=(const LMatrix& rhs) {
143  if (this->rows_ == rhs.rows_ && this->cols_ == rhs.cols_) {
144  // performs matrix addition and stores the result in A
145  for (int i = 0; i < rows_; i++) {
146  for (int j = 0; j < cols_; j++) {
147  M_[i][j] += rhs.M_[i][j];
148  }
149  }
150  }
151
152  return *this; // returns the new matrix "A = A + B"
153 }
154
155 // overloads the assignment operator "A = A - B" for matrix objects
156 // Note: if the matrix dimensions do not match, then A is returned unchanged
157 const LMatrix& LMatrix::operator-=(const LMatrix& rhs) {
158  if (this->rows_ == rhs.rows_ && this->cols_ == rhs.cols_) {
159  // performs matrix subtraction and stores the result in A
160  for (int i = 0; i < rows_; i++) {
161  for (int j = 0; j < cols_; j++) {
162  M_[i][j] -= rhs.M_[i][j];
163  }
164  }
165  }
166
167  return *this; // returns the new matrix "A = A - B"
168 }
169
170 // overloads the assignment operator "A = A * B" for matrix objects
171 //
172 // Note: This function will perform matrix multiplication only if the number
173 // of columns in A is equal to the number of rows in B. If the matrices
174 // cannot be multipled due to having incorrect dimensions for matrix
175 // multiplication to be defined, then A will be returned unchanged
176 const LMatrix& LMatrix::operator*=(const LMatrix& rhs) {
177  long double aij; // stores temporary matrix element for row i and column j
178
179  // creates a new matrix called temp with the number of rows of A and the
180  // number of columns of B
181  LMatrix temp(this->rows_, rhs.cols_);
182
183  // performs matrix multiplication if the number of columns of A equals the
184  // number of rows of B
185  if (this->cols_ == rhs.rows_) {
186  // multiplies row vector i of A by the column vector j of B and stores
187  // the result in the new matrix temp as element aij
188  for (int i = 0; i < rows_; i++) {
189  for (int j = 0; j < rhs.cols_; j++) {
190  aij = 0; // resets the sum for the next element
191
192  // sums the product of the row vector from A and the column vector
193  // from B using the formula aij = ai1*b1j + ai2*b2j + ...
194  for (int k = 0; k < rows_; k++) {
195  aij = aij + M_[i][k] * rhs.M_[k][j];
196  }
197
198  temp.M_[i][j] = aij; // copies element aij into temp
199  }
200  }
201
202  *this = temp; // sets the new matrix "A = A * B"
203  }
204
205  return *this; // returns the new matrix "A = A * B"
206 }
207
208 // overloads the arithmetic operator "A + B" for matrix objects
209
210 LMatrix operator+(const LMatrix& lhs, const LMatrix& rhs) {
211  LMatrix ans(lhs);
212  ans += rhs;
213  return ans; // returns the resulting matrix A + B
214 }
215
216 // overloads the arithmetic operator "A - B" for matrix objects
217
218 LMatrix operator-(const LMatrix& lhs, const LMatrix& rhs) {
219  LMatrix ans(lhs);
220  ans -= rhs;
221  return ans; // returns the resulting matrix A - B
222 }
223
224 // overloads the arithmetic operator "A * B" for matrix objects
225
226 LMatrix operator*(const LMatrix& lhs, const LMatrix& rhs) {
227  LMatrix ans(lhs);
228  ans *= rhs;
229  return ans; // returns the resulting matrix A * B
230 }
231
232 // friend of the Matrix class that performs scalar multiplication k * A
233
234 LMatrix operator*(const long double k, const LMatrix& A) {
235  LMatrix ans(A); // copies A into a new matrix called ans
236  long double aij;
237
238  // multiplies every element in the matrix A by the scalar k
239  for (int i = 0; i < A.rows_; i++) {
240  for (int j = 0; j < A.cols_; j++) {
241  aij = A(i + 1, j + 1) * k;
242  ans.M_[i][j] = aij;
243  }
244  }
245
246  return ans; // returns the resulting matrix k * A
247 }
248
249 // friend of the Matrix class that performs scalar multiplication A * k
250
251 LMatrix operator*(const LMatrix& A, const long double k) {
252  LMatrix ans(A); // copies A into a new matrix called ans
253  long double aij;
254
255  // multiplies every element in the matrix A by the scalar k
256  for (int i = 0; i < A.rows_; i++) {
257  for (int j = 0; j < A.cols_; j++) {
258  aij = A(i + 1, j + 1) * k;
259  ans.M_[i][j] = aij;
260  }
261  }
262
263  return ans; // returns the resulting matrix A * k
264 }
265
266 // friend of the Matrix class that calculates powers of a square matrix A ^ k
267 // Note: if the matrix is not square, then A is returned unchanged
268
269 LMatrix operator^(const LMatrix& A, const int k) {
270  LMatrix ans(A); // copies A into a new matrix called ans
271
272  // multiplies A by itself k times if the matrix is square
273  if (A.cols_ == A.rows_) {
274  for (int i = 0; i < k - 1; i++) {
275  ans *= A;
276  }
277  }
278
279  return ans; // returns the resulting matrix A ^ k
280 }
281
282 // creates and returns an nxn identity matrix I
283
285  LMatrix I(n, n);
286
287  for (int i = 0; i < n ; i++) {
288  I.SetElement(i + 1, i + 1, 1);
289  }
290
291  return I; // returns the nxn identity matrix
292 }
293
294 } // namespace cyclus
void SetElement(int i, int j, long double aij)
Definition: l_matrix.cc:74
void Print() const
Definition: l_matrix.cc:94
const LMatrix & operator-=(const LMatrix &rhs)
Definition: l_matrix.cc:157
const long double & operator()(int i, int j) const
Definition: l_matrix.cc:69
LMatrix operator+(const LMatrix &lhs, const LMatrix &rhs)
Definition: l_matrix.cc:210
LMatrix operator-(const LMatrix &lhs, const LMatrix &rhs)
Definition: l_matrix.cc:218
int NumRows() const
Definition: l_matrix.cc:58
const LMatrix & operator*=(const LMatrix &rhs)
Definition: l_matrix.cc:176
void AddRow(std::vector< long double > row)
Definition: l_matrix.cc:85
LMatrix identity(int n)
Definition: l_matrix.cc:284
taken directly from OsiSolverInterface.cpp on 2/17/14 from https://projects.coin-or.org/Osi/browser/trunk.
Definition: agent.cc:14
friend LMatrix operator^(const LMatrix &A, const int k)
Definition: l_matrix.cc:269
const LMatrix & operator+=(const LMatrix &rhs)
Definition: l_matrix.cc:142
const LMatrix & operator=(const LMatrix &rhs)
Definition: l_matrix.cc:118
int NumCols() const
Definition: l_matrix.cc:63
friend LMatrix operator*(const long double k, const LMatrix &A)
Definition: l_matrix.cc:234