diff --git a/include/MatrixMath.h b/include/MatrixMath.h index b6009094..a9622d37 100644 --- a/include/MatrixMath.h +++ b/include/MatrixMath.h @@ -296,7 +296,10 @@ template std::size_t max_cols (std::vector > co } return cols; }; + /// Some shortcuts and regularly needed operations +template std::size_t num_rows ( std::vector const& in){ return in.size(); } +template std::size_t num_cols ( std::vector const& in){ return 1; } template std::size_t num_rows (std::vector > const& in){ return in.size(); } template std::size_t num_cols (std::vector > const& in){ if (num_rows(in)>0) { diff --git a/src/MatrixMath.cpp b/src/MatrixMath.cpp index 4d868ad4..5636c58d 100644 --- a/src/MatrixMath.cpp +++ b/src/MatrixMath.cpp @@ -13,403 +13,29 @@ namespace CoolProp{ - - - - - - - - - - - - - - -///* -//Owe a debt of gratitude to http://sole.ooz.ie/en - very clear treatment of GJ -//*/ -//template void swap_rows(std::vector > *A, size_t row1, size_t row2) -//{ -// for (size_t col = 0; col < (*A)[0].size(); col++){ -// std::swap((*A)[row1][col],(*A)[row2][col]); -// } -//} -//template void subtract_row_multiple(std::vector > *A, size_t row, T multiple, size_t pivot_row) -//{ -// for (size_t col = 0; col < (*A)[0].size(); col++){ -// (*A)[row][col] -= multiple*(*A)[pivot_row][col]; -// } -//} -//template void divide_row_by(std::vector > *A, size_t row, T value) -//{ -// for (size_t col = 0; col < (*A)[0].size(); col++){ -// (*A)[row][col] /= value; -// } -//} -// -//template size_t get_pivot_row(std::vector > *A, size_t col) -//{ -// int index = col; -// T max = 0, val; -// -// for (size_t row = col; row < (*A).size(); row++) -// { -// val = (*A)[row][col]; -// if (fabs(val) > max) -// { -// max = fabs(val); -// index = row; -// } -// } -// return index; -//} -// -// -//template std::vector > linsolve_Gauss_Jordan(std::vector > const& A, std::vector > const& B) { -// std::vector > AB; -// std::vector > X; -// size_t pivot_row; -// T pivot_element; -// -// size_t NrowA = num_rows(A); -// size_t NrowB = num_rows(B); -// size_t NcolA = num_cols(A); -// size_t NcolB = num_cols(B); -// -// if (NrowA!=NrowB) throw ValueError(format("You have to provide matrices with the same number of rows: %d is not %d. ",NrowA,NrowB)); -// -// AB.resize(NrowA, std::vector(NcolA+NcolB, 0)); -// X.resize(NrowA, std::vector(NcolB, 0)); -// -// // Build the augmented matrix -// for (size_t row = 0; row < NrowA; row++){ -// for (size_t col = 0; col < NcolA; col++){ -// AB[row][col] = A[row][col]; -// } -// for (size_t col = NcolA; col < NcolA+NcolB; col++){ -// AB[row][col] = B[row][col-NcolA]; -// } -// } -// -// for (size_t col = 0; col < NcolA; col++){ -// // Find the pivot value -// pivot_row = get_pivot_row(&AB, col); -// -// if (fabs(AB[pivot_row][col]) < 10*DBL_EPSILON){ throw ValueError(format("Zero occurred in row %d, the matrix is singular. ",pivot_row));} -// -// if (pivot_row>=col){ -// // Swap pivot row and current row -// swap_rows(&AB, col, pivot_row); -// } -// // Get the pivot element -// pivot_element = AB[col][col]; -// // Divide the pivot row by the pivot element -// divide_row_by(&AB,col,pivot_element); -// -// if (col < NrowA-1) -// { -// // All the rest of the rows, subtract the value of the [r][c] combination -// for (size_t row = col + 1; row < NrowA; row++) -// { -// subtract_row_multiple(&AB,row,AB[row][col],col); -// } -// } -// } -// for (int col = NcolA - 1; col > 0; col--) -// { -// for (int row = col - 1; row >=0; row--) -// { -// subtract_row_multiple(&AB,row,AB[row][col],col); -// } -// } -// // Set the output value -// for (size_t row = 0; row < NrowA; row++){ -// for (size_t col = 0; col < NcolB; col++){ -// X[row][col] = AB[row][NcolA+col]; -// } -// } -// return X; -//} -// -// -////std::vector > linsolve_Gauss_Jordan_reimpl(std::vector > const& A, std::vector > const& B) { -//// std::vector > AB; -//// std::vector > X; -//// size_t pivot_row; -//// double pivot_element; -//// double tmp_element; -//// -//// size_t NrowA = num_rows(A); -//// size_t NrowB = num_rows(B); -//// size_t NcolA = num_cols(A); -//// size_t NcolB = num_cols(B); -//// -//// if (NrowA!=NrowB) throw ValueError(format("You have to provide matrices with the same number of rows: %d is not %d. ",NrowA,NrowB)); -//// -//// AB.resize(NrowA, std::vector(NcolA+NcolB, 0)); -//// X.resize(NrowA, std::vector(NcolB, 0)); -//// -//// // Build the augmented matrix -//// for (size_t row = 0; row < NrowA; row++){ -//// for (size_t col = 0; col < NcolA; col++){ -//// AB[row][col] = A[row][col]; -//// } -//// for (size_t col = NcolA; col < NcolA+NcolB; col++){ -//// AB[row][col] = B[row][col-NcolA]; -//// } -//// } -//// -//// for (size_t col = 0; col < NcolA; col++){ -//// // Find the pivot row -//// pivot_row = 0; -//// pivot_element = 0.0; -//// for (size_t row = col; row < NrowA; row++){ -//// tmp_element = fabs(AB[row][col]); -//// if (tmp_element>pivot_element) { -//// pivot_element = tmp_element; -//// pivot_row = row; -//// } -//// } -//// // Check for errors -//// if (AB[pivot_row][col]<1./_HUGE) throw ValueError(format("Zero occurred in row %d, the matrix is singular. ",pivot_row)); -//// // Swap the rows -//// if (pivot_row>col) { -//// for (size_t colInt = 0; colInt < NcolA; colInt++){ -//// std::swap(AB[pivot_row][colInt],AB[pivot_row][colInt]); -//// } -//// } -//// // Process the entries below current element -//// for (size_t row = col; row < NrowA; row++){ -//// // Entries to the right of current element (until end of A) -//// for (size_t colInt = col+1; colInt < NcolA; colInt++){ -//// // All entries in augmented matrix -//// for (size_t colFull = col; colFull < NcolA+NcolB; colFull++){ -//// AB[colInt][colFull] -= AB[col][colFull] * AB[colInt][col] / AB[col][col]; -//// } -//// AB[colInt][col] = 0.0; -//// } -//// } -//// } -//// return AB; -////} -// -// -// -// -// -// -//template std::vector > linsolve(std::vector > const& A, std::vector > const& B){ -// return linsolve_Gauss_Jordan(A, B); -//} -// -//template std::vector linsolve(std::vector > const& A, std::vector const& b){ -// std::vector > B; -// for (size_t i = 0; i < b.size(); i++){ -// B.push_back(std::vector(1,b[i])); -// } -// B = linsolve(A, B); -// B[0].resize(B.size(),0.0); -// for (size_t i = 1; i < B.size(); i++){ -// B[0][i] = B[i][0]; -// } -// return B[0]; -//} -// -// -///// Some shortcuts and regularly needed operations -//template std::size_t num_rows (std::vector > const& in){ return in.size(); } -//template std::size_t num_cols (std::vector > const& in){ -// if (num_rows(in)>0) { -// if (is_squared(in)) { -// return in[0].size(); -// } else { -// return max_cols(in); -// } -// } else { -// return 0; -// } -//} -//template std::size_t max_cols (std::vector > const& in){ -// std::size_t cols = 0; -// std::size_t col = 0; -// for (std::size_t i = 0; i < in.size(); i++) { -// col = in[i].size(); -// if (cols std::vector get_row(std::vector< std::vector > const& in, size_t row) { return in[row]; } -//template std::vector get_col(std::vector< std::vector > const& in, size_t col) { -// std::size_t sizeX = in.size(); -// if (sizeX<1) throw ValueError(format("You have to provide values, a vector length of %d is not valid. ",sizeX)); -// size_t sizeY = in[0].size(); -// if (sizeY<1) throw ValueError(format("You have to provide values, a vector length of %d is not valid. ",sizeY)); -// std::vector out; -// for (std::size_t i = 0; i < sizeX; i++) { -// sizeY = in[i].size(); -// if (sizeY-1 bool is_squared(std::vector > const& in){ -// std::size_t cols = max_cols(in); -// if (cols!=num_rows(in)) { return false;} -// else { -// for (std::size_t i = 0; i < in.size(); i++) { -// if (cols!=in[i].size()) {return false; } -// } -// } -// return true; -//} -//template std::vector > make_squared(std::vector > const& in){ -// std::size_t cols = max_cols(in); -// std::size_t rows = num_rows(in); -// std::size_t maxVal = 0; -// std::vector > out; -// std::vector tmp; -// -// if (cols>rows) {maxVal = cols; } -// else {maxVal = rows; } -// out.clear(); -// for (std::size_t i = 0; i < in.size(); i++) { -// tmp.clear(); -// for (std::size_t j = 0; j < in[i].size(); j++) { -// tmp.push_back(in[i][j]); -// } -// while (maxVal>tmp.size()) { -// tmp.push_back(0.0); -// } -// out.push_back(tmp); -// } -// // Check rows -// tmp.clear(); -// tmp.resize(maxVal,0.0); -// while (maxVal>out.size()) { -// out.push_back(tmp); -// } -// return out; -//} -// -//template T multiply( std::vector const& a, std::vector const& b){ -// return dot_product(a,b); -// -//} -//template std::vector multiply(std::vector > const& A, std::vector const& b){ -// std::vector > B; -// for (size_t i = 0; i < b.size(); i++){ -// B.push_back(std::vector(1,b[i])); -// } -// B = multiply(A, B); -// B[0].resize(B.size(),0.0); -// for (size_t i = 1; i < B.size(); i++){ -// B[0][i] = B[i][0]; -// } -// return B[0]; -//} -// -//template std::vector > multiply(std::vector > const& A, std::vector > const& B){ -// if (num_cols(A) != num_rows(B)){ -// 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))); -// } -// size_t rows = num_rows(A); -// size_t cols = num_cols(B); -// T tmp; -// std::vector > outVec; -// std::vector tmpVec; -// outVec.clear(); -// for (size_t i = 0; i < rows; i++){ -// tmpVec.clear(); -// for (size_t j = 0; j < cols; j++){ -// tmp = 0.0; -// for (size_t k = 0; k < num_cols(A); k++){ -// tmp += A[i][k] * B[k][j]; -// } -// tmpVec.push_back(tmp); -// } -// outVec.push_back(tmpVec); -// } -// return outVec; -//} -// -//template T dot_product(std::vector const& a, std::vector const& b){ -// if (a.size()==b.size()){ -// return std::inner_product(a.begin(), a.end(), b.begin(), 0.0); -// } -// throw ValueError(format("You have to provide vectors with the same length: %d is not equal to %d. ",a.size(),b.size())); -//} -// -//template std::vector cross_product(std::vector const& a, std::vector const& b){ -// throw NotImplementedError("The cross product function has not been implemented, yet"); -//} -// -//template std::vector< std::vector > transpose(std::vector > const& in){ -// size_t sizeX = in.size(); -// if (sizeX<1) throw ValueError(format("You have to provide values, a vector length of %d is not a valid. ",sizeX)); -// size_t sizeY = in[0].size(); -// size_t sizeYOld = sizeY; -// if (sizeY<1) throw ValueError(format("You have to provide values, a vector length of %d is not a valid. ",sizeY)); -// std::vector< std::vector > out(sizeY,std::vector(sizeX)); -// for (size_t i = 0; i < sizeX; ++i){ -// sizeY = in[i].size(); -// if (sizeY!=sizeYOld) throw ValueError(format("You have to provide a rectangular matrix: %d is not equal to %d. ",sizeY,sizeYOld)); -// for (size_t j = 0; j < sizeY; ++j){ -// out[j][i] = in[i][j]; -// } -// } -// return out; -//} -// -//template std::vector< std::vector > invert(std::vector > const& in){ -// 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))); -// std::vector > identity; -// // Build the identity matrix -// size_t dim = num_rows(in); -// identity.resize(dim, std::vector(dim, 0)); -// for (size_t row = 0; row < dim; row++){ -// identity[row][row] = 1.0; -// } -// return linsolve(in,identity); -//} -// -//template std::string vec_to_string( T const& a){ -// std::stringstream out; -// out << format("[ %7.3f ]",a); -// return out.str(); -//} -// -//template std::string vec_to_string( std::vector const& a) { -// return vec_to_string(a,"%7.3g"); -//} -//template std::string vec_to_string( std::vector const& a, const char *fmt) { -// if (a.size()<1) { -// return std::string(""); -// } else { -// std::stringstream out; -// out << format("[ "); -// out << format(fmt,a[0]); -// for (size_t j = 1; j < a.size(); j++) { -// out << ", "; -// out << format(fmt,a[j]); -// } -// out << " ]"; -// return out.str(); -// } -//} -// -//template std::string vec_to_string(std::vector > const& A) { -// return vec_to_string(A, "%7.3g"); -//} -// -//template std::string vec_to_string(std::vector > const& A, const char *fmt) { -// std::stringstream out; -// for (size_t j = 0; j < A.size(); j++) { -// out << vec_to_string(A[j], fmt); -// } -// return out.str(); -//} - - }; /* namespace CoolProp */ + + + + +#ifdef ENABLE_CATCH +#include +#include "catch.hpp" + +TEST_CASE("Internal consistency checks and example use cases for MatrixMath.h","[MatrixMath]") +{ + /// Test case for "SylthermXLT" by "Dow Chemicals" + std::vector cHeat; + cHeat.clear(); + cHeat.push_back(+1.1562261074E+03); + cHeat.push_back(+2.0994549103E+00); + cHeat.push_back(+7.7175381057E-07); + cHeat.push_back(-3.7008444051E-20); + + SECTION("Eigen::Vector from std::vector") { + Eigen::Matrix matrix; + CoolProp::convert(cHeat, matrix); + } +} + +#endif /* CATCH_ENABLED */ diff --git a/src/Tests/eigenTest.cpp b/src/Tests/eigenTest.cpp index 66803074..ec211b7d 100644 --- a/src/Tests/eigenTest.cpp +++ b/src/Tests/eigenTest.cpp @@ -32,7 +32,7 @@ std::cout << "Evaluation of the polynomial at " << input << std::endl; std::cout << eval << std::endl; double vec0 = 0.1; -std::vector vec1(2,0.2); +std::vector vec1(2,0.44); std::vector< std::vector > vec2; vec2.push_back(std::vector(2,0.2)); vec2.push_back(std::vector(2,0.3)); @@ -56,6 +56,12 @@ CoolProp::convert(vec2, mat2); CoolProp::convert(mat2, vec); std::cout << CoolProp::vec_to_string(vec) << std::endl; +Eigen::Matrix mat1; +CoolProp::convert(vec1, mat1); +std::vector vec3; +CoolProp::convert(mat1, vec); +std::cout << CoolProp::vec_to_string(vec) << std::endl; + //std::vector< std::vector > vec(vec2); //CoolProp::convert(mat,vec);