zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Matrix: Switch n() and m() so that now m() = number of rows and n() = number of columns

Browse Source 
For consistency with current references and other numerics libraries
This commit is contained in:
Henry Weller
2016-03-10 22:53:09 +00:00
parent a018e83e75
commit 97489df33d
18 changed files with 227 additions and 227 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
src/OpenFOAM/dimensionSet/dimensionSets.C
View File

@ -234,17 +234,17 @@ Foam::dimensionSets::dimensionSets
valid_ = true; valid_ = true;
// Determine conversion from basic units to write units // Determine conversion from basic units to write units
for (label rowI = 0; rowI < conversion_.n(); rowI++) for (label rowI = 0; rowI < conversion_.m(); rowI++)
{ {
scalar* row = conversion_[rowI]; scalar* row = conversion_[rowI];
for (label columnI = 0; columnI < conversion_.m(); columnI++) for (label columnI = 0; columnI < conversion_.n(); columnI++)
{ {
const dimensionedScalar& dSet = units_[columnI]; const dimensionedScalar& dSet = units_[columnI];
row[columnI] = dSet.dimensions()[rowI]; row[columnI] = dSet.dimensions()[rowI];
} }
} }
conversionPivots_.setSize(conversion_.n()); conversionPivots_.setSize(conversion_.m());
LUDecompose(conversion_, conversionPivots_); LUDecompose(conversion_, conversionPivots_);
} }
} }

4
src/OpenFOAM/matrices/DiagonalMatrix/DiagonalMatrix.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -38,7 +38,7 @@ template<class Type>
template<class Form> template<class Form>
Foam::DiagonalMatrix<Type>::DiagonalMatrix(const Matrix<Form, Type>& a) Foam::DiagonalMatrix<Type>::DiagonalMatrix(const Matrix<Form, Type>& a)
: :
List<Type>(min(a.n(), a.m())) List<Type>(min(a.m(), a.n()))
{ {
forAll(*this, i) forAll(*this, i)
{ {

14
src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -43,7 +43,7 @@ Foam::LUscalarMatrix::LUscalarMatrix(const scalarSquareMatrix& matrix)
: :
scalarSquareMatrix(matrix), scalarSquareMatrix(matrix),
comm_(Pstream::worldComm), comm_(Pstream::worldComm),
pivotIndices_(n()) pivotIndices_(m())
{ {
LUDecompose(*this, pivotIndices_); LUDecompose(*this, pivotIndices_);
} }
@ -144,8 +144,8 @@ Foam::LUscalarMatrix::LUscalarMatrix
if (Pstream::master(comm_)) if (Pstream::master(comm_))
{ {
label nRows = n(); label nRows = m();
label nColumns = m(); label nColumns = n();
if (debug) if (debug)
{ {
@ -180,7 +180,7 @@ Foam::LUscalarMatrix::LUscalarMatrix
Pout<< endl; Pout<< endl;
} }
pivotIndices_.setSize(n()); pivotIndices_.setSize(m());
LUDecompose(*this, pivotIndices_); LUDecompose(*this, pivotIndices_);
} }
} }
@ -387,10 +387,10 @@ void Foam::LUscalarMatrix::convert
void Foam::LUscalarMatrix::printDiagonalDominance() const void Foam::LUscalarMatrix::printDiagonalDominance() const
{ {
for (label i=0; i<n(); i++) for (label i=0; i<m(); i++)
{ {
scalar sum = 0.0; scalar sum = 0.0;
for (label j=0; j<n(); j++) for (label j=0; j<m(); j++)
{ {
if (i != j) if (i != j)
{ {

4
src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -32,7 +32,7 @@ void Foam::LUscalarMatrix::solve(Field<Type>& sourceSol) const
{ {
if (Pstream::parRun()) if (Pstream::parRun())
{ {
Field<Type> completeSourceSol(n()); Field<Type> completeSourceSol(m());
if (Pstream::master(comm_)) if (Pstream::master(comm_))
{ {

114
src/OpenFOAM/matrices/Matrix/Matrix.C
View File

@ -96,9 +96,9 @@ Foam::Matrix<Form, Type>::Matrix(const label n, const label m, const Type& a)
{ {
Type* v = v_[0]; Type* v = v_[0];
label nm = nRows_*nCols_; label mn = nRows_*nCols_;
for (label i=0; i<nm; i++) for (label i=0; i<mn; i++)
{ {
v[i] = a; v[i] = a;
} }
@ -119,8 +119,8 @@ Foam::Matrix<Form, Type>::Matrix(const Matrix<Form, Type>& a)
Type* v = v_[0]; Type* v = v_[0];
const Type* av = a.v_[0]; const Type* av = a.v_[0];
label nm = nRows_*nCols_; label mn = nRows_*nCols_;
for (label i=0; i<nm; i++) for (label i=0; i<mn; i++)
{ {
v[i] = av[i]; v[i] = av[i];
} }
@ -162,11 +162,11 @@ template<class Form, class Type>
Form Foam::Matrix<Form, Type>::T() const Form Foam::Matrix<Form, Type>::T() const
{ {
const Matrix<Form, Type>& A = *this; const Matrix<Form, Type>& A = *this;
Form At(m(), n()); Form At(n(), m());
for (label i=0; i<n(); i++) for (label i=0; i<m(); i++)
{ {
for (label j=0; j<m(); j++) for (label j=0; j<n(); j++)
{ {
At[j][i] = A[i][j]; At[j][i] = A[i][j];
} }
@ -185,8 +185,8 @@ void Foam::Matrix<Form, Type>::operator=(const Type& t)
{ {
Type* v = v_[0]; Type* v = v_[0];
label nm = nRows_*nCols_; label mn = nRows_*nCols_;
for (label i=0; i<nm; i++) for (label i=0; i<mn; i++)
{ {
v[i] = t; v[i] = t;
} }
@ -200,7 +200,7 @@ void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& a)
if (this == &a) if (this == &a)
{ {
FatalErrorInFunction FatalErrorInFunction
<< "attempted assignment to self" << "attempted assigmnent to self"
<< abort(FatalError); << abort(FatalError);
} }
@ -217,8 +217,8 @@ void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& a)
Type* v = v_[0]; Type* v = v_[0];
const Type* av = a.v_[0]; const Type* av = a.v_[0];
label nm = nRows_*nCols_; label mn = nRows_*nCols_;
for (label i=0; i<nm; i++) for (label i=0; i<mn; i++)
{ {
v[i] = av[i]; v[i] = av[i];
} }
@ -231,14 +231,14 @@ void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& a)
template<class Form, class Type> template<class Form, class Type>
const Type& Foam::max(const Matrix<Form, Type>& a) const Type& Foam::max(const Matrix<Form, Type>& a)
{ {
label nm = a.n()*a.m(); label mn = a.m()*a.n();
if (nm) if (mn)
{ {
label curMaxI = 0; label curMaxI = 0;
const Type* v = a[0]; const Type* v = a[0];
for (label i=1; i<nm; i++) for (label i=1; i<mn; i++)
{ {
if (v[i] > v[curMaxI]) if (v[i] > v[curMaxI])
{ {
@ -263,14 +263,14 @@ const Type& Foam::max(const Matrix<Form, Type>& a)
template<class Form, class Type> template<class Form, class Type>
const Type& Foam::min(const Matrix<Form, Type>& a) const Type& Foam::min(const Matrix<Form, Type>& a)
{ {
label nm = a.n()*a.m(); label mn = a.m()*a.n();
if (nm) if (mn)
{ {
label curMinI = 0; label curMinI = 0;
const Type* v = a[0]; const Type* v = a[0];
for (label i=1; i<nm; i++) for (label i=1; i<mn; i++)
{ {
if (v[i] < v[curMinI]) if (v[i] < v[curMinI])
{ {
@ -297,15 +297,15 @@ const Type& Foam::min(const Matrix<Form, Type>& a)
template<class Form, class Type> template<class Form, class Type>
Form Foam::operator-(const Matrix<Form, Type>& a) Form Foam::operator-(const Matrix<Form, Type>& a)
{ {
Form na(a.n(), a.m()); Form na(a.m(), a.n());
if (a.n() && a.m()) if (a.m() && a.n())
{ {
Type* nav = na[0]; Type* nav = na[0];
const Type* av = a[0]; const Type* av = a[0];
label nm = a.n()*a.m(); label mn = a.m()*a.n();
for (label i=0; i<nm; i++) for (label i=0; i<mn; i++)
{ {
nav[i] = -av[i]; nav[i] = -av[i];
} }
@ -318,30 +318,30 @@ Form Foam::operator-(const Matrix<Form, Type>& a)
template<class Form, class Type> template<class Form, class Type>
Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b) Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
{ {
if (a.n() != b.n())
{
FatalErrorInFunction
<< "attempted add matrices with different number of rows: "
<< a.n() << ", " << b.n()
<< abort(FatalError);
}
if (a.m() != b.m()) if (a.m() != b.m())
{ {
FatalErrorInFunction FatalErrorInFunction
<< "attempted add matrices with different number of columns: " << "attempted add matrices with different number of rows: "
<< a.m() << ", " << b.m() << a.m() << ", " << b.m()
<< abort(FatalError); << abort(FatalError);
} }
Form ab(a.n(), a.m()); if (a.n() != b.n())
{
FatalErrorInFunction
<< "attempted add matrices with different number of columns: "
<< a.n() << ", " << b.n()
<< abort(FatalError);
}
Form ab(a.m(), a.n());
Type* abv = ab[0]; Type* abv = ab[0];
const Type* av = a[0]; const Type* av = a[0];
const Type* bv = b[0]; const Type* bv = b[0];
label nm = a.n()*a.m(); label mn = a.m()*a.n();
for (label i=0; i<nm; i++) for (label i=0; i<mn; i++)
{ {
abv[i] = av[i] + bv[i]; abv[i] = av[i] + bv[i];
} }
@ -353,30 +353,30 @@ Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
template<class Form, class Type> template<class Form, class Type>
Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b) Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
{ {
if (a.n() != b.n())
{
FatalErrorInFunction
<< "attempted add matrices with different number of rows: "
<< a.n() << ", " << b.n()
<< abort(FatalError);
}
if (a.m() != b.m()) if (a.m() != b.m())
{ {
FatalErrorInFunction FatalErrorInFunction
<< "attempted add matrices with different number of columns: " << "attempted add matrices with different number of rows: "
<< a.m() << ", " << b.m() << a.m() << ", " << b.m()
<< abort(FatalError); << abort(FatalError);
} }
Form ab(a.n(), a.m()); if (a.n() != b.n())
{
FatalErrorInFunction
<< "attempted add matrices with different number of columns: "
<< a.n() << ", " << b.n()
<< abort(FatalError);
}
Form ab(a.m(), a.n());
Type* abv = ab[0]; Type* abv = ab[0];
const Type* av = a[0]; const Type* av = a[0];
const Type* bv = b[0]; const Type* bv = b[0];
label nm = a.n()*a.m(); label mn = a.m()*a.n();
for (label i=0; i<nm; i++) for (label i=0; i<mn; i++)
{ {
abv[i] = av[i] - bv[i]; abv[i] = av[i] - bv[i];
} }
@ -388,15 +388,15 @@ Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
template<class Form, class Type> template<class Form, class Type>
Form Foam::operator*(const scalar s, const Matrix<Form, Type>& a) Form Foam::operator*(const scalar s, const Matrix<Form, Type>& a)
{ {
Form sa(a.n(), a.m()); Form sa(a.m(), a.n());
if (a.n() && a.m()) if (a.m() && a.n())
{ {
Type* sav = sa[0]; Type* sav = sa[0];
const Type* av = a[0]; const Type* av = a[0];
label nm = a.n()*a.m(); label mn = a.m()*a.n();
for (label i=0; i<nm; i++) for (label i=0; i<mn; i++)
{ {
sav[i] = s*av[i]; sav[i] = s*av[i];
} }
@ -409,24 +409,24 @@ Form Foam::operator*(const scalar s, const Matrix<Form, Type>& a)
template<class Form, class Type> template<class Form, class Type>
Form Foam::operator*(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b) Form Foam::operator*(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
{ {
if (a.m() != b.n()) if (a.n() != b.m())
{ {
FatalErrorInFunction FatalErrorInFunction
<< "attempted to multiply incompatible matrices:" << nl << "attempted to multiply incompatible matrices:" << nl
<< "Matrix A : " << a.n() << " rows, " << a.m() << " columns" << nl << "Matrix A : " << a.m() << " rows, " << a.n() << " columns" << nl
<< "Matrix B : " << b.n() << " rows, " << b.m() << " columns" << nl << "Matrix B : " << b.m() << " rows, " << b.n() << " columns" << nl
<< "In order to multiply matrices, columns of A must equal " << "In order to multiply matrices, columns of A must equal "
<< "rows of B" << "rows of B"
<< abort(FatalError); << abort(FatalError);
} }
Form ab(a.n(), b.m(), scalar(0)); Form ab(a.m(), b.n(), scalar(0));
for (label i = 0; i < ab.n(); i++) for (label i = 0; i < ab.m(); i++)
{ {
for (label j = 0; j < ab.m(); j++) for (label j = 0; j < ab.n(); j++)
{ {
for (label l = 0; l < b.n(); l++) for (label l = 0; l < b.m(); l++)
{ {
ab[i][j] += a[i][l]*b[l][j]; ab[i][j] += a[i][l]*b[l][j];
} }

4
src/OpenFOAM/matrices/Matrix/Matrix.H
View File

@ -125,10 +125,10 @@ public:
// Access // Access
//- Return the number of rows //- Return the number of rows
inline label n() const; inline label m() const;
//- Return the number of columns //- Return the number of columns
inline label m() const; inline label n() const;
//- Return the number of elements in matrix (n*m) //- Return the number of elements in matrix (n*m)
inline label size() const; inline label size() const;

4
src/OpenFOAM/matrices/Matrix/MatrixI.H
View File

@ -53,14 +53,14 @@ inline const Foam::Matrix<Form, Type>& Foam::Matrix<Form, Type>::null()
//- Return the number of rows //- Return the number of rows
template<class Form, class Type> template<class Form, class Type>
inline Foam::label Foam::Matrix<Form, Type>::n() const inline Foam::label Foam::Matrix<Form, Type>::m() const
{ {
return nRows_; return nRows_;
} }
template<class Form, class Type> template<class Form, class Type>
inline Foam::label Foam::Matrix<Form, Type>::m() const inline Foam::label Foam::Matrix<Form, Type>::n() const
{ {
return nCols_; return nCols_;
} }

38
src/OpenFOAM/matrices/Matrix/MatrixIO.C
View File

@ -62,7 +62,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
M.nRows_ = firstToken.labelToken(); M.nRows_ = firstToken.labelToken();
M.nCols_ = readLabel(is); M.nCols_ = readLabel(is);
label nm = M.nRows_*M.nCols_; label mn = M.nRows_*M.nCols_;
// Read list contents depending on data format // Read list contents depending on data format
if (is.format() == IOstream::ASCII || !contiguous<Type>()) if (is.format() == IOstream::ASCII || !contiguous<Type>())
@ -70,7 +70,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
// Read beginning of contents // Read beginning of contents
char listDelimiter = is.readBeginList("Matrix"); char listDelimiter = is.readBeginList("Matrix");
if (nm) if (mn)
{ {
M.allocate(); M.allocate();
Type* v = M.v_[0]; Type* v = M.v_[0];
@ -80,11 +80,11 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
label k = 0; label k = 0;
// loop over rows // loop over rows
for (label i=0; i<M.n(); i++) for (label i=0; i<M.m(); i++)
{ {
listDelimiter = is.readBeginList("MatrixRow"); listDelimiter = is.readBeginList("MatrixRow");
for (label j=0; j<M.m(); j++) for (label j=0; j<M.n(); j++)
{ {
is >> v[k++]; is >> v[k++];
@ -109,7 +109,7 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
"reading the single entry" "reading the single entry"
); );
for (label i=0; i<nm; i++) for (label i=0; i<mn; i++)
{ {
v[i] = element; v[i] = element;
} }
@ -121,12 +121,12 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
} }
else else
{ {
if (nm) if (mn)
{ {
M.allocate(); M.allocate();
Type* v = M.v_[0]; Type* v = M.v_[0];
is.read(reinterpret_cast<char*>(v), nm*sizeof(Type)); is.read(reinterpret_cast<char*>(v), mn*sizeof(Type));
is.fatalCheck is.fatalCheck
( (
@ -151,24 +151,24 @@ Foam::Istream& Foam::operator>>(Istream& is, Matrix<Form, Type>& M)
template<class Form, class Type> template<class Form, class Type>
Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M) Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
{ {
label nm = M.nRows_*M.nCols_; label mn = M.nRows_*M.nCols_;
os << M.n() << token::SPACE << M.m(); os << M.m() << token::SPACE << M.n();
// Write list contents depending on data format // Write list contents depending on data format
if (os.format() == IOstream::ASCII || !contiguous<Type>()) if (os.format() == IOstream::ASCII || !contiguous<Type>())
{ {
if (nm) if (mn)
{ {
bool uniform = false; bool uniform = false;
const Type* v = M.v_[0]; const Type* v = M.v_[0];
if (nm > 1 && contiguous<Type>()) if (mn > 1 && contiguous<Type>())
{ {
uniform = true; uniform = true;
for (label i=0; i< nm; i++) for (label i=0; i< mn; i++)
{ {
if (v[i] != v[0]) if (v[i] != v[0])
{ {
@ -189,7 +189,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
// Write end of contents delimiter // Write end of contents delimiter
os << token::END_BLOCK; os << token::END_BLOCK;
} }
else if (nm < 10 && contiguous<Type>()) else if (mn < 10 && contiguous<Type>())
{ {
// Write size of list and start contents delimiter // Write size of list and start contents delimiter
os << token::BEGIN_LIST; os << token::BEGIN_LIST;
@ -197,12 +197,12 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
label k = 0; label k = 0;
// loop over rows // loop over rows
for (label i=0; i< M.n(); i++) for (label i=0; i< M.m(); i++)
{ {
os << token::BEGIN_LIST; os << token::BEGIN_LIST;
// Write row // Write row
for (label j=0; j< M.m(); j++) for (label j=0; j< M.n(); j++)
{ {
if (j > 0) os << token::SPACE; if (j > 0) os << token::SPACE;
os << v[k++]; os << v[k++];
@ -222,12 +222,12 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
label k = 0; label k = 0;
// loop over rows // loop over rows
for (label i=0; i< M.n(); i++) for (label i=0; i< M.m(); i++)
{ {
os << nl << token::BEGIN_LIST; os << nl << token::BEGIN_LIST;
// Write row // Write row
for (label j=0; j< M.m(); j++) for (label j=0; j< M.n(); j++)
{ {
os << nl << v[k++]; os << nl << v[k++];
} }
@ -246,9 +246,9 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const Matrix<Form, Type>& M)
} }
else else
{ {
if (nm) if (mn)
{ {
os.write(reinterpret_cast<const char*>(M.v_[0]), nm*sizeof(Type)); os.write(reinterpret_cast<const char*>(M.v_[0]), mn*sizeof(Type));
} }
} }

8
src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2013 OpenFOAM Foundation \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -37,7 +37,7 @@ Foam::scalar Foam::detDecomposed
{ {
scalar diagProduct = 1.0; scalar diagProduct = 1.0;
for (label i = 0; i < matrix.n(); ++i) for (label i = 0; i < matrix.m(); ++i)
{ {
diagProduct *= matrix[i][i]; diagProduct *= matrix[i][i];
} }
@ -51,7 +51,7 @@ Foam::scalar Foam::det(const SquareMatrix<Type>& matrix)
{ {
SquareMatrix<Type> matrixTmp = matrix; SquareMatrix<Type> matrixTmp = matrix;
labelList pivotIndices(matrix.n()); labelList pivotIndices(matrix.m());
label sign; label sign;
LUDecompose(matrixTmp, pivotIndices, sign); LUDecompose(matrixTmp, pivotIndices, sign);
@ -62,7 +62,7 @@ Foam::scalar Foam::det(const SquareMatrix<Type>& matrix)
template<class Type> template<class Type>
Foam::scalar Foam::det(SquareMatrix<Type>& matrix) Foam::scalar Foam::det(SquareMatrix<Type>& matrix)
{ {
labelList pivotIndices(matrix.n()); labelList pivotIndices(matrix.m());
label sign; label sign;
LUDecompose(matrix, pivotIndices, sign); LUDecompose(matrix, pivotIndices, sign);

8
src/OpenFOAM/matrices/SymmetricSquareMatrix/SymmetricSquareMatrix.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -33,9 +33,9 @@ Foam::SymmetricSquareMatrix<Type> Foam::invDecomposed
const SymmetricSquareMatrix<Type>& matrix const SymmetricSquareMatrix<Type>& matrix
) )
{ {
SymmetricSquareMatrix<Type> inv(matrix.n(), matrix.n(), 0.0); SymmetricSquareMatrix<Type> inv(matrix.m(), matrix.m(), 0.0);
for (label i = 0; i < matrix.n(); ++i) for (label i = 0; i < matrix.m(); ++i)
{ {
inv[i][i] = 1.0/matrix[i][i]; inv[i][i] = 1.0/matrix[i][i];
@ -75,7 +75,7 @@ Foam::scalar Foam::detDecomposed(const SymmetricSquareMatrix<Type>& matrix)
{ {
scalar diagProduct = 1.0; scalar diagProduct = 1.0;
for (label i = 0; i < matrix.n(); ++i) for (label i = 0; i < matrix.m(); ++i)
{ {
diagProduct *= matrix[i][i]; diagProduct *= matrix[i][i];
} }

113
src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -33,39 +33,39 @@ License
Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition) Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
: :
U_(A), U_(A),
V_(A.m(), A.m()), V_(A.n(), A.n()),
S_(A.m()), S_(A.n()),
VSinvUt_(A.m(), A.n()), VSinvUt_(A.n(), A.m()),
nZeros_(0) nZeros_(0)
{ {
// SVDcomp to find U_, V_ and S_ - the singular values // SVDcomp to find U_, V_ and S_ - the singular values
const label Um = U_.m();
const label Un = U_.n(); const label Un = U_.n();
const label Um = U_.m();
scalarList rv1(Um); scalarList rv1(Un);
scalar g = 0; scalar g = 0;
scalar scale = 0; scalar scale = 0;
scalar s = 0; scalar s = 0;
scalar anorm = 0; scalar anorm = 0;
label l = 0; label l = 0;
for (label i = 0; i < Um; i++) for (label i = 0; i < Un; i++)
{ {
l = i+2; l = i+2;
rv1[i] = scale*g; rv1[i] = scale*g;
g = s = scale = 0; g = s = scale = 0;
if (i < Un) if (i < Um)
{ {
for (label k = i; k < Un; k++) for (label k = i; k < Um; k++)
{ {
scale += mag(U_[k][i]); scale += mag(U_[k][i]);
} }
if (scale != 0) if (scale != 0)
{ {
for (label k = i; k < Un; k++) for (label k = i; k < Um; k++)
{ {
U_[k][i] /= scale; U_[k][i] /= scale;
s += U_[k][i]*U_[k][i]; s += U_[k][i]*U_[k][i];
@ -76,22 +76,22 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
scalar h = f*g - s; scalar h = f*g - s;
U_[i][i] = f - g; U_[i][i] = f - g;
for (label j = l-1; j < Um; j++) for (label j = l-1; j < Un; j++)
{ {
s = 0; s = 0;
for (label k = i; k < Un; k++) for (label k = i; k < Um; k++)
{ {
s += U_[k][i]*U_[k][j]; s += U_[k][i]*U_[k][j];
} }
f = s/h; f = s/h;
for (label k = i; k < A.n(); k++) for (label k = i; k < A.m(); k++)
{ {
U_[k][j] += f*U_[k][i]; U_[k][j] += f*U_[k][i];
} }
} }
for (label k = i; k < Un; k++) for (label k = i; k < Um; k++)
{ {
U_[k][i] *= scale; U_[k][i] *= scale;
} }
@ -102,16 +102,16 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
g = s = scale = 0; g = s = scale = 0;
if (i+1 <= Un && i != Um) if (i+1 <= Um && i != Un)
{ {
for (label k = l-1; k < Um; k++) for (label k = l-1; k < Un; k++)
{ {
scale += mag(U_[i][k]); scale += mag(U_[i][k]);
} }
if (scale != 0) if (scale != 0)
{ {
for (label k=l-1; k < Um; k++) for (label k=l-1; k < Un; k++)
{ {
U_[i][k] /= scale; U_[i][k] /= scale;
s += U_[i][k]*U_[i][k]; s += U_[i][k]*U_[i][k];
@ -122,25 +122,25 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
scalar h = f*g - s; scalar h = f*g - s;
U_[i][l-1] = f - g; U_[i][l-1] = f - g;
for (label k = l-1; k < Um; k++) for (label k = l-1; k < Un; k++)
{ {
rv1[k] = U_[i][k]/h; rv1[k] = U_[i][k]/h;
} }
for (label j = l-1; j < Un; j++) for (label j = l-1; j < Um; j++)
{ {
s = 0; s = 0;
for (label k = l-1; k < Um; k++) for (label k = l-1; k < Un; k++)
{ {
s += U_[j][k]*U_[i][k]; s += U_[j][k]*U_[i][k];
} }
for (label k = l-1; k < Um; k++) for (label k = l-1; k < Un; k++)
{ {
U_[j][k] += s*rv1[k]; U_[j][k] += s*rv1[k];
} }
} }
for (label k = l-1; k < Um; k++) for (label k = l-1; k < Un; k++)
{ {
U_[i][k] *= scale; U_[i][k] *= scale;
} }
@ -150,33 +150,33 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
anorm = max(anorm, mag(S_[i]) + mag(rv1[i])); anorm = max(anorm, mag(S_[i]) + mag(rv1[i]));
} }
for (label i = Um-1; i >= 0; i--) for (label i = Un-1; i >= 0; i--)
{ {
if (i < Um-1) if (i < Un-1)
{ {
if (g != 0) if (g != 0)
{ {
for (label j = l; j < Um; j++) for (label j = l; j < Un; j++)
{ {
V_[j][i] = (U_[i][j]/U_[i][l])/g; V_[j][i] = (U_[i][j]/U_[i][l])/g;
} }
for (label j=l; j < Um; j++) for (label j=l; j < Un; j++)
{ {
s = 0; s = 0;
for (label k = l; k < Um; k++) for (label k = l; k < Un; k++)
{ {
s += U_[i][k]*V_[k][j]; s += U_[i][k]*V_[k][j];
} }
for (label k = l; k < Um; k++) for (label k = l; k < Un; k++)
{ {
V_[k][j] += s*V_[k][i]; V_[k][j] += s*V_[k][i];
} }
} }
} }
for (label j = l; j < Um;j++) for (label j = l; j < Un;j++)
{ {
V_[i][j] = V_[j][i] = 0.0; V_[i][j] = V_[j][i] = 0.0;
} }
@ -187,12 +187,12 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
l = i; l = i;
} }
for (label i = min(Um, Un) - 1; i >= 0; i--) for (label i = min(Un, Um) - 1; i >= 0; i--)
{ {
l = i+1; l = i+1;
g = S_[i]; g = S_[i];
for (label j = l; j < Um; j++) for (label j = l; j < Un; j++)
{ {
U_[i][j] = 0.0; U_[i][j] = 0.0;
} }
@ -201,30 +201,30 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
{ {
g = 1.0/g; g = 1.0/g;
for (label j = l; j < Um; j++) for (label j = l; j < Un; j++)
{ {
s = 0; s = 0;
for (label k = l; k < Un; k++) for (label k = l; k < Um; k++)
{ {
s += U_[k][i]*U_[k][j]; s += U_[k][i]*U_[k][j];
} }
scalar f = (s/U_[i][i])*g; scalar f = (s/U_[i][i])*g;
for (label k = i; k < Un; k++) for (label k = i; k < Um; k++)
{ {
U_[k][j] += f*U_[k][i]; U_[k][j] += f*U_[k][i];
} }
} }
for (label j = i; j < Un; j++) for (label j = i; j < Um; j++)
{ {
U_[j][i] *= g; U_[j][i] *= g;
} }
} }
else else
{ {
for (label j = i; j < Un; j++) for (label j = i; j < Um; j++)
{ {
U_[j][i] = 0.0; U_[j][i] = 0.0;
} }
@ -233,22 +233,22 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
++U_[i][i]; ++U_[i][i];
} }
for (label k = Um-1; k >= 0; k--) for (label k = Un-1; k >= 0; k--)
{ {
for (label its = 0; its < 35; its++) for (label its = 0; its < 35; its++)
{ {
bool flag = true; bool flag = true;
label nm; label mn;
for (l = k; l >= 0; l--) for (l = k; l >= 0; l--)
{ {
nm = l-1; mn = l-1;
if (mag(rv1[l]) + anorm == anorm) if (mag(rv1[l]) + anorm == anorm)
{ {
flag = false; flag = false;
break; break;
} }
if (mag(S_[nm]) + anorm == anorm) break; if (mag(S_[mn]) + anorm == anorm) break;
} }
if (flag) if (flag)
@ -269,11 +269,11 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
c = g*h; c = g*h;
s = -f*h; s = -f*h;
for (label j = 0; j < Un; j++) for (label j = 0; j < Um; j++)
{ {
scalar y = U_[j][nm]; scalar y = U_[j][mn];
scalar z = U_[j][i]; scalar z = U_[j][i];
U_[j][nm] = y*c + z*s; U_[j][mn] = y*c + z*s;
U_[j][i] = z*c - y*s; U_[j][i] = z*c - y*s;
} }
} }
@ -287,21 +287,21 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
{ {
S_[k] = -z; S_[k] = -z;
for (label j = 0; j < Um; j++) V_[j][k] = -V_[j][k]; for (label j = 0; j < Un; j++) V_[j][k] = -V_[j][k];
} }
break; break;
} }
if (its == 34) if (its == 34)
{ {
WarningInFunction WarningInFunction
<< "no convergence in 35 SVD iterations" << "No convergence in 35 SVD iterations"
<< endl; << endl;
} }
scalar x = S_[l]; scalar x = S_[l];
nm = k-1; mn = k-1;
scalar y = S_[nm]; scalar y = S_[mn];
g = rv1[nm]; g = rv1[mn];
scalar h = rv1[k]; scalar h = rv1[k];
scalar f = ((y - z)*(y + z) + (g - h)*(g + h))/(2.0*h*y); scalar f = ((y - z)*(y + z) + (g - h)*(g + h))/(2.0*h*y);
g = sqrtSumSqr(f, scalar(1)); g = sqrtSumSqr(f, scalar(1));
@ -309,7 +309,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
scalar c = 1.0; scalar c = 1.0;
s = 1.0; s = 1.0;
for (label j = l; j <= nm; j++) for (label j = l; j <= mn; j++)
{ {
label i = j + 1; label i = j + 1;
g = rv1[i]; g = rv1[i];
@ -325,7 +325,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
h = y*s; h = y*s;
y *= c; y *= c;
for (label jj = 0; jj < Um; jj++) for (label jj = 0; jj < Un; jj++)
{ {
x = V_[jj][j]; x = V_[jj][j];
z = V_[jj][i]; z = V_[jj][i];
@ -344,7 +344,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
f = c*g + s*y; f = c*g + s*y;
x = c*y - s*g; x = c*y - s*g;
for (label jj=0; jj < Un; jj++) for (label jj=0; jj < Um; jj++)
{ {
y = U_[jj][j]; y = U_[jj][j];
z = U_[jj][i]; z = U_[jj][i];
@ -358,7 +358,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
} }
} }
// zero singular values that are less than minCondition*maxS // Zero singular values that are less than minCondition*maxS
const scalar minS = minCondition*S_[findMax(S_)]; const scalar minS = minCondition*S_[findMax(S_)];
forAll(S_, i) forAll(S_, i)
{ {
@ -370,17 +370,18 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
} }
} }
// now multiply out to find the pseudo inverse of A, VSinvUt_ // Now multiply out to find the pseudo inverse of A, VSinvUt_
multiply(VSinvUt_, V_, inv(S_), U_.T()); multiply(VSinvUt_, V_, inv(S_), U_.T());
// test SVD // test SVD
/*scalarRectangularMatrix SVDA(A.n(), A.m()); /*
scalarRectangularMatrix SVDA(A.m(), A.n());
multiply(SVDA, U_, S_, transpose(V_)); multiply(SVDA, U_, S_, transpose(V_));
scalar maxDiff = 0; scalar maxDiff = 0;
scalar diff = 0; scalar diff = 0;
for (label i = 0; i < A.n(); i++) for (label i = 0; i < A.m(); i++)
{ {
for (label j = 0; j < A.m(); j++) for (label j = 0; j < A.n(); j++)
{ {
diff = mag(A[i][j] - SVDA[i][j]); diff = mag(A[i][j] - SVDA[i][j]);
if (diff > maxDiff) maxDiff = diff; if (diff > maxDiff) maxDiff = diff;

3
src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -45,7 +45,6 @@ namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class SVD Declaration Class SVD Declaration
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/

64
src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -46,17 +46,17 @@ void Foam::LUDecompose
label& sign label& sign
) )
{ {
label n = matrix.n(); label m = matrix.m();
scalar vv[n]; scalar vv[m];
sign = 1; sign = 1;
for (label i=0; i<n; i++) for (label i=0; i<m; i++)
{ {
scalar largestCoeff = 0.0; scalar largestCoeff = 0.0;
scalar temp; scalar temp;
const scalar* __restrict__ matrixi = matrix[i]; const scalar* __restrict__ matrixi = matrix[i];
for (label j=0; j<n; j++) for (label j=0; j<m; j++)
{ {
if ((temp = mag(matrixi[j])) > largestCoeff) if ((temp = mag(matrixi[j])) > largestCoeff)
{ {
@ -73,7 +73,7 @@ void Foam::LUDecompose
vv[i] = 1.0/largestCoeff; vv[i] = 1.0/largestCoeff;
} }
for (label j=0; j<n; j++) for (label j=0; j<m; j++)
{ {
scalar* __restrict__ matrixj = matrix[j]; scalar* __restrict__ matrixj = matrix[j];
@ -92,7 +92,7 @@ void Foam::LUDecompose
label iMax = 0; label iMax = 0;
scalar largestCoeff = 0.0; scalar largestCoeff = 0.0;
for (label i=j; i<n; i++) for (label i=j; i<m; i++)
{ {
scalar* __restrict__ matrixi = matrix[i]; scalar* __restrict__ matrixi = matrix[i];
scalar sum = matrixi[j]; scalar sum = matrixi[j];
@ -118,7 +118,7 @@ void Foam::LUDecompose
{ {
scalar* __restrict__ matrixiMax = matrix[iMax]; scalar* __restrict__ matrixiMax = matrix[iMax];
for (label k=0; k<n; k++) for (label k=0; k<m; k++)
{ {
Swap(matrixj[k], matrixiMax[k]); Swap(matrixj[k], matrixiMax[k]);
} }
@ -132,11 +132,11 @@ void Foam::LUDecompose
matrixj[j] = SMALL; matrixj[j] = SMALL;
} }
if (j != n-1) if (j != m-1)
{ {
scalar rDiag = 1.0/matrixj[j]; scalar rDiag = 1.0/matrixj[j];
for (label i=j+1; i<n; i++) for (label i=j+1; i<m; i++)
{ {
matrix[i][j] *= rDiag; matrix[i][j] *= rDiag;
} }
@ -148,7 +148,7 @@ void Foam::LUDecompose
void Foam::LUDecompose(scalarSymmetricSquareMatrix& matrix) void Foam::LUDecompose(scalarSymmetricSquareMatrix& matrix)
{ {
// Store result in upper triangular part of matrix // Store result in upper triangular part of matrix
label size = matrix.n(); label size = matrix.m();
// Set upper triangular parts to zero. // Set upper triangular parts to zero.
for (label j = 0; j < size; j++) for (label j = 0; j < size; j++)
@ -205,32 +205,32 @@ void Foam::multiply
const scalarRectangularMatrix& C const scalarRectangularMatrix& C
) )
{ {
if (A.m() != B.n()) if (A.n() != B.m())
{ {
FatalErrorInFunction FatalErrorInFunction
<< "A and B must have identical inner dimensions but A.m = " << "A and B must have identical inner dimensions but A.n = "
<< A.m() << " and B.n = " << B.n() << A.n() << " and B.m = " << B.m()
<< abort(FatalError); << abort(FatalError);
} }
if (B.m() != C.n()) if (B.n() != C.m())
{ {
FatalErrorInFunction FatalErrorInFunction
<< "B and C must have identical inner dimensions but B.m = " << "B and C must have identical inner dimensions but B.n = "
<< B.m() << " and C.n = " << C.n() << B.n() << " and C.m = " << C.m()
<< abort(FatalError); << abort(FatalError);
} }
ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0)); ans = scalarRectangularMatrix(A.m(), C.n(), scalar(0));
for (label i = 0; i < A.n(); i++) for (label i = 0; i < A.m(); i++)
{ {
for (label g = 0; g < C.m(); g++) for (label g = 0; g < C.n(); g++)
{ {
for (label l = 0; l < C.n(); l++) for (label l = 0; l < C.m(); l++)
{ {
scalar ab = 0; scalar ab = 0;
for (label j = 0; j < A.m(); j++) for (label j = 0; j < A.n(); j++)
{ {
ab += A[i][j]*B[j][l]; ab += A[i][j]*B[j][l];
} }
@ -249,29 +249,29 @@ void Foam::multiply
const scalarRectangularMatrix& C const scalarRectangularMatrix& C
) )
{ {
if (A.m() != B.size()) if (A.n() != B.size())
{ {
FatalErrorInFunction FatalErrorInFunction
<< "A and B must have identical inner dimensions but A.m = " << "A and B must have identical inner dimensions but A.n = "
<< A.m() << " and B.n = " << B.size() << A.n() << " and B.m = " << B.size()
<< abort(FatalError); << abort(FatalError);
} }
if (B.size() != C.n()) if (B.size() != C.m())
{ {
FatalErrorInFunction FatalErrorInFunction
<< "B and C must have identical inner dimensions but B.m = " << "B and C must have identical inner dimensions but B.n = "
<< B.size() << " and C.n = " << C.n() << B.size() << " and C.m = " << C.m()
<< abort(FatalError); << abort(FatalError);
} }
ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0)); ans = scalarRectangularMatrix(A.m(), C.n(), scalar(0));
for (label i = 0; i < A.n(); i++) for (label i = 0; i < A.m(); i++)
{ {
for (label g = 0; g < C.m(); g++) for (label g = 0; g < C.n(); g++)
{ {
for (label l = 0; l < C.n(); l++) for (label l = 0; l < C.m(); l++)
{ {
ans[i][g] += C[l][g] * A[i][l]*B[l]; ans[i][g] += C[l][g] * A[i][l]*B[l];
} }

50
src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -36,16 +36,16 @@ void Foam::solve
List<Type>& sourceSol List<Type>& sourceSol
) )
{ {
label n = tmpMatrix.n(); label m = tmpMatrix.m();
// Elimination // Elimination
for (label i=0; i<n; i++) for (label i=0; i<m; i++)
{ {
label iMax = i; label iMax = i;
scalar largestCoeff = mag(tmpMatrix[iMax][i]); scalar largestCoeff = mag(tmpMatrix[iMax][i]);
// Swap entries around to find a good pivot // Swap entries around to find a good pivot
for (label j=i+1; j<n; j++) for (label j=i+1; j<m; j++)
{ {
if (mag(tmpMatrix[j][i]) > largestCoeff) if (mag(tmpMatrix[j][i]) > largestCoeff)
{ {
@ -58,7 +58,7 @@ void Foam::solve
{ {
//Info<< "Pivoted on " << i << " " << iMax << endl; //Info<< "Pivoted on " << i << " " << iMax << endl;
for (label k=i; k<n; k++) for (label k=i; k<m; k++)
{ {
Swap(tmpMatrix[i][k], tmpMatrix[iMax][k]); Swap(tmpMatrix[i][k], tmpMatrix[iMax][k]);
} }
@ -74,11 +74,11 @@ void Foam::solve
} }
// Reduce to upper triangular form // Reduce to upper triangular form
for (label j=i+1; j<n; j++) for (label j=i+1; j<m; j++)
{ {
sourceSol[j] -= sourceSol[i]*(tmpMatrix[j][i]/tmpMatrix[i][i]); sourceSol[j] -= sourceSol[i]*(tmpMatrix[j][i]/tmpMatrix[i][i]);
for (label k=n-1; k>=i; k--) for (label k=m-1; k>=i; k--)
{ {
tmpMatrix[j][k] -= tmpMatrix[j][k] -=
tmpMatrix[i][k]*tmpMatrix[j][i]/tmpMatrix[i][i]; tmpMatrix[i][k]*tmpMatrix[j][i]/tmpMatrix[i][i];
@ -87,11 +87,11 @@ void Foam::solve
} }
// Back-substitution // Back-substitution
for (label j=n-1; j>=0; j--) for (label j=m-1; j>=0; j--)
{ {
Type ntempvec = pTraits<Type>::zero; Type ntempvec = pTraits<Type>::zero;
for (label k=j+1; k<n; k++) for (label k=j+1; k<m; k++)
{ {
ntempvec += tmpMatrix[j][k]*sourceSol[k]; ntempvec += tmpMatrix[j][k]*sourceSol[k];
} }
@ -123,11 +123,11 @@ void Foam::LUBacksubstitute
List<Type>& sourceSol List<Type>& sourceSol
) )
{ {
label n = luMatrix.n(); label m = luMatrix.m();
label ii = 0; label ii = 0;
for (label i=0; i<n; i++) for (label i=0; i<m; i++)
{ {
label ip = pivotIndices[i]; label ip = pivotIndices[i];
Type sum = sourceSol[ip]; Type sum = sourceSol[ip];
@ -149,12 +149,12 @@ void Foam::LUBacksubstitute
sourceSol[i] = sum; sourceSol[i] = sum;
} }
for (label i=n-1; i>=0; i--) for (label i=m-1; i>=0; i--)
{ {
Type sum = sourceSol[i]; Type sum = sourceSol[i];
const scalar* __restrict__ luMatrixi = luMatrix[i]; const scalar* __restrict__ luMatrixi = luMatrix[i];
for (label j=i+1; j<n; j++) for (label j=i+1; j<m; j++)
{ {
sum -= luMatrixi[j]*sourceSol[j]; sum -= luMatrixi[j]*sourceSol[j];
} }
@ -171,11 +171,11 @@ void Foam::LUBacksubstitute
List<Type>& sourceSol List<Type>& sourceSol
) )
{ {
label n = luMatrix.n(); label m = luMatrix.m();
label ii = 0; label ii = 0;
for (label i=0; i<n; i++) for (label i=0; i<m; i++)
{ {
Type sum = sourceSol[i]; Type sum = sourceSol[i];
const scalar* __restrict__ luMatrixi = luMatrix[i]; const scalar* __restrict__ luMatrixi = luMatrix[i];
@ -195,12 +195,12 @@ void Foam::LUBacksubstitute
sourceSol[i] = sum/luMatrixi[i]; sourceSol[i] = sum/luMatrixi[i];
} }
for (label i=n-1; i>=0; i--) for (label i=m-1; i>=0; i--)
{ {
Type sum = sourceSol[i]; Type sum = sourceSol[i];
const scalar* __restrict__ luMatrixi = luMatrix[i]; const scalar* __restrict__ luMatrixi = luMatrix[i];
for (label j=i+1; j<n; j++) for (label j=i+1; j<m; j++)
{ {
sum -= luMatrixi[j]*sourceSol[j]; sum -= luMatrixi[j]*sourceSol[j];
} }
@ -217,7 +217,7 @@ void Foam::LUsolve
List<Type>& sourceSol List<Type>& sourceSol
) )
{ {
labelList pivotIndices(matrix.n()); labelList pivotIndices(matrix.m());
LUDecompose(matrix, pivotIndices); LUDecompose(matrix, pivotIndices);
LUBacksubstitute(matrix, pivotIndices, sourceSol); LUBacksubstitute(matrix, pivotIndices, sourceSol);
} }
@ -243,21 +243,21 @@ void Foam::multiply
const Matrix<Form, Type>& B const Matrix<Form, Type>& B
) )
{ {
if (A.m() != B.n()) if (A.n() != B.m())
{ {
FatalErrorInFunction FatalErrorInFunction
<< "A and B must have identical inner dimensions but A.m = " << "A and B must have identical inner dimensions but A.n = "
<< A.m() << " and B.n = " << B.n() << A.n() << " and B.m = " << B.m()
<< abort(FatalError); << abort(FatalError);
} }
ans = Matrix<Form, Type>(A.n(), B.m(), scalar(0)); ans = Matrix<Form, Type>(A.m(), B.n(), scalar(0));
for (label i = 0; i < A.n(); i++) for (label i = 0; i < A.m(); i++)
{ {
for (label j = 0; j < B.m(); j++) for (label j = 0; j < B.n(); j++)
{ {
for (label l = 0; l < B.n(); l++) for (label l = 0; l < B.m(); l++)
{ {
ans[i][j] += A[i][l]*B[l][j]; ans[i][j] += A[i][l]*B[l][j];
} }

4
src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -106,7 +106,7 @@ void Foam::simpleMatrix<Type>::operator=(const simpleMatrix<Type>& m)
<< abort(FatalError); << abort(FatalError);
} }
if (n() != m.n()) if (m() != m.m())
{ {
FatalErrorInFunction FatalErrorInFunction
<< "Different size matrices" << "Different size matrices"

6
src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.C
View File

@ -120,7 +120,7 @@ void Foam::CentredFitSnGradData<Polynomial>::calcFit
} }
// Additional weighting for constant and linear terms // Additional weighting for constant and linear terms
for (label i = 0; i < B.n(); i++) for (label i = 0; i < B.m(); i++)
{ {
B[i][0] *= wts[0]; B[i][0] *= wts[0];
B[i][1] *= wts[0]; B[i][1] *= wts[0];
@ -172,13 +172,13 @@ void Foam::CentredFitSnGradData<Polynomial>::calcFit
wts[0] *= 10; wts[0] *= 10;
wts[1] *= 10; wts[1] *= 10;
for (label j = 0; j < B.m(); j++) for (label j = 0; j < B.n(); j++)
{ {
B[0][j] *= 10; B[0][j] *= 10;
B[1][j] *= 10; B[1][j] *= 10;
} }
for (label i = 0; i < B.n(); i++) for (label i = 0; i < B.m(); i++)
{ {
B[i][0] *= 10; B[i][0] *= 10;
B[i][1] *= 10; B[i][1] *= 10;

8
src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -194,7 +194,7 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit
} }
// Additional weighting for constant and linear terms // Additional weighting for constant and linear terms
for (label i = 0; i < B.n(); i++) for (label i = 0; i < B.m(); i++)
{ {
B[i][0] *= wts[0]; B[i][0] *= wts[0];
B[i][1] *= wts[0]; B[i][1] *= wts[0];
@ -267,13 +267,13 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::calcFit
wts[1] *= 10; wts[1] *= 10;
} }
for (label j = 0; j < B.m(); j++) for (label j = 0; j < B.n(); j++)
{ {
B[0][j] *= 10; B[0][j] *= 10;
B[1][j] *= 10; B[1][j] *= 10;
} }
for (label i = 0; i < B.n(); i++) for (label i = 0; i < B.m(); i++)
{ {
B[i][0] *= 10; B[i][0] *= 10;
B[i][1] *= 10; B[i][1] *= 10;

2
src/thermophysicalModels/radiation/radiationModels/viewFactor/viewFactor.C
View File

@ -232,7 +232,7 @@ void Foam::radiation::viewFactor::initialise()
) )
); );
pivotIndices_.setSize(CLU_().n()); pivotIndices_.setSize(CLU_().m());
} }
} }
} }