CoolProp  6.6.1dev
An open-source fluid property and humid air property database
MatrixMath.h
Go to the documentation of this file.
1 #ifndef MATRIXMATH_H
2 #define MATRIXMATH_H
3 
4 #include "CoolPropTools.h"
5 #include "Exceptions.h"
6 
7 #include <vector>
8 #include <string>
9 #include <numeric> // inner_product
10 #include <sstream>
11 #include "float.h"
12 
13 #include "Eigen/Core"
14 
16 
22 template <size_t dimcount, typename T>
23 struct VectorNd
24 {
25  typedef std::vector<typename VectorNd<dimcount - 1, T>::type> type;
26 };
27 template <typename T>
28 struct VectorNd<0, T>
29 {
30  typedef T type;
31 };
32 
33 namespace CoolProp {
34 
36 template <class T>
37 std::size_t num_rows(std::vector<T> const& in) {
38  return in.size();
39 }
40 template <class T>
41 std::size_t num_rows(std::vector<std::vector<T>> const& in) {
42  return in.size();
43 }
44 
45 template <class T>
46 std::size_t max_cols(std::vector<std::vector<T>> const& in) {
47  std::size_t cols = 0;
48  std::size_t col = 0;
49  for (std::size_t i = 0; i < in.size(); i++) {
50  col = in[i].size();
51  if (cols < col) {
52  cols = col;
53  }
54  }
55  return cols;
56 };
57 template <class T>
58 bool is_squared(std::vector<std::vector<T>> const& in) {
59  std::size_t cols = max_cols(in);
60  if (cols != num_rows(in)) {
61  return false;
62  } else {
63  for (std::size_t i = 0; i < in.size(); i++) {
64  if (cols != in[i].size()) {
65  return false;
66  }
67  }
68  }
69  return true;
70 };
71 
72 template <class T>
73 std::size_t num_cols(std::vector<T> const& in) {
74  return 1;
75 }
76 template <class T>
77 std::size_t num_cols(std::vector<std::vector<T>> const& in) {
78  if (num_rows(in) > 0) {
79  if (is_squared(in)) {
80  return in[0].size();
81  } else {
82  return max_cols(in);
83  }
84  } else {
85  return 0;
86  }
87 };
88 
90 
96 template <typename T>
97 std::vector<T> eigen_to_vec1D(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& coefficients, int axis = 0) {
98  std::vector<T> result;
99  size_t r = coefficients.rows(), c = coefficients.cols();
100  if (axis == 0) {
101  if (c != 1) throw ValueError(format("Your matrix has the wrong dimensions: %d,%d", r, c));
102  result.resize(r);
103  for (size_t i = 0; i < r; ++i) {
104  result[i] = coefficients(i, 0);
105  }
106  } else if (axis == 1) {
107  if (r != 1) throw ValueError(format("Your matrix has the wrong dimensions: %d,%d", r, c));
108  result.resize(c);
109  for (size_t i = 0; i < c; ++i) {
110  result[i] = coefficients(0, i);
111  }
112  } else {
113  throw ValueError(format("You have to provide axis information: %d is not valid. ", axis));
114  }
115  return result;
116 }
118 template <class T>
119 std::vector<std::vector<T>> eigen_to_vec(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& coefficients) {
120  // Eigen uses columns as major axis, this might be faster than the row iteration.
121  // However, the 2D vector stores things differently, no idea what is faster...
122  std::vector<std::vector<T>> result;
123  size_t r = coefficients.rows(), c = coefficients.cols();
124  result.resize(r, std::vector<T>(c, 0)); // extends vector if necessary
125  for (size_t i = 0; i < r; ++i) {
126  result[i].resize(c, 0);
127  for (size_t j = 0; j < c; ++j) {
128  result[i][j] = coefficients(i, j);
129  }
130  }
131  return result;
132 }
133 
135 template <class T>
136 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> vec_to_eigen(const std::vector<std::vector<T>>& coefficients) {
137  size_t nRows = num_rows(coefficients), nCols = num_cols(coefficients);
138  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> result(nRows, nCols);
139  for (size_t i = 0; i < nCols; ++i) {
140  for (size_t j = 0; j < nRows; ++j) {
141  result(j, i) = coefficients[j][i];
142  }
143  }
144  return result;
145 }
150 template <class T>
151 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> vec_to_eigen(const std::vector<T>& coefficients, int axis = 0) {
152  size_t nRows = num_rows(coefficients);
153  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> result;
154  if (axis == 0)
155  result.resize(nRows, 1);
156  else if (axis == 1)
157  result.resize(1, nRows);
158  else
159  throw ValueError(format("You have to provide axis information: %d is not valid. ", axis));
160  for (size_t i = 0; i < nRows; ++i) {
161  if (axis == 0) result(i, 0) = coefficients[i];
162  if (axis == 1) result(0, i) = coefficients[i];
163  }
164  return result;
165 }
167 template <class T>
168 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> vec_to_eigen(const T& coefficient) {
169  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> result = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>(1, 1);
170  result(0, 0) = coefficient;
171  return result;
172 }
173 
175 
180 template <class T>
181 Eigen::Matrix<T, Eigen::Dynamic, 1> makeColVector(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& matrix) {
182  std::size_t r = matrix.rows();
183  std::size_t c = matrix.cols();
184  Eigen::Matrix<T, Eigen::Dynamic, 1> vector;
185  if (r == 1 && c >= 1) { // Check passed, matrix can be transformed
186  vector = matrix.transpose().block(0, 0, c, r);
187  } else if (r >= 1 && c == 1) { // Check passed, matrix can be transformed
188  vector = matrix.block(0, 0, r, c);
189  } else { // Check failed, throw error
190  throw ValueError(format("Your matrix (%d,%d) cannot be converted into a vector (x,1).", r, c));
191  }
192  return vector;
193 }
194 template <class T>
195 Eigen::Matrix<T, Eigen::Dynamic, 1> makeVector(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& matrix) {
196  return makeColVector(matrix);
197 }
198 template <class T>
199 Eigen::Matrix<T, 1, Eigen::Dynamic> makeRowVector(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& matrix) {
200  std::size_t r = matrix.rows();
201  std::size_t c = matrix.cols();
202  Eigen::Matrix<T, 1, Eigen::Dynamic> vector;
203  if (r == 1 && c >= 1) { // Check passed, matrix can be transformed
204  vector = matrix.block(0, 0, r, c);
205  } else if (r >= 1 && c == 1) { // Check passed, matrix can be transformed
206  vector = matrix.transpose().block(0, 0, c, r);
207  } else { // Check failed, throw error
208  throw ValueError(format("Your matrix (%d,%d) cannot be converted into a vector (1,x).", r, c));
209  }
210  return vector;
211 }
212 
214 
217 template <class T>
218 void removeRow(Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& matrix, std::size_t rowToRemove) {
219  //template<class T> void removeRow(Eigen::MatrixXd& matrix, std::size_t rowToRemove){
220  //void removeRow(Eigen::MatrixXd& matrix, unsigned int rowToRemove){
221  //template <typename Derived> void removeRow(Eigen::MatrixBase<Derived> &matrix, std::size_t rowToRemove){
222  std::size_t numRows = matrix.rows() - 1;
223  std::size_t numCols = matrix.cols();
224  if (rowToRemove < numRows) {
225  matrix.block(rowToRemove, 0, numRows - rowToRemove, numCols) = matrix.block(rowToRemove + 1, 0, numRows - rowToRemove, numCols);
226  } else {
227  if (rowToRemove > numRows) {
228  throw ValueError(format("Your matrix does not have enough rows, %d is not greater or equal to %d.", numRows, rowToRemove));
229  }
230  // Do nothing, resize removes the last row
231  }
232  matrix.conservativeResize(numRows, numCols);
233 }
234 
235 template <class T>
236 void removeColumn(Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& matrix, std::size_t colToRemove) {
237  //template<class T> void removeColumn(Eigen::MatrixXd& matrix, std::size_t colToRemove){
238  //void removeColumn(Eigen::MatrixXd& matrix, unsigned int colToRemove){
239  //template <typename Derived> void removeColumn(Eigen::MatrixBase<Derived> &matrix, std::size_t colToRemove){
240  std::size_t numRows = matrix.rows();
241  std::size_t numCols = matrix.cols() - 1;
242  if (colToRemove < numCols) {
243  matrix.block(0, colToRemove, numRows, numCols - colToRemove) = matrix.block(0, colToRemove + 1, numRows, numCols - colToRemove);
244  } else {
245  if (colToRemove > numCols) {
246  throw ValueError(format("Your matrix does not have enough columns, %d is not greater or equal to %d.", numCols, colToRemove));
247  }
248  // Do nothing, resize removes the last column
249  }
250  matrix.conservativeResize(numRows, numCols);
251 }
252 
254 //template <class T> Eigen::Matrix<T, Eigen::Dynamic,Eigen::Dynamic> convert(const std::vector<std::vector<T> > &coefficients){
255 // size_t nRows = num_rows(coefficients), nCols = num_cols(coefficients);
256 // Eigen::MatrixBase<T> result(nRows,nCols);
257 // for (size_t i = 0; i < nCols; ++i) {
258 // for (size_t j = 0; j < nRows; ++j) {
259 // result(j,i) = coefficients[j][i];
260 // }
261 // }
262 // return result;
263 //}
264 //
266 //template <class T, int R, int C> void convert(const std::vector<std::vector<T> > &coefficients, Eigen::Matrix<T,R,C> &result){
267 // size_t nRows = num_rows(coefficients), nCols = num_cols(coefficients);
268 // //Eigen::MatrixBase<T> result(nRows,nCols);
269 // for (size_t i = 0; i < nCols; ++i) {
270 // for (size_t j = 0; j < nRows; ++j) {
271 // result(j,i) = coefficients[j][i];
272 // }
273 // }
274 // //return result;
275 //}
276 
277 //
278 //template <class T> void convert(const std::vector<std::vector<T> > &coefficients, Eigen::MatrixBase<T> &result){
279 // size_t nRows = num_rows(coefficients), nCols = num_cols(coefficients);
280 // //Eigen::MatrixBase<T> result;
281 // //if ((R!=nRows) || (C!=nCols))
282 // result.resize(nRows,nCols);
283 // for (size_t i = 0; i < nCols; ++i) {
284 // for (size_t j = 0; j < nRows; ++j) {
285 // result(j,i) = coefficients[j][i];
286 // }
287 // }
288 // //return result;
289 //}
290 
291 //template <class Derived>
292 //inline void func1(MatrixBase<Derived> &out_mat ){
293 // // Do something then return a matrix
294 // out_mat = ...
295 //}
296 
297 //template <class Derived>
298 //Eigen::Matrix<class Derived::Scalar, Derived::RowsAtCompileTime, Derived::ColsAtCompileTime>
299 //Multiply(const Eigen::MatrixBase<DerivedA>& p1,
300 // const Eigen::MatrixBase<DerivedB>& p2)
301 //{
302 // return (p1 * p2);
303 //}
304 //
305 //
306 //template <typename DerivedA, typename DerivedB>
307 //Eigen::Matrix<typename DerivedA::Scalar, DerivedA::RowsAtCompileTime, DerivedB::ColsAtCompileTime>
308 //Multiply(const Eigen::MatrixBase<DerivedA>& p1,
309 // const Eigen::MatrixBase<DerivedB>& p2)
310 //{
311 // return (p1 * p2);
312 //}
313 
315 static const char* stdFmt = "%8.3f";
316 
318 template <class T>
319 std::string vec_to_string(const std::vector<T>& a, const char* fmt) {
320  if (a.size() < 1) return std::string("");
321  std::stringstream out;
322  out << "[ " << format(fmt, a[0]);
323  for (size_t j = 1; j < a.size(); j++) {
324  out << ", " << format(fmt, a[j]);
325  }
326  out << " ]";
327  return out.str();
328 };
329 template <class T>
330 std::string vec_to_string(const std::vector<T>& a) {
331  return vec_to_string(std::vector<double>(a.begin(), a.end()), stdFmt);
332 };
334 inline std::string stringvec_to_string(const std::vector<std::string>& a) {
335  if (a.size() < 1) return std::string("");
336  std::stringstream out;
337  out << "[ " << format("%s", a[0].c_str());
338  for (size_t j = 1; j < a.size(); j++) {
339  out << ", " << format("%s", a[j].c_str());
340  }
341  out << " ]";
342  return out.str();
343 };
344 
346 template <class T>
347 std::string vec_to_string(const T& a, const char* fmt) {
348  std::vector<T> vec;
349  vec.push_back(a);
350  return vec_to_string(vec, fmt);
351 };
352 template <class T>
353 std::string vec_to_string(const T& a) {
354  return vec_to_string((double)a, stdFmt);
355 };
356 
358 template <class T>
359 std::string vec_to_string(const std::vector<std::vector<T>>& A, const char* fmt) {
360  if (A.size() < 1) return std::string("");
361  std::stringstream out;
362  out << "[ " << vec_to_string(A[0], fmt);
363  for (size_t j = 1; j < A.size(); j++) {
364  out << ", " << std::endl << " " << vec_to_string(A[j], fmt);
365  }
366  out << " ]";
367  return out.str();
368 };
369 template <class T>
370 std::string vec_to_string(const std::vector<std::vector<T>>& A) {
371  return vec_to_string(A, stdFmt);
372 };
373 
375 template <class T>
376 std::string mat_to_string(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& A, const char* fmt) {
377  //std::string mat_to_string(const Eigen::MatrixXd &A, const char *fmt) {
378  std::size_t r = A.rows();
379  std::size_t c = A.cols();
380  if ((r < 1) || (c < 1)) return std::string("");
381  std::stringstream out;
382  out << "[ ";
383  if (r == 1) {
384  out << format(fmt, A(0, 0));
385  for (size_t j = 1; j < c; j++) {
386  out << ", " << format(fmt, A(0, j));
387  }
388  } else {
389  out << mat_to_string(Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>(A.row(0)), fmt);
390  for (size_t i = 1; i < r; i++) {
391  out << ", " << std::endl << " " << mat_to_string(Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>(A.row(i)), fmt);
392  }
393  }
394  out << " ]";
395  return out.str();
396 };
397 template <class T>
398 std::string mat_to_string(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& A) {
399  //std::string vec_to_string(const Eigen::MatrixXd &A) {
400  return mat_to_string(A, stdFmt);
401 };
402 
404 //static const char* stdFmt = "%8.3f";
405 //
407 //template<class T> std::string vec_to_string(const std::vector<T> &a, const char *fmt) {
408 // if (a.size()<1) return std::string("");
409 // std::stringstream out;
410 // out << "[ " << format(fmt,a[0]);
411 // for (size_t j = 1; j < a.size(); j++) {
412 // out << ", " << format(fmt, a[j]);
413 // }
414 // out << " ]";
415 // return out.str();
416 //};
417 //template<class T> std::string vec_to_string(const std::vector<T> &a) {
418 // return vec_to_string(a, stdFmt);
419 //};
420 //
422 //template<class T> std::string vec_to_string(const T &a, const char *fmt) {
423 // std::vector<T> vec;
424 // vec.push_back(a);
425 // return vec_to_string(vec, fmt);
426 //};
427 //template<class T> std::string vec_to_string(const T &a) {
428 // return vec_to_string(a, stdFmt);
429 //};
430 //
432 //template<class T> std::string vec_to_string(const std::vector<std::vector<T> > &A, const char *fmt) {
433 // if (A.size()<1) return std::string("");
434 // std::stringstream out;
435 // out << "[ " << vec_to_string(A[0], fmt);
436 // for (size_t j = 1; j < A.size(); j++) {
437 // out << ", " << std::endl << " " << vec_to_string(A[j], fmt);
438 // }
439 // out << " ]";
440 // return out.str();
441 //};
442 //template<class T> std::string vec_to_string(const std::vector<std::vector<T> > &A) {
443 // return vec_to_string(A, stdFmt);
444 //};
445 //
447 //template <class T> std::string mat_to_string(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> &A, const char *fmt) {
449 // std::size_t r = A.rows();
450 // std::size_t c = A.cols();
451 // if ((r<1)||(c<1)) return std::string("");
452 // std::stringstream out;
453 // out << "[ ";
454 // if (r==1) {
455 // out << format(fmt, A(0,0));
456 // for (size_t j = 1; j < c; j++) {
457 // out << ", " << format(fmt, A(0,j));
458 // }
459 // } else {
460 // out << mat_to_string(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>(A.row(0)), fmt);
461 // for (size_t i = 1; i < r; i++) {
462 // out << ", " << std::endl << " " << mat_to_string(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>(A.row(i)), fmt);
463 // }
464 // }
465 // out << " ]";
466 // return out.str();
467 //};
468 //template <class T> std::string mat_to_string(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> &A) {
470 // return mat_to_string(A, stdFmt);
471 //};
472 
474 //template<class T> std::string vec_to_string(const T &a){
475 // return vec_to_string(a, stdFmt);
476 // std::stringstream out;
477 // out << format("[ %7.3f ]",a);
478 // return out.str();
479 //};
480 //template<class T> std::string vec_to_string(const VectorNd<0, T> &a){
481 // return vec_to_string(a, stdFmt);
482 //};
483 //template<class T> std::string vec_to_string(const VectorNd<0, T> &a, const char *fmt) {
484 // VectorNd<1, T> vec;
485 // vec.push_back(a);
486 // return vec_to_string(vec, fmt);
487 //};
488 //
490 //template<class T> std::string vec_to_string(const VectorNd<1, T> &a) {
491 // return vec_to_string(a, stdFmt);
492 //};
493 //template<class T> std::string vec_to_string(const VectorNd<1, T> &a, const char *fmt) {
494 // if (a.size()<1) {
495 // return std::string("");
496 // } else {
497 // std::stringstream out;
498 // out << "[ ";
499 // out << format(fmt,a[0]);
500 // for (size_t j = 1; j < a.size(); j++) {
501 // out << ", ";
502 // out << format(fmt,a[j]);
503 // }
504 // out << " ]";
505 // return out.str();
506 // }
507 //};
508 //
510 //template<class T> std::string vec_to_string(const VectorNd<2, T> &A) {
511 // return vec_to_string(A, stdFmt);
512 //}
513 //template<class T> std::string vec_to_string(const VectorNd<2, T> &A, const char *fmt) {
514 // if (A.size()<1) return std::string("");
515 // std::stringstream out;
516 // out << "[ " << format(fmt,A[0]);
517 // for (size_t j = 1; j < A.size(); j++) {
518 // out << ", " << std::endl << " " << vec_to_string(A[j], fmt);
519 // }
520 // out << " ]";
521 // return out.str();
522 //}
523 
525 //template<class T> std::vector<T> linsolve(std::vector<std::vector<T> > const& A, std::vector<T> const& b);
526 //template<class T> std::vector<std::vector<T> > linsolve(std::vector<std::vector<T> > const& A, std::vector<std::vector<T> > const& B);
527 //
529 //template<class T> std::size_t num_rows (std::vector<std::vector<T> > const& in);
530 //template<class T> std::size_t num_cols (std::vector<std::vector<T> > const& in);
531 //template<class T> std::size_t max_cols (std::vector<std::vector<T> > const& in);
532 //template<class T> std::vector<T> get_row (std::vector<std::vector<T> > const& in, size_t row);
533 //template<class T> std::vector<T> get_col (std::vector<std::vector<T> > const& in, size_t col);
534 //template<class T> bool is_squared(std::vector<std::vector<T> > const& in);
535 //template<class T> std::vector<std::vector<T> > make_squared(std::vector<std::vector<T> > const& in);
536 //
538 //template<class T> T multiply( std::vector<T> const& A, std::vector<T> const& B);
539 //template<class T> std::vector<T> multiply(std::vector<std::vector<T> > const& A, std::vector<T> const& B);
540 //template<class T> std::vector<std::vector<T> > multiply(std::vector<std::vector<T> > const& A, std::vector<std::vector<T> > const& B);
541 //
542 //template<class T> T dot_product(std::vector<T> const& a, std::vector<T> const& b);
543 //template<class T> std::vector<T> cross_product(std::vector<T> const& a, std::vector<T> const& b);
544 //
545 //template<class T> std::vector<std::vector<T> > transpose(std::vector<std::vector<T> > const& in);
546 //template<class T> std::vector<std::vector<T> > invert(std::vector<std::vector<T> > const& in);
547 //
548 //template<class T> std::string vec_to_string( T const& a);
549 //template<class T> std::string vec_to_string( std::vector<T> const& a);
550 //template<class T> std::string vec_to_string(std::vector<std::vector<T> > const& A);
551 //
552 //template<class T> std::string vec_to_string( std::vector<T> const& a, const char *fmt);
553 //template<class T> std::string vec_to_string(std::vector<std::vector<T> > const& A, const char *fmt);
554 
555 /*
556 Owe a debt of gratitude to http://sole.ooz.ie/en - very clear treatment of GJ
557 */
558 template <typename T>
559 void swap_rows(std::vector<std::vector<T>>* A, size_t row1, size_t row2) {
560  for (size_t col = 0; col < (*A)[0].size(); col++) {
561  std::swap((*A)[row1][col], (*A)[row2][col]);
562  }
563 };
564 template <typename T>
565 void subtract_row_multiple(std::vector<std::vector<T>>* A, size_t row, T multiple, size_t pivot_row) {
566  for (size_t col = 0; col < (*A)[0].size(); col++) {
567  (*A)[row][col] -= multiple * (*A)[pivot_row][col];
568  }
569 };
570 template <typename T>
571 void divide_row_by(std::vector<std::vector<T>>* A, size_t row, T value) {
572  for (size_t col = 0; col < (*A)[0].size(); col++) {
573  (*A)[row][col] /= value;
574  }
575 };
576 
577 template <typename T>
578 size_t get_pivot_row(std::vector<std::vector<T>>* A, size_t col) {
579  std::size_t index = col;
580  T max = 0, val;
581 
582  for (size_t row = col; row < (*A).size(); row++) {
583  val = (*A)[row][col];
584  if (std::abs(val) > max) {
585  max = std::abs(val);
586  index = row;
587  }
588  }
589  return index;
590 };
591 
592 template <typename T>
593 std::vector<std::vector<T>> linsolve_Gauss_Jordan(std::vector<std::vector<T>> const& A, std::vector<std::vector<T>> const& B) {
594  std::vector<std::vector<T>> AB;
595  std::vector<std::vector<T>> X;
596  size_t pivot_row;
597  T pivot_element;
598 
599  size_t NrowA = num_rows(A);
600  size_t NrowB = num_rows(B);
601  size_t NcolA = num_cols(A);
602  size_t NcolB = num_cols(B);
603 
604  if (NrowA != NrowB) throw ValueError(format("You have to provide matrices with the same number of rows: %d is not %d. ", NrowA, NrowB));
605 
606  AB.resize(NrowA, std::vector<T>(NcolA + NcolB, 0));
607  X.resize(NrowA, std::vector<T>(NcolB, 0));
608 
609  // Build the augmented matrix
610  for (size_t row = 0; row < NrowA; row++) {
611  for (size_t col = 0; col < NcolA; col++) {
612  AB[row][col] = A[row][col];
613  }
614  for (size_t col = NcolA; col < NcolA + NcolB; col++) {
615  AB[row][col] = B[row][col - NcolA];
616  }
617  }
618 
619  for (size_t col = 0; col < NcolA; col++) {
620  // Find the pivot value
621  pivot_row = get_pivot_row(&AB, col);
622 
623  if (std::abs(AB[pivot_row][col]) < 10 * DBL_EPSILON) {
624  throw ValueError(format("Zero occurred in row %d, the matrix is singular. ", pivot_row));
625  }
626 
627  if (pivot_row >= col) {
628  // Swap pivot row and current row
629  swap_rows(&AB, col, pivot_row);
630  }
631  // Get the pivot element
632  pivot_element = AB[col][col];
633  // Divide the pivot row by the pivot element
634  divide_row_by(&AB, col, pivot_element);
635 
636  if (col < NrowA - 1) {
637  // All the rest of the rows, subtract the value of the [r][c] combination
638  for (size_t row = col + 1; row < NrowA; row++) {
639  subtract_row_multiple(&AB, row, AB[row][col], col);
640  }
641  }
642  }
643  for (std::size_t col = NcolA - 1; col > 0; col--) {
644  for (int row = static_cast<int>(col) - 1; row >= 0; row--) {
645  subtract_row_multiple(&AB, row, AB[row][col], col);
646  }
647  }
648  // Set the output value
649  for (size_t row = 0; row < NrowA; row++) {
650  for (size_t col = 0; col < NcolB; col++) {
651  X[row][col] = AB[row][NcolA + col];
652  }
653  }
654  return X;
655 };
656 
657 //std::vector<std::vector<double> > linsolve_Gauss_Jordan_reimpl(std::vector<std::vector<double> > const& A, std::vector<std::vector<double> > const& B) {
658 // std::vector<std::vector<double> > AB;
659 // std::vector<std::vector<double> > X;
660 // size_t pivot_row;
661 // double pivot_element;
662 // double tmp_element;
663 //
664 // size_t NrowA = num_rows(A);
665 // size_t NrowB = num_rows(B);
666 // size_t NcolA = num_cols(A);
667 // size_t NcolB = num_cols(B);
668 //
669 // if (NrowA!=NrowB) throw ValueError(format("You have to provide matrices with the same number of rows: %d is not %d. ",NrowA,NrowB));
670 //
671 // AB.resize(NrowA, std::vector<double>(NcolA+NcolB, 0));
672 // X.resize(NrowA, std::vector<double>(NcolB, 0));
673 //
674 // // Build the augmented matrix
675 // for (size_t row = 0; row < NrowA; row++){
676 // for (size_t col = 0; col < NcolA; col++){
677 // AB[row][col] = A[row][col];
678 // }
679 // for (size_t col = NcolA; col < NcolA+NcolB; col++){
680 // AB[row][col] = B[row][col-NcolA];
681 // }
682 // }
683 //
684 // for (size_t col = 0; col < NcolA; col++){
685 // // Find the pivot row
686 // pivot_row = 0;
687 // pivot_element = 0.0;
688 // for (size_t row = col; row < NrowA; row++){
689 // tmp_element = std::abs(AB[row][col]);
690 // if (tmp_element>pivot_element) {
691 // pivot_element = tmp_element;
692 // pivot_row = row;
693 // }
694 // }
695 // // Check for errors
696 // if (AB[pivot_row][col]<1./_HUGE) throw ValueError(format("Zero occurred in row %d, the matrix is singular. ",pivot_row));
697 // // Swap the rows
698 // if (pivot_row>col) {
699 // for (size_t colInt = 0; colInt < NcolA; colInt++){
700 // std::swap(AB[pivot_row][colInt],AB[pivot_row][colInt]);
701 // }
702 // }
703 // // Process the entries below current element
704 // for (size_t row = col; row < NrowA; row++){
705 // // Entries to the right of current element (until end of A)
706 // for (size_t colInt = col+1; colInt < NcolA; colInt++){
707 // // All entries in augmented matrix
708 // for (size_t colFull = col; colFull < NcolA+NcolB; colFull++){
709 // AB[colInt][colFull] -= AB[col][colFull] * AB[colInt][col] / AB[col][col];
710 // }
711 // AB[colInt][col] = 0.0;
712 // }
713 // }
714 // }
715 // return AB;
716 //}
717 
718 template <class T>
719 std::vector<std::vector<T>> linsolve(std::vector<std::vector<T>> const& A, std::vector<std::vector<T>> const& B) {
720  return linsolve_Gauss_Jordan(A, B);
721 };
722 
723 template <class T>
724 std::vector<T> linsolve(std::vector<std::vector<T>> const& A, std::vector<T> const& b) {
725  std::vector<std::vector<T>> B;
726  for (size_t i = 0; i < b.size(); i++) {
727  B.push_back(std::vector<T>(1, b[i]));
728  }
729  B = linsolve(A, B);
730  B[0].resize(B.size(), 0.0);
731  for (size_t i = 1; i < B.size(); i++) {
732  B[0][i] = B[i][0];
733  }
734  return B[0];
735 };
736 
737 template <class T>
738 std::vector<T> get_row(std::vector<std::vector<T>> const& in, size_t row) {
739  return in[row];
740 };
741 template <class T>
742 std::vector<T> get_col(std::vector<std::vector<T>> const& in, size_t col) {
743  std::size_t sizeX = in.size();
744  if (sizeX < 1) throw ValueError(format("You have to provide values, a vector length of %d is not valid. ", sizeX));
745  size_t sizeY = in[0].size();
746  if (sizeY < 1) throw ValueError(format("You have to provide values, a vector length of %d is not valid. ", sizeY));
747  std::vector<T> out;
748  for (std::size_t i = 0; i < sizeX; i++) {
749  sizeY = in[i].size();
750  if (sizeY - 1 < col)
751  throw ValueError(format("Your matrix does not have enough entries in row %d, last index %d is less than %d. ", i, sizeY - 1, col));
752  out.push_back(in[i][col]);
753  }
754  return out;
755 };
756 
757 template <class T>
758 std::vector<std::vector<T>> make_squared(std::vector<std::vector<T>> const& in) {
759  std::size_t cols = max_cols(in);
760  std::size_t rows = num_rows(in);
761  std::size_t maxVal = 0;
762  std::vector<std::vector<T>> out;
763  std::vector<T> tmp;
764 
765  if (cols > rows) {
766  maxVal = cols;
767  } else {
768  maxVal = rows;
769  }
770  out.clear();
771  for (std::size_t i = 0; i < in.size(); i++) {
772  tmp.clear();
773  for (std::size_t j = 0; j < in[i].size(); j++) {
774  tmp.push_back(in[i][j]);
775  }
776  while (maxVal > tmp.size()) {
777  tmp.push_back(0.0);
778  }
779  out.push_back(tmp);
780  }
781  // Check rows
782  tmp.clear();
783  tmp.resize(maxVal, 0.0);
784  while (maxVal > out.size()) {
785  out.push_back(tmp);
786  }
787  return out;
788 };
789 
790 template <class T>
791 T multiply(std::vector<T> const& a, std::vector<T> const& b) {
792  return dot_product(a, b);
793 };
794 template <class T>
795 std::vector<T> multiply(std::vector<std::vector<T>> const& A, std::vector<T> const& b) {
796  std::vector<std::vector<T>> B;
797  for (size_t i = 0; i < b.size(); i++) {
798  B.push_back(std::vector<T>(1, b[i]));
799  }
800  B = multiply(A, B);
801  B[0].resize(B.size(), 0.0);
802  for (size_t i = 1; i < B.size(); i++) {
803  B[0][i] = B[i][0];
804  }
805  return B[0];
806 }
807 
808 template <class T>
809 std::vector<std::vector<T>> multiply(std::vector<std::vector<T>> const& A, std::vector<std::vector<T>> const& B) {
810  if (num_cols(A) != num_rows(B)) {
811  throw ValueError(format("You have to provide matrices with the same columns and rows: %d is not equal to %d. ", num_cols(A), num_rows(B)));
812  }
813  size_t rows = num_rows(A);
814  size_t cols = num_cols(B);
815  T tmp;
816  std::vector<std::vector<T>> outVec;
817  std::vector<T> tmpVec;
818  outVec.clear();
819  for (size_t i = 0; i < rows; i++) {
820  tmpVec.clear();
821  for (size_t j = 0; j < cols; j++) {
822  tmp = 0.0;
823  for (size_t k = 0; k < num_cols(A); k++) {
824  tmp += A[i][k] * B[k][j];
825  }
826  tmpVec.push_back(tmp);
827  }
828  outVec.push_back(tmpVec);
829  }
830  return outVec;
831 };
832 
833 template <class T>
834 T dot_product(std::vector<T> const& a, std::vector<T> const& b) {
835  if (a.size() == b.size()) {
836  return std::inner_product(a.begin(), a.end(), b.begin(), 0.0);
837  }
838  throw ValueError(format("You have to provide vectors with the same length: %d is not equal to %d. ", a.size(), b.size()));
839 };
840 
841 template <class T>
842 std::vector<T> cross_product(std::vector<T> const& a, std::vector<T> const& b) {
843  throw NotImplementedError("The cross product function has not been implemented, yet");
844 };
845 
846 template <class T>
847 std::vector<std::vector<T>> transpose(std::vector<std::vector<T>> const& in) {
848  size_t sizeX = in.size();
849  if (sizeX < 1) throw ValueError(format("You have to provide values, a vector length of %d is not a valid. ", sizeX));
850  size_t sizeY = in[0].size();
851  size_t sizeYOld = sizeY;
852  if (sizeY < 1) throw ValueError(format("You have to provide values, a vector length of %d is not a valid. ", sizeY));
853  std::vector<std::vector<T>> out(sizeY, std::vector<T>(sizeX));
854  for (size_t i = 0; i < sizeX; ++i) {
855  sizeY = in[i].size();
856  if (sizeY != sizeYOld) throw ValueError(format("You have to provide a rectangular matrix: %d is not equal to %d. ", sizeY, sizeYOld));
857  for (size_t j = 0; j < sizeY; ++j) {
858  out[j][i] = in[i][j];
859  }
860  }
861  return out;
862 };
863 
864 template <class T>
865 std::vector<std::vector<T>> invert(std::vector<std::vector<T>> const& in) {
866  if (!is_squared(in)) throw ValueError(format("Only square matrices can be inverted: %d is not equal to %d. ", num_rows(in), num_cols(in)));
867  std::vector<std::vector<T>> identity;
868  // Build the identity matrix
869  size_t dim = num_rows(in);
870  identity.resize(dim, std::vector<T>(dim, 0));
871  for (size_t row = 0; row < dim; row++) {
872  identity[row][row] = 1.0;
873  }
874  return linsolve(in, identity);
875 };
876 
877 inline void removeRow(Eigen::MatrixXd& matrix, unsigned int rowToRemove) {
878  unsigned int numRows = static_cast<unsigned int>(matrix.rows()) - 1;
879  unsigned int numCols = static_cast<unsigned int>(matrix.cols());
880 
881  if (rowToRemove <= numRows)
882  matrix.block(rowToRemove, 0, numRows - rowToRemove, numCols) = matrix.block(rowToRemove + 1, 0, numRows - rowToRemove, numCols);
883  else {
884  throw ValueError(format("Trying to remove row index [%d] greater than max index [%d] ", rowToRemove, numRows));
885  }
886  matrix.conservativeResize(numRows, numCols);
887 };
888 
889 inline void removeColumn(Eigen::MatrixXd& matrix, unsigned int colToRemove) {
890  unsigned int numRows = static_cast<unsigned int>(matrix.rows());
891  unsigned int numCols = static_cast<unsigned int>(matrix.cols()) - 1;
892 
893  if (colToRemove <= numCols)
894  matrix.block(0, colToRemove, numRows, numCols - colToRemove) = matrix.block(0, colToRemove + 1, numRows, numCols - colToRemove);
895  else {
896  throw ValueError(format("Trying to remove column index [%d] greater than max index [%d] ", colToRemove, numCols));
897  }
898  matrix.conservativeResize(numRows, numCols);
899 };
900 template <typename Derived>
901 inline Eigen::MatrixXd minor_matrix(const Eigen::MatrixBase<Derived>& A, std::size_t i, std::size_t j) {
902  Eigen::MatrixXd Am = A;
903  removeRow(Am, static_cast<unsigned int>(i));
904  removeColumn(Am, static_cast<unsigned int>(j));
905  return Am;
906 };
907 
908 template <typename Derived>
909 static Eigen::MatrixXd adjugate(const Eigen::MatrixBase<Derived>& A) {
910  std::size_t N = A.rows();
911  if (N == 1) {
912  Eigen::MatrixXd Aadj(1, 1);
913  Aadj << 1;
914  return Aadj;
915  }
916  Eigen::MatrixXd Aadj(N, N);
917  for (std::size_t i = 0; i < N; ++i) {
918  for (std::size_t j = 0; j < N; ++j) {
919  int negative_1_to_the_i_plus_j = ((i + j) % 2 == 0) ? 1 : -1;
920  Aadj(i, j) = negative_1_to_the_i_plus_j * minor_matrix(A, j, i).determinant();
921  }
922  }
923  return Aadj;
924 }
925 
926 template <typename Derived>
927 static Eigen::MatrixXd adjugate_derivative(const Eigen::MatrixBase<Derived>& A, const Eigen::MatrixBase<Derived>& dAdt) {
928  std::size_t N = A.rows();
929  Eigen::MatrixXd Aadj(N, N);
930  for (std::size_t i = 0; i < N; ++i) {
931  for (std::size_t j = 0; j < N; ++j) {
932  int negative_1_to_the_i_plus_j = ((i + j) % 2 == 0) ? 1 : -1;
933  Eigen::MatrixXd mm = minor_matrix(A, j, i);
934  Aadj(i, j) = negative_1_to_the_i_plus_j * (adjugate(minor_matrix(A, j, i)) * minor_matrix(dAdt, j, i)).trace();
935  }
936  }
937  return Aadj;
938 }
939 
940 }; /* namespace CoolProp */
941 #endif