CoolProp 6.8.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
22template <size_t dimcount, typename T>
23struct VectorNd
24{
25 typedef std::vector<typename VectorNd<dimcount - 1, T>::type> type;
26};
27template <typename T>
28struct VectorNd<0, T>
29{
30 typedef T type;
31};
32
33namespace CoolProp {
34
36template <class T>
37std::size_t num_rows(std::vector<T> const& in) {
38 return in.size();
39}
40template <class T>
41std::size_t num_rows(std::vector<std::vector<T>> const& in) {
42 return in.size();
43}
44
45template <class T>
46std::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};
57template <class T>
58bool 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
72template <class T>
73std::size_t num_cols(std::vector<T> const& in) {
74 return 1;
75}
76template <class T>
77std::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
96template <typename T>
97std::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}
118template <class T>
119std::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
135template <class T>
136Eigen::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}
150template <class T>
151Eigen::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}
167template <class T>
168Eigen::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
180template <class T>
181Eigen::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}
194template <class T>
195Eigen::Matrix<T, Eigen::Dynamic, 1> makeVector(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& matrix) {
196 return makeColVector(matrix);
197}
198template <class T>
199Eigen::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
217template <class T>
218void 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
235template <class T>
236void 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
315static const char* stdFmt = "%8.3f";
316
318template <class T>
319std::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};
329template <class T>
330std::string vec_to_string(const std::vector<T>& a) {
331 return vec_to_string(std::vector<double>(a.begin(), a.end()), stdFmt);
332};
334inline 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
346template <class T>
347std::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};
352template <class T>
353std::string vec_to_string(const T& a) {
354 return vec_to_string((double)a, stdFmt);
355};
356
358template <class T>
359std::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};
369template <class T>
370std::string vec_to_string(const std::vector<std::vector<T>>& A) {
371 return vec_to_string(A, stdFmt);
372};
373
375template <class T>
376std::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};
397template <class T>
398std::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/*
556Owe a debt of gratitude to http://sole.ooz.ie/en - very clear treatment of GJ
557*/
558template <typename T>
559void 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};
564template <typename T>
565void 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};
570template <typename T>
571void 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
577template <typename T>
578size_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
592template <typename T>
593std::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
718template <class T>
719std::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
723template <class T>
724std::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
737template <class T>
738std::vector<T> get_row(std::vector<std::vector<T>> const& in, size_t row) {
739 return in[row];
740};
741template <class T>
742std::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
757template <class T>
758std::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
790template <class T>
791T multiply(std::vector<T> const& a, std::vector<T> const& b) {
792 return dot_product(a, b);
793};
794template <class T>
795std::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
808template <class T>
809std::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
833template <class T>
834T 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
841template <class T>
842std::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
846template <class T>
847std::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
864template <class T>
865std::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
877inline 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
889inline 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};
900template <typename Derived>
901inline 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
908template <typename Derived>
909static 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
926template <typename Derived>
927static 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
Generated on Mon Jun 2 2025 05:21:52 for CoolProp by doxygen 1.9.4