diff --git a/BuildIssues.txt b/BuildIssues.txt index e995ea5105..a05b5301c3 100644 --- a/BuildIssues.txt +++ b/BuildIssues.txt @@ -118,34 +118,4 @@ your `~/.spack/packages.yaml` file: It appears that spack will otherwise ignore any paraview+qt version and attempt to install a paraview~qt version instead. ---------------------------- -Building on Darwin (Mac-OS) ---------------------------- - -Support for Darwin is incomplete, but has been provisioned for. - -The following are typical (as of yet) unresolved issues. - -* Scotch, ptscotch: - - The librt linkage is required for Linux, but not for Darwin. - - Current resolution: - Edit or patch - src/parallel/decompose/ptscotchDecomp/Make/options - src/parallel/decompose/scotchDecomp/Make/options - - to remove the '-lrt' library - -* CGAL: - - ThirdParty CGAL will normally need to be compiled without mpfr/gmp. - This should be done manually prior to building OpenFOAM or other - ThirdParty. Eg, - - cd $WM_THIRD_PARTY_DIR - ./makeCGAL gmp-none mpfr-none - - The erroneous references to gmp/mpfr library can be directly removed - from the wmake/rules/General/CGAL, but it is more advisable to - override them instead in the wmake/rules/darwin64Clang/CGAL file. - -- diff --git a/applications/test/00-dummy/dummy/dummyLib.C b/applications/test/00-dummy/dummy/dummyLib.C index 383c431de5..9963595c01 100644 --- a/applications/test/00-dummy/dummy/dummyLib.C +++ b/applications/test/00-dummy/dummy/dummyLib.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2018 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2018-2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -30,6 +30,9 @@ License #if defined WM_SP # define PRECISION "SP" # define SCALAR_SIZE (8*sizeof(float)) +#elif defined(WM_SPDP) +# define PRECISION "SPDP" +# define SCALAR_SIZE (8*sizeof(float)) #elif defined WM_DP # define PRECISION "DP" # define SCALAR_SIZE (8*sizeof(double)) @@ -38,9 +41,23 @@ License # define SCALAR_SIZE (8*sizeof(long double)) #endif +// Test additional exported symbols +#ifdef _WIN32 + #define defineWindowsLibEntryPoint(libName) \ + extern "C" void lib_##libName##_entry_point() {} +#else + #define defineWindowsLibEntryPoint(libName) /* Nothing */ +#endif + // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // +// The 'extern C' export is independent of namespace +namespace Foam +{ + defineWindowsLibEntryPoint(dummyLib); +} + const std::string Foam::Detail::dummyLib::arch(WM_ARCH); const std::string Foam::Detail::dummyLib::compiler(WM_COMPILER); diff --git a/applications/test/Matrix/Test-Matrix.C b/applications/test/Matrix/Test-Matrix.C index 7c9cb61195..dbccf225fa 100644 --- a/applications/test/Matrix/Test-Matrix.C +++ b/applications/test/Matrix/Test-Matrix.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | + \\ / A nd | Copyright (C) 2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- | Copyright (C) 2011-2016 OpenFOAM Foundation @@ -25,162 +25,946 @@ License \*---------------------------------------------------------------------------*/ -#include "scalarMatrices.H" +#include "MatrixTools.H" #include "LUscalarMatrix.H" #include "LLTMatrix.H" -#include "QRMatrix.H" -#include "vector.H" -#include "tensor.H" -#include "IFstream.H" +#include "Random.H" +#include "SortList.H" +#include using namespace Foam; +using namespace Foam::MatrixTools; -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -// Main program: +#define isEqual MatrixTools::equal + +void horizontalLine() +{ + Info<< "+---------+---------+---------+---------+---------+" << nl; +} + + +// Copy values into matrix +template +void assignMatrix +( + Matrix& mat, + std::initializer_list::cmptType> list +) +{ + const label nargs = list.size(); + + if (nargs != mat.size()) + { + FatalErrorInFunction + << "Mismatch in matrix dimension (" + << mat.m() << ", " + << mat.n() << ") and number of args (" << nargs << ')' << nl + << exit(FatalError); + } + + std::copy(list.begin(), list.end(), mat.begin()); +} + + +template +MatrixType makeRandomMatrix +( + const labelPair& dims, + Random& rndGen +) +{ + MatrixType mat(dims); + + std::generate + ( + mat.begin(), + mat.end(), + [&]{return rndGen.GaussNormal();} + ); + + return mat; +} + + +template +Ostream& operator<<(Ostream& os, const ConstMatrixBlock& mat) +{ + return MatrixTools::printMatrix(os, mat); +} + + +template +Ostream& operator<<(Ostream& os, const MatrixBlock& mat) +{ + return MatrixTools::printMatrix(os, mat); +} + + +// * * * * * * * * * * * * * * * Main Program * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { - SquareMatrix hmm(3); + typedef SquareMatrix SMatrix; + typedef SquareMatrix SCMatrix; + typedef RectangularMatrix RMatrix; + typedef RectangularMatrix RCMatrix; - hmm(0, 0) = -3.0; - hmm(0, 1) = 10.0; - hmm(0, 2) = -4.0; - hmm(1, 0) = 2.0; - hmm(1, 1) = 3.0; - hmm(1, 2) = 10.0; - hmm(2, 0) = 2.0; - hmm(2, 1) = 6.0; - hmm(2, 2) = 1.0; - - //Info<< hmm << endl << hmm - 2.0*(-hmm) << endl; - Info<< max(hmm) << endl; - Info<< min(hmm) << endl; - - SquareMatrix hmm2(3, I); - - hmm = hmm2; - - Info<< hmm << endl; - - //SquareMatrix hmm3(Sin); - - //Info<< hmm3 << endl; - - SquareMatrix hmm4; - - hmm4 = hmm2; - - Info<< hmm4 << endl; - - SquareMatrix hmm5; - - hmm4 = hmm5; - Info<< hmm5 << endl; + Random rndGen(1234); + #if 1 { - RectangularMatrix rm1(5, 6, 3.1); - rm1(0, 1) = 4.5; - RectangularMatrix rm1b(rm1.block(2, 2, 0, 0)); - Info<< "rm1b = " << rm1b << endl; + horizontalLine(); + + Info<< "## Matrix Constructors & Assignments:" << nl << nl; + + Info<< "# SquareMatrix examples:" << nl; + { + const label mRows = 2; + const label nCols = 2; + Matrix MS0(mRows, nCols); + Matrix MS1(mRows, nCols, Zero); + Matrix MS2(mRows, nCols, 17.44); + Matrix MS3(MS2); + + SMatrix S0(mRows); + SMatrix S1(mRows, Zero); + SMatrix S2(mRows, I); + SMatrix S3(mRows, 17.44); + SMatrix S4(2, Zero); + SMatrix S5(labelPair{mRows, nCols}); + SMatrix S6(mRows, nCols, Zero); + SMatrix S7({mRows, nCols}, Zero); + SMatrix S8({mRows, nCols}, 17.44); + assignMatrix + ( + S8, + { + 1, 2, // 0th row + 3, 4 // 1st row and so on + } + ); + + S1 = Zero; + S2 = 13.13; + S3 = I; + } + + Info<< "# RectangularMatrix examples:" << nl; + { + const label mRows = 2; + const label nCols = 3; + Matrix MR0(mRows, nCols); + Matrix MR1(mRows, nCols, Zero); + Matrix MR2(mRows, nCols, 17.44); + + RMatrix R0(mRows, nCols); + RMatrix R1(labelPair{mRows, nCols}); + RMatrix R2(mRows, nCols, Zero); + RMatrix R3({mRows, nCols}, Zero); + RMatrix R4(mRows, nCols, -17.44); + RMatrix R5({mRows, nCols}, -17.44); + RMatrix R6(3, 4, Zero); + RMatrix R7({3, 4}, Zero); + assignMatrix + ( + R7, + { + 1, 2, 3, 4, + 5, 6, 7, 8.8, + 9, 10, 11, 12 + } + ); + R0 = Zero; + R2 = -13.13; + } + + Info<< "# SquareMatrix examples:" << nl; + { + const label mRows = 2; + const label nCols = 2; + Matrix MS0(mRows, nCols); + Matrix MS1(mRows, nCols, Zero); + Matrix MS2(mRows, nCols, complex(17.44, 44.17)); + Matrix MS3(MS2); + + SCMatrix S0(mRows); + SCMatrix S1(mRows, Zero); + SCMatrix S2(mRows, I); + SCMatrix S3(mRows, complex(17.44,0)); + SCMatrix S4(2, Zero); + SCMatrix S5(labelPair{mRows, nCols}); + SCMatrix S6(mRows, nCols, Zero); + SCMatrix S7({mRows, nCols}, Zero); + SCMatrix S8({mRows, nCols}, complex(17.44,0)); + assignMatrix + ( + S8, + { + complex(1,0), complex(2,2), + complex(4,2), complex(5,0) + } + ); + + S1 = Zero; + S2 = complex(13.13,0); + S3 = I; + } + + Info<< nl; + horizontalLine(); } + #endif + + #if 1 { - scalarSymmetricSquareMatrix symmMatrix(3, Zero); + horizontalLine(); - symmMatrix(0, 0) = 4; - symmMatrix(1, 0) = 12; - symmMatrix(1, 1) = 37; - symmMatrix(2, 0) = -16; - symmMatrix(2, 1) = -43; - symmMatrix(2, 2) = 98; + Info<< "## Access Functions:" << nl; - Info<< "Symmetric Square Matrix = " << symmMatrix << endl; + SMatrix S(3, Zero); + assignMatrix + ( + S, + { + -3, 11, -4, + 2, 3, 10, + 2, 6, 1 + } + ); + S(0,1) = 10; - Info<< "Inverse = " << inv(symmMatrix) << endl; - Info<< "Determinant = " << det(symmMatrix) << endl; + Info<< "# SquareMatrix example:" << nl; + printMatrix(Info, S) << nl; - scalarSymmetricSquareMatrix symmMatrix2(symmMatrix); - LUDecompose(symmMatrix2); + const label mRows = S.m(); + const label nCols = S.n(); + const label size = S.size(); + const labelPair sizes = S.sizes(); + const bool isEmpty = S.empty(); + const long constPointer = long(S.cdata()); + long pointer = long(S.data()); - Info<< "Inverse = " << invDecomposed(symmMatrix2) << endl; - Info<< "Determinant = " << detDecomposed(symmMatrix2) << endl; + Info + << "Number of rows =" << tab << mRows << nl + << "Number of columns = " << tab << nCols << nl + << "Number of elements = " << tab << size << nl + << "Number of rows/columns = " << tab << sizes << nl + << "Matrix is empty = " << tab << isEmpty << nl + << "Constant pointer = " << tab << constPointer << nl + << "Pointer = " << tab << pointer << nl + << nl; + + horizontalLine(); + } + #endif + + + #if 1 + { + horizontalLine(); + + Info<< "## Block Access Functions:" << nl; + + RMatrix R(4, 3, Zero); + assignMatrix + ( + R, + { + -3, 11, -4, + 2, 3, 10, + 2, 6, 1, + 1, 2, 3 + } + ); + + Info<< "# RectangularMatrix example:" << nl; + printMatrix(Info, R) << nl; + + // Indices begin at 0 + const label columnI = 1; + const label rowI = 2; + const label colStartI = 1; + const label rowStartI = 1; + const label szRows = 2; + const label szCols = 2; + + Info + << "column: " << R.subColumn(columnI) << nl + << "sub-column: " << R.subColumn(columnI, rowStartI) << nl + << "sub-sub-column: " << R.subColumn(columnI, rowStartI, szRows) << nl + << "row: " << R.subRow(rowI) << nl + << "sub-row: " << R.subRow(rowI, colStartI) << nl + << "sub-sub-row: " << R.subRow(rowI, colStartI, szCols) << nl + << "sub-block: " << R.subMatrix(rowStartI, colStartI) << nl + << "sub-block with row size: " << R.subMatrix(rowStartI, colStartI, szRows) << nl + << "sub-block with row|col size: " << R.subMatrix(rowStartI, colStartI, szRows, szCols) << nl; + + horizontalLine(); + } + #endif + + + #if 1 + { + horizontalLine(); + + Info<< "## Edit Functions:" << nl; + + RMatrix R(2, 2, Zero); + assignMatrix + ( + R, + { + 1, 2, + 3, 4 + } + ); + + Info<< "# Before Matrix.clear():" << nl << R << nl; + R.clear(); + Info<< "# After Matrix.clear():" << nl << R << nl; + + Info<< "# Before Matrix.resize(m, n):" << nl << R; + R.resize(3, 2); + Info<< "# After Matrix.resize(m, n):" << nl << R << nl << nl; + + // Steal matrix contents + Info<< "# Before Matrix.release():" << nl << R << nl; + List newOwner(R.release()); + Info<< "# After Matrix.release():" << nl << R + << "& corresponding List content:" << nl << newOwner << nl << nl; + + R.resize(3, 2); + assignMatrix + ( + R, + { + 1, 2, + 5, 6, + 9e-19, 1e-17 + } + ); + + Info<< "# Before Matrix.round():" << nl << R << nl; + R.round(); + Info<< "# After Matrix.round():" << nl << R << nl << nl; + + Info<< "# Before Matrix.T() (Transpose):" << nl << R << nl; + RMatrix RT = R.T(); + Info<< "# After Matrix.T():" << nl << RT << nl << nl; + + RCMatrix RC(3, 2, Zero); + assignMatrix + ( + RC, + { + complex(1, 0), complex(2,5), + complex(4, 3), complex(1,-1), + complex(1, 2), complex(2,9) + } + ); + + Info<< "# Before Matrix.T() (Hermitian transpose):" << nl << RC << nl; + RCMatrix RCH = RC.T(); + Info<< "# After Matrix.T():" << nl << RCH << nl << nl; + + Info<< "# Diagonal and trace:" << nl; + RMatrix R1(5, 7, 3.4); + { + List diagR = R1.diag(); + scalar traceR = R1.trace(); + Info<< "RectangularMatrix:" + << nl << R1 << nl + << "Major diagonal of RectangularMatrix:" + << nl << diagR << nl + << "Trace of RectangularMatrix:" + << nl << traceR << nl << nl; + } + + Info<< "# List assignment to the main diagonal of Matrix:" << nl; + { + SMatrix S1(5, I); + Info<< "Before List assignment:" << nl; + printMatrix(Info, S1); + List lst(5, 2.0); + S1.diag(lst); + Info<< "After List assignment:" << nl; + printMatrix(Info, S1); + } + + Info<< "# Diagonal sort of Matrix:" << nl; + { + auto S1 + ( + makeRandomMatrix + ( + {3, 3}, rndGen + ) + ); + + List diag = S1.diag(); + SortList incrDiag(diag); + incrDiag.sort(); + + Info<< "Diagonal list:" << diag << nl + << "Ascending diagonal list: " << incrDiag << nl << nl; + } + + horizontalLine(); + } + #endif + + + #if 1 + { + horizontalLine(); + + Info<< "## Transpose Verifications:" << nl; + + const label numberOfTests = 10; + + for (label i = 0; i < numberOfTests; ++i) + { + const label mRows = rndGen.position(1,10); + const label nCols = rndGen.position(1,10); + + Info << "# SquareMatrix example:" << nl; + { + Info<< "A(mxm) where m = " << mRows << nl; + auto A(makeRandomMatrix({mRows, mRows}, rndGen)); + auto B(makeRandomMatrix({mRows, mRows}, rndGen)); + + Info<< nl << "# (A.T()).T() = A:" << nl; + isEqual((A.T()).T(), A, true); + + Info<< nl << "# (A + B).T() = A.T() + B.T():" << nl; + isEqual((A + B).T(), A.T() + B.T(), true); + + Info<< nl << "# (A*B).T() = B.T()*A.T():" << nl; + isEqual((A*B).T(), B.T()*A.T(), true); + Info<< nl << nl; + } + + Info << "# RectangularMatrix example:" << nl; + { + Info<< "A(mxn) where" + << " m = " << mRows << "," + << " n = " << nCols << nl; + auto A(makeRandomMatrix({mRows, nCols}, rndGen)); + auto B(makeRandomMatrix({mRows, nCols}, rndGen)); + auto C(makeRandomMatrix({nCols, mRows}, rndGen)); + + Info<< nl << "# (A.T()).T() = A:" << nl; + isEqual((A.T()).T(), A, true); + + Info<< nl << "# (A + B).T() = A.T() + B.T():" << nl; + isEqual((A + B).T(), A.T() + B.T(), true); + + Info<< nl << "# (A*C).T() = C.T()*C.T():" << nl; + isEqual((A*C).T(), C.T()*A.T(), true); + Info<< nl << nl; + } + + Info << "# SquareMatrix example:" << nl; + { + Info<< "A(mxm) where m = " << mRows << nl; + SCMatrix A(mRows); + SCMatrix B(mRows); + + for (auto& val : A) + { + val.Re() = rndGen.GaussNormal(); + val.Im() = rndGen.GaussNormal(); + } + for (auto& val : B) + { + val.Re() = rndGen.GaussNormal(); + val.Im() = rndGen.GaussNormal(); + } + + Info<< nl << "# (A.T()).T() = A:" << nl; + isEqual((A.T()).T(), A, true); + + Info<< nl << "# (A + B).T() = A.T() + B.T():" << nl; + isEqual((A + B).T(), A.T() + B.T(), true); + + Info<< nl << "# (A*B).T() = B.T()*A.T():" << nl; + isEqual((A*B).T(), B.T()*A.T(), true); + Info<< nl << nl; + } + } + + horizontalLine(); + } + #endif + + + #if 1 + { + horizontalLine(); + + Info<< "## Scalar Arithmetic Operations:" << nl; + + if (true) + { + SMatrix S1(3, Zero); + assignMatrix + ( + S1, + { + 4.1, 12.5, -16.3, + -192.3, -9.1, -3.0, + 1.0, 5.02, -4.4 + } + ); + SMatrix S2 = S1; + SMatrix S3(3, Zero); + + Info<< "SquareMatrix S1 = " << S1 << nl; + Info<< "SquareMatrix S2 = " << S2 << nl; + Info<< "SquareMatrix S3 = " << S3 << nl; + + S3 = S1*S2; + Info<< "S1*S2 = " << S3 << nl; + + S3 = S1 - S2; + Info<< "S1 - S2 = " << S3 << nl; + + S3 = S1 + S2; + Info<< "S1 + S2 = " << S3 << nl; + +// S3 *= S1; // Not Implemented + + S3 -= S1; + Info<< "S3 -= S1; S3 = " << S3 << nl; + + S3 += S1; + Info<< "S3 += S1; S3 = " << S3 << nl; + + + Info<< nl << "# Scalar broadcasting:" << nl; + + S3 *= 5.0; + Info<< "S3 *= 5.0; S3 = " << S3 << nl; + + S3 /= 5.0; + Info<< "S3 /= 5.0; S3 = " << S3 << nl; + + S3 -= 1.0; + Info<< "S3 -= 1.0; S3 = " << S3 << nl; + + S3 += 1.0; + Info<< "S3 += 1.0; S3 = " << S3 << nl; + + + Info<< nl << "# Operations between different matrix types:" << nl; + RMatrix R1 = S1; + Info<< "RectangularMatrix R1 = " << R1 << nl; + + // RectangularMatrix*SquareMatrix returns RectangularMatrix + // S3 = S3*R1; // Not implemented + R1 = S3*R1; + Info<< "R1 = S3*R1; R1 = " << R1 << nl; + + S3 = S3 - R1; + Info<< "S3 = S3 - R1; S3 = " << S3 << nl; + + S3 = S3 + R1; + Info<< "S3 = S3 + R1; S3 = " << S3 << nl; + + R1 = S3 + R1; + Info<< "R1 = S3 + R1; R1 = " << R1 << nl; + + + Info<< nl << "# Extrema operations:" << nl; + + Info<< "min(S3) = " << min(S3) << nl + << "max(S3) = " << max(S3) << nl + << "minMax(S3) = " << minMax(S3) << nl; + } + + if (true) + { + Info<< nl << "# Operations with ListTypes:" << nl; + + SMatrix A(3, Zero); + assignMatrix + ( + A, + { + 4.1, 12.5, -16.3, + -192.3, -9.1, -3.0, + 1.0, 5.02, -4.4 + } + ); + const List lst({1, 2, 3}); + RMatrix colVector(3, 1, Zero); + assignMatrix + ( + colVector, + {1, 2, 3} + ); + RMatrix rowVector(1, 3, Zero); + assignMatrix + ( + rowVector, + {1, 2, 3} + ); + + Info<< "A = " << A << nl; + Info<< "colVector = " << colVector << nl; + Info<< "rowVector = " << rowVector << nl; + + const Field field1(A.Amul(lst)); + const Field field2(A*lst); + const Field field3(A.Tmul(lst)); + const Field field4(lst*A); + + Info + << "Field1 = A*lst = A.Amul(lst):" << nl << field1 << nl + << "Field2 = A*lst:" << nl << field2 << nl + << "A*colVector:" << nl << A*colVector << nl + << "Field3 = lst*A = A.Tmul(lst):" << nl << field3 << nl + << "Field4 = lst*A:" << nl << field4 << nl + << "rowVector*A:" << nl << rowVector*A << nl + << nl; + } + + if (true) + { + Info<< nl << "# Implicit inner/outer products:" << nl; + + const scalar r = 2.0; + RMatrix A(3, 2, Zero); + assignMatrix + ( + A, + { + 1, 2, + 3, 4, + 5, 6 + } + ); + const RMatrix B(A); + const RMatrix C(B); + + Info<< nl << "# Inner product:" << nl; + { + Info<< "- Explicit vs Implicit => A.T()*B == A&B:" << nl; + isEqual(A.T()*B, A&B, true); + + Info<< "- Commutative => A&B == B&A:" << nl; + isEqual(A&B, B&A, true); + + Info<< "- Distributive => A&(B+C) == A&B + A&C:" << nl; + isEqual(A&(B + C), (A&B) + (A&C), true); + + Info<< "- Bilinear => A&(rB+C) == r(A&B) + (A&C):" << nl; + isEqual(A&(r*B + C), r*(A&B) + (A&C), true); + + Info<< "- Scalar multiplication => (rA)&(rB) == rr(A&B):" << nl; + isEqual((r*A)&(r*B), r*r*(A&B), true); + } + + Info<< nl << "# Outer product:" << nl; + { + Info<< "- Explicit vs Implicit => A*B.T() == A^B:" << nl; + isEqual(A*B.T(), A^B, true); + + Info<< "- Commutative => A^B == B^A:" << nl; + isEqual(A^B, B^A, true); + + Info<< "- Distributive => A^(B+C) == A^B + A^C:" << nl; + isEqual(A^(B + C), (A^B) + (A^C), true); + + Info<< "- Scalar multiplication => r*(A^B) == (r*A)^B:" << nl; + isEqual(r*(A^B), (r*A)^B, true); + } + } + + Info<< nl; + horizontalLine(); + } + #endif + + + #if 1 + { + horizontalLine(); + + Info<< "## Complex Arithmetic Operations:" << nl; + + if (true) + { + SCMatrix S1(3, Zero); + assignMatrix + ( + S1, + { + complex(4.1, 1.0), complex(12.5, -1), complex(-16.3,-3), + complex(-192.3, 5), complex(-9.1, 3), complex(-3.0, 1), + complex(1.0, 3), complex(5.02, 0.3), complex(-4.4, 1) + } + ); + SCMatrix S2 = S1; + SCMatrix S3(3, Zero); + + Info<< "SquareMatrix S1 = " << S1 << nl; + Info<< "SquareMatrix S2 = " << S2 << nl; + Info<< "SquareMatrix S3 = " << S3 << nl; + + S3 = S1*S2; + Info<< "S1*S2 = " << S3 << nl; + + S3 = S1 - S2; + Info<< "S1 - S2 = " << S3 << nl; + + S3 = S1 + S2; + Info<< "S1 + S2 = " << S3 << nl; + +// S3 *= S1; // Not Implemented + + S3 -= S1; + Info<< "S3 -= S1; S3 = " << S3 << nl; + + S3 += S1; + Info<< "S3 += S1; S3 = " << S3 << nl; + + + Info<< nl << "# Complex broadcasting:" << nl; + + S3 *= complex(5, 0); + Info<< "S3 *= complex(5, 0); S3 = " << S3 << nl; + + S3 /= complex(5, 0); + Info<< "S3 /= complex(5, 0); S3 = " << S3 << nl; + + S3 -= complex(1, 0); + Info<< "S3 -= complex(1, 0); S3 = " << S3 << nl; + + S3 += complex(1, 0); + Info<< "S3 += complex(1, 0); S3 = " << S3 << nl; + + + Info<< nl << "# Operations between different matrix types:" << nl; + RCMatrix R1 = S1; + Info<< "RectangularMatrix R1 = " << R1 << nl; + + R1 = S3*R1; + Info<< "R1 = S3*R1; R1 = " << R1 << nl; + + S3 = S3 - R1; + Info<< "S3 = S3 - R1; S3 = " << S3 << nl; + + S3 = S3 + R1; + Info<< "S3 = S3 + R1:" << S3 << nl; + + // Extrama operations // Not Implemented + } + + if (true) + { + Info<< nl << "# Operations with ListTypes:" << nl; + + SCMatrix A(3, Zero); + assignMatrix + ( + A, + { + complex(4.1, 1.0), complex(12.5, -1), complex(-16.3,-3), + complex(-192.3, 5), complex(-9.1, 3), complex(-3.0, 1), + complex(1.0, 3), complex(5.02, 0.3), complex(-4.4, 1) + } + ); + const List lst({complex(1,1), complex(2,2), complex(3,3)}); + RCMatrix colVector(3, 1, Zero); + assignMatrix + ( + colVector, + {complex(1,1), complex(2,2), complex(3,3)} + ); + RCMatrix rowVector(1, 3, Zero); + assignMatrix + ( + rowVector, + {complex(1,1), complex(2,2), complex(3,3)} + ); + + Info<< "A = " << A << nl; + Info<< "colVector = " << colVector << nl; + Info<< "rowVector = " << rowVector << nl; + + const Field field1(A.Amul(lst)); + const Field field2(A*lst); + const Field field3(A.Tmul(lst)); + const Field field4(lst*A); + + Info + << "Field1 = A*lst = A.Amul(lst):" << nl << field1 << nl + << "Field2 = A*lst:" << nl << field2 << nl + << "A*colVector:" << nl << A*colVector << nl + << "Field3 = lst*A = A.Tmul(lst):" << nl << field3 << nl + << "Field4 = lst*A:" << nl << field4 << nl + << "rowVector*A:" << nl << rowVector*A << nl + << nl; + } + + #if 1 + { + Info<< nl << "# Implicit inner/outer products:" << nl; + + const complex r(2.0,1.0); + RCMatrix A(3, 2, Zero); + assignMatrix + ( + A, + { + complex(1,0), complex(2,1), + complex(3,-0.3), complex(4,-1), + complex(0.5,-0.1), complex(6,0.4) + } + ); + const RCMatrix B(A); + const RCMatrix C(B); + + Info<< nl << "# Inner product:" << nl; + { + Info<< "- Explicit vs Implicit => A.T()*B == A&B:" << nl; + isEqual(A.T()*B, A&B, true); + + Info<< "- Commutative => A&B == B&A:" << nl; + isEqual(A&B, B&A, true); + + Info<< "- Distributive => A&(B+C) == A&B + A&C:" << nl; + isEqual(A&(B + C), (A&B) + (A&C), true); + + Info<< "- Bilinear => A&(rB+C) == r(A&B) + (A&C):" << nl; + isEqual(A&(r*B + C), r*(A&B) + (A&C), true); + } + + Info<< nl << "# Outer product:" << nl; + { + Info<< "- Explicit vs Implicit => A*B.T() == A^B:" << nl; + isEqual(A*B.T(), A^B, true); + + Info<< "- Commutative => A^B == B^A:" << nl; + isEqual(A^B, B^A, true); + + Info<< "- Distributive => A^(B+C) == A^B + A^C:" << nl; + isEqual(A^(B + C), (A^B) + (A^C), true); + + Info<< "- Scalar multiplication => r*(A^B) == (r*A)^B:" << nl; + isEqual(r*(A^B), (r*A)^B, true); + } + } + #endif + + Info<< nl; + horizontalLine(); + } + #endif + + + #if 1 + { + horizontalLine(); + + Info<< "## SymmetricSquareMatrix algorithms:" << nl; + + scalarSymmetricSquareMatrix symm(3, Zero); + + symm(0, 0) = 4; + symm(1, 0) = 12; + symm(1, 1) = 37; + symm(2, 0) = -16; + symm(2, 1) = -43; + symm(2, 2) = 98; + + Info<< "SymmetricSquareMatrix = " << nl << symm << nl + << "Inverse = " << nl << inv(symm) << nl + << "Determinant = " << det(symm) << nl; + + scalarSymmetricSquareMatrix symm2(symm); + LUDecompose(symm2); + Info<< "LU decomposition:" << nl + << "Inverse = " << nl << invDecomposed(symm2) << nl + << "Determinant = " << detDecomposed(symm2) << nl; scalarDiagonalMatrix rhs(3, 0); rhs[0] = 1; rhs[1] = 2; rhs[2] = 3; - LUsolve(symmMatrix, rhs); - - Info<< "Decomposition = " << symmMatrix << endl; - Info<< "Solution = " << rhs << endl; + LUsolve(symm, rhs); + Info<< "Solving linear system through LU decomposition:" << nl + << "Decomposition = " << nl << symm << nl + << "Solution = " << rhs << nl; } + #endif - scalarSquareMatrix squareMatrix(3, Zero); - - squareMatrix(0, 0) = 4; - squareMatrix(0, 1) = 12; - squareMatrix(0, 2) = -16; - squareMatrix(1, 0) = 12; - squareMatrix(1, 1) = 37; - squareMatrix(1, 2) = -43; - squareMatrix(2, 0) = -16; - squareMatrix(2, 1) = -43; - squareMatrix(2, 2) = 98; - - Info<< nl << "Square Matrix = " << squareMatrix << endl; - - const scalarField source(3, 1); - + #if 1 { + scalarSquareMatrix squareMatrix(3, Zero); + + scalarSquareMatrix S(3, Zero); + assignMatrix + ( + S, + { + 4, 12, -16, + 12, 37, -43, + -16, -43, 98 + } + ); + + const scalarField source(3, 1); + + Info<< nl << "Square Matrix = " << S << endl; + + if (true) { - scalarSquareMatrix sm(squareMatrix); - Info<< "det = " << det(sm) << endl; + { + scalarSquareMatrix S2(S); + Info<< "Determinant = " << det(S2) << nl; + } + + scalarSquareMatrix S2(S); + labelList rhs(3, 0); + label sign; + LUDecompose(S2, rhs, sign); + + Info<< "LU decomposition = " << S2 << nl + << "Pivots = " << rhs << nl + << "Sign = " << sign << nl + << "Determinant = " << detDecomposed(S2, sign) << nl; } - scalarSquareMatrix sm(squareMatrix); - labelList rhs(3, 0); - label sign; - LUDecompose(sm, rhs, sign); + if (true) + { + LUscalarMatrix LU(S); + scalarField x(LU.solve(source)); + scalarSquareMatrix inv(3); + LU.inv(inv); + Info<< "LU.solve(source) residual: " << (S*x - source) + << nl << "LU inv " << inv << nl + << "LU inv*S " << (inv*S) << nl; + } - Info<< "Decomposition = " << sm << endl; - Info<< "Pivots = " << rhs << endl; - Info<< "Sign = " << sign << endl; - Info<< "det = " << detDecomposed(sm, sign) << endl; + if (true) + { + LLTMatrix LLT(S); + scalarField x(LLT.solve(source)); + Info<< "LLT solve residual " << (S*x - source) << nl; + } + + horizontalLine(); } + #endif - { - LUscalarMatrix LU(squareMatrix); - scalarField x(LU.solve(source)); - Info<< "LU solve residual " << (squareMatrix*x - source) << endl; - scalarSquareMatrix inv(3); - LU.inv(inv); - Info<< "LU inv " << inv << endl; - Info<< "LU inv*squareMatrix " << (inv*squareMatrix) << endl; - } - - { - LLTMatrix LLT(squareMatrix); - scalarField x(LLT.solve(source)); - Info<< "LLT solve residual " << (squareMatrix*x - source) << endl; - } - - { - QRMatrix QR(squareMatrix); - scalarField x(QR.solve(source)); - - Info<< "QR solve residual " - << (squareMatrix*x - source) << endl; - - Info<< "QR inverse solve residual " - << (x - QR.inv()*source) << endl; - - Info<< "QR inv *squareMatrix " << (QR.inv()*squareMatrix) << endl; - } - - Info<< "\nEnd\n" << endl; + Info<< nl << "End" << nl; return 0; } diff --git a/applications/test/complex/Test-complex.C b/applications/test/complex/Test-complex.C index 938f63fb8e..57551c5060 100644 --- a/applications/test/complex/Test-complex.C +++ b/applications/test/complex/Test-complex.C @@ -24,11 +24,12 @@ License Application Description - Some tests for complex numbers + Tests for complex numbers \*---------------------------------------------------------------------------*/ #include "argList.H" +#include "complex.H" #include "complexFields.H" #include "ops.H" #include "ListOps.H" @@ -41,15 +42,14 @@ void print1(const complex& z) } -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -// Main program: +// * * * * * * * * * * * * * * * Main Program * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { Info<< "complex() : " << complex() << nl << "complex(zero) : " << complex(Zero) << nl - << "complex::zero : " << complex::zero << nl - << "complex::one : " << complex::one << nl + << "pTraits::zero : " << pTraits::zero << nl + << "pTraits::one : " << pTraits::one << nl << "complex(scalar) : " << complex(3.14519) << nl << nl; @@ -57,7 +57,6 @@ int main(int argc, char *argv[]) Info<< "std::complex : " << c1 << nl; Info<< "sin: " << std::sin(c1) << nl; - Info<< "complexVector::zero : " << complexVector::zero << nl << "complexVector::one : " << complexVector::one << nl << nl; @@ -98,6 +97,7 @@ int main(int argc, char *argv[]) c.Im() *= 5; } + Info<< "sumProd: " << sumProd(fld1, fld2) << nl; fld1 *= 10; Info<< "scalar multiply: " << flatOutput(fld1) << nl; @@ -113,11 +113,140 @@ int main(int argc, char *argv[]) // Info<< "operator * : " << (fld1 * fld2) << nl; // Info<< "operator / : " << (fld1 / fld2) << nl; - Info<< "operator / : " << (fld1 / 2) << nl; + // Info<< "operator / : " << (fld1 / 2) << nl; // Info<< "operator / : " << (fld1 / fld2) << nl; - Info<< "sqrt : " << sqrt(fld1) << nl; + // Info<< "sqrt : " << sqrt(fld1) << nl; // Info<< "pow(2) : " << pow(fld1, 2) << nl; + #if 1 + Info<< nl << "## Elementary complex-complex arithmetic operations:" << nl; + { + const complex a(6, 1); + complex b = a; + + Info << "# Compound assignment operations:" << nl; + + Info<< "a = " << a << ", b = " << b << nl; + + // Addition + b += a; + Info<< "b += a:" << tab << "b =" << b << nl; + + // Subtraction + b -= a; + Info<< "b -= a:" << tab << "b =" << b << nl; + + // Multiplication + b *= a; + Info<< "b *= a:" << tab << "b =" << b << nl; + + // Division + b /= a; + Info<< "b /= a:" << tab << "b =" << b << nl; + } + #endif + + + #if 1 + Info<< nl << "## Elementary complex-scalar arithmetic operations:" << nl; + { + const scalar a = 5; + complex b(6, 1); + + Info << "# Non-assignment operations:" << nl; + + Info<< "(scalar) a = " << a << ", b = " << b << nl; + + // Addition + b = a + b; + Info<< "b = a + b: " << tab << b << nl; + + b = b + a; + Info<< "b = b + a: " << tab << b << nl; + + // Subtraction + b = a - b; + Info<< "b = a - b: " << tab << b << nl; + + b = b - a; + Info<< "b = b - a: " << tab << b << nl; + + // Multiplication + b = a*b; + Info<< "b = a*b: " << tab << b << nl; + + b = b*a; + Info<< "b = b*a: " << tab << b << nl; + + // Division + b = a/b; + Info<< "b = a/b = scalar(a)/b = complex(a)/b:" << tab << b << nl; + + b = b/a; + Info<< "b = b/a: " << tab << b << nl; + + + Info << "# Compound assignment operations:" << nl; + + Info<< "(scalar) a = " << a << ", b = " << b << nl; + + // Addition: complex+scalar + b += a; + Info<< "b += a (only real part):" << tab << b << nl; + + // Subtraction: complex-scalar + b -= a; + Info<< "b -= a (only real part):" << tab << b << nl; + + // Multiplication: complex*scalar + b *= a; + Info<< "b *= a (real and imag parts):" << tab << b << nl; + + // Division: complex/scalar + b /= a; + Info<< "b /= a (real and imag parts):" << tab << b << nl; + + } + #endif + + + #if 1 + Info<< nl << "## Other mathematical expressions:" << nl; + { + const complex a(4.3, -3.14); + const complex b(0, -4.3); + const complex c(-4.3, 0); + + Info<< "a = " << a << ", b = " << b << ", c = " << c << nl; + + // Square-root + Info<< "sqrt(a) = " << Foam::sqrt(a) << ", " + << "sqrt(b) = " << Foam::sqrt(b) << ", " + << "sqrt(c) = " << Foam::sqrt(c) << nl; + + // Square + Info<< "sqr(a) = " << sqr(a) << ", " + << "sqr(b) = " << sqr(b) << ", " + << "sqr(c) = " << sqr(c) << nl; + + // n^th power + Info<< "pow(a, -1) = " << pow(a, -1) << ", " + << "pow(b, -1) = " << pow(b, -1) << ", " + << "pow(c, -1) = " << pow(c, -1) << nl; + + // Exponential + Info<< "exp(a) = " << exp(a) << ", " + << "exp(b) = " << exp(b) << ", " + << "exp(c) = " << exp(c) << nl; + + // Natural logarithm + Info<< "log(a) = " << log(a) << ", " + << "log(b) = " << log(b) << ", " + << "log(c) = " << log(c) << nl; + + } + #endif + // Make some changes { @@ -150,6 +279,7 @@ int main(int argc, char *argv[]) // Info<< "min/max = " << MinMax(fld1) << nl; Info<< "\nEnd\n" << endl; + return 0; } diff --git a/applications/test/dictionary2/Test-dictionary2.C b/applications/test/dictionary2/Test-dictionary2.C index d66e0d0590..b05c18c3c3 100644 --- a/applications/test/dictionary2/Test-dictionary2.C +++ b/applications/test/dictionary2/Test-dictionary2.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2017-2018 OpenCFD Ltd. + \\ / A nd | Copyright (C) 2017-2019 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -34,6 +34,8 @@ Description #include "IOobject.H" #include "IFstream.H" #include "dictionary.H" +#include "ops.H" +#include "scalarRange.H" #include "stringOps.H" using namespace Foam; @@ -108,6 +110,42 @@ scalar try_getScalar(const dictionary& dict, const word& k) } +// Try with getCheck +template +scalar try_getCheckScalar +( + const dictionary& dict, + const word& k, + const Predicate& pred +) +{ + scalar val(-GREAT); + + const bool throwingIOError = FatalIOError.throwExceptions(); + const bool throwingError = FatalError.throwExceptions(); + + try + { + val = dict.getCheck(k, pred); + Info<< "getCheck(" << k << ") = " << val << nl; + } + catch (const Foam::IOerror& err) + { + Info<< "getCheck(" << k << ") Caught FatalIOError " + << err << nl << endl; + } + catch (const Foam::error& err) + { + Info<< "getCheck(" << k << ") Caught FatalError " + << err << nl << endl; + } + FatalError.throwExceptions(throwingError); + FatalIOError.throwExceptions(throwingIOError); + + return val; +} + + // Try with *entry (from findEntry) and get scalar try_getScalar(const entry* eptr, const word& k) { @@ -311,6 +349,7 @@ int main(int argc, char *argv[]) IStringStream ( "good 3.14159;\n" + "negative -3.14159;\n" "empty;\n" // "bad text;\n" // always fails // "bad 3.14159 1234;\n" // fails for readScalar @@ -338,6 +377,26 @@ int main(int argc, char *argv[]) try_getScalar(dict2, "empty"); } + + // With getCheck + { + Info<< nl << "Test some input with getCheck()" << nl; + + try_getCheckScalar(dict2, "good", scalarRange::gt0()); + try_getCheckScalar(dict2, "negative", scalarRange::gt0()); + + try_getCheckScalar(dict2, "good", greaterOp1(0)); + try_getCheckScalar(dict2, "negative", greaterOp1(0)); + + Info<< nl << "with lambda" << nl; + try_getCheckScalar + ( + dict2, + "good", + [](const scalar x) { return x > 0; } + ); + } + // With findEntry and get { Info<< nl diff --git a/applications/test/regex1/Test-regex1.C b/applications/test/regex1/Test-regex1.C index ca5bfb60cd..9bc38d64b2 100644 --- a/applications/test/regex1/Test-regex1.C +++ b/applications/test/regex1/Test-regex1.C @@ -27,12 +27,16 @@ Description \*---------------------------------------------------------------------------*/ #include "argList.H" +#include "IOobject.H" #include "IOstreams.H" #include "IFstream.H" #include "Switch.H" +#include "SubStrings.H" #include "regExpCxx.H" +#ifndef _WIN32 #include "regExpPosix.H" +#endif using namespace Foam; @@ -83,6 +87,7 @@ static Ostream& operator<<(Ostream& os, const regExpCxx::results_type& sm) // Simple output of match groups +#ifndef _WIN32 static Ostream& operator<<(Ostream& os, const regExpPosix::results_type& sm) { for (std::smatch::size_type i = 1; i < sm.size(); ++i) @@ -92,6 +97,7 @@ static Ostream& operator<<(Ostream& os, const regExpPosix::results_type& sm) return os; } +#endif template @@ -209,7 +215,6 @@ void generalTests() } - template void testExpressions(const UList& tests) { @@ -293,11 +298,13 @@ int main(int argc, char *argv[]) "Test C++11 regular expressions" ); + #ifndef _WIN32 argList::addBoolOption ( "posix", "Test POSIX regular expressions" ); + #endif argList::addArgument("file"); argList::addArgument("..."); @@ -306,6 +313,17 @@ int main(int argc, char *argv[]) #include "setRootCase.H" + if (std::is_same::value) + { + Info<<"Foam::regExp uses C++11 regex" << nl << nl; + } + #ifndef _WIN32 + if (std::is_same::value) + { + Info<<"Foam::regExp uses POSIX regex" << nl << nl; + } + #endif + if (!args.count({"cxx", "posix"})) { Info<< "Specified one or more of -cxx, -posix" << nl; @@ -321,10 +339,12 @@ int main(int argc, char *argv[]) generalTests(); } + #ifndef _WIN32 if (args.found("posix")) { generalTests(); } + #endif } for (label argi = 1; argi < args.size(); ++argi) @@ -339,10 +359,12 @@ int main(int argc, char *argv[]) testExpressions(tests); } + #ifndef _WIN32 if (args.found("posix")) { testExpressions(tests); } + #endif } Info<< "\nDone" << nl << endl; diff --git a/applications/test/regex1/testRegexps2 b/applications/test/regex1/testRegexps2 new file mode 100644 index 0000000000..d74f8cec15 --- /dev/null +++ b/applications/test/regex1/testRegexps2 @@ -0,0 +1,23 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: v1812 | +| \\ / A nd | Web: www.OpenFOAM.com | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// Pattern, String +( + ( true "(U|k|epsilon)" "U" ) + ( false "(U|k|epsilon)" "alpha" ) + ( true "ab.*" "abc" ) + ( true ".*" "abc" ) + + ( true "div\(phi,alpha.*)" "div(phi,alpha.gas)" ) // quoting error + ( true "div\(phi,alpha.*\)" "div(phi,alpha.gas)" ) + ( true "div\(phi,alpha\..*\)" "div(phi,alpha.gas)" ) +) + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/applications/test/scalarPredicates/Test-scalarPredicates.C b/applications/test/scalarPredicates/Test-scalarPredicates.C index 22338e9174..587ca38b0d 100644 --- a/applications/test/scalarPredicates/Test-scalarPredicates.C +++ b/applications/test/scalarPredicates/Test-scalarPredicates.C @@ -36,6 +36,8 @@ Description #include "FlatOutput.H" #include "Tuple2.H" #include "StringStream.H" +#include "ops.H" +#include "bitSet.H" using namespace Foam; @@ -44,7 +46,7 @@ void doTest(const scalarList& values, const predicates::scalars& accept) { // Also tests that output is suppressed Info<<"Have: " << accept.size() << " predicates" << accept << endl; - Info<<"values: " << flatOutput(values) << endl; + Info<<"values: " << flatOutput(values) << endl; for (const scalar& value : values) { @@ -60,6 +62,30 @@ void doTest(const scalarList& values, const predicates::scalars& accept) } +template +void testPredicate(const scalarList& values, const Predicate& pred) +{ + bitSet matches; + + label i=0; + + for (const scalar& value : values) + { + if (pred(value)) + { + matches.set(i); + } + + ++i; + } + + IndirectList matched(values, matches.toc()); + + Info<< "matched: " << flatOutput(matched.addressing()) + << " = " << flatOutput(matched) << nl; +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Main program: @@ -149,6 +175,16 @@ int main(int argc, char *argv[]) } + Info<< nl << "Test with ops" << nl; + Info<<"values: " << flatOutput(values) << endl; + { + testPredicate(values, lessOp1(10)); + testPredicate(values, greaterOp1(100)); + + // Also with dissimilar type + testPredicate(values, lessEqOp1