SquareMatrix.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2013-2016 OpenFOAM Foundation
9 Copyright (C) 2019-2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "SquareMatrix.H"
30#include "RectangularMatrix.H"
31#include "labelList.H"
32
33// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
34
35template<class Type>
36template<class CompOp>
38(
39 CompOp& compare
40) const
41{
42 List<label> p(this->m());
43 std::iota(p.begin(), p.end(), 0);
44 std::sort
45 (
46 p.begin(),
47 p.end(),
48 [&](label i, label j){ return compare((*this)(i,i), (*this)(j,j)); }
49 );
50
51 return p;
52}
53
54
55template<class Type>
57{
58 #ifdef FULLDEBUG
59 if (this->m() != p.size())
60 {
62 << "Attempt to column-reorder according to an uneven list: " << nl
63 << "SquareMatrix diagonal size = " << this->m() << nl
64 << "Permutation list size = " << p.size() << nl
65 << abort(FatalError);
66 }
67 #endif
68
69 SquareMatrix<Type> reordered(this->sizes());
70
71 label j = 0;
72 for (const label i : p)
73 {
74 reordered.subColumn(j) = this->subColumn(i);
75 ++j;
76 }
77
78 this->transfer(reordered);
79}
80
81
82// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
83
84template<class Type>
85template<class AnyType>
87{
88 Matrix<SquareMatrix<Type>, Type>::operator=(Zero);
89
90 for (label i = 0; i < this->n(); ++i)
91 {
92 this->operator()(i, i) = pTraits<Type>::one;
93 }
94}
95
96
97// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
98
99namespace Foam
100{
101
102// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
103
104//- Return the determinant of the LU decomposed SquareMatrix
105template<class Type>
107(
108 const SquareMatrix<Type>& matrix,
109 const label sign
110)
111{
112 Type diagProduct = pTraits<Type>::one;
113
114 for (label i = 0; i < matrix.m(); ++i)
115 {
116 diagProduct *= matrix(i, i);
117 }
118
119 return sign*diagProduct;
120}
121
122
123//- Return the determinant of SquareMatrix
124template<class Type>
125scalar det(const SquareMatrix<Type>& matrix)
126{
127 SquareMatrix<Type> matrixTmp = matrix;
128
129 labelList pivotIndices(matrix.m());
130 label sign;
131 LUDecompose(matrixTmp, pivotIndices, sign);
132
133 return detDecomposed(matrixTmp, sign);
134}
135
136
137//- Return the SquareMatrix det and the LU decomposition in the original matrix
138template<class Type>
139scalar det(SquareMatrix<Type>& matrix)
140{
141 labelList pivotIndices(matrix.m());
142 label sign;
143 LUDecompose(matrix, pivotIndices, sign);
144
145 return detDecomposed(matrix, sign);
146}
147
148
149//- Return Matrix column-reordered according to
150//- a given permutation labelList
151template<class Type>
153(
154 const SquareMatrix<Type>& mat,
155 const List<label>& p
156)
157{
158 #ifdef FULLDEBUG
159 if (mat.m() != p.size())
160 {
162 << "Attempt to column-reorder according to an uneven list: " << nl
163 << "SquareMatrix diagonal size = " << mat.m() << nl
164 << "Permutation list size = " << p.size() << nl
165 << abort(FatalError);
166 }
167 #endif
168
169 SquareMatrix<Type> reordered(mat.sizes());
170
171 label j = 0;
172 for (const label i : p)
173 {
174 reordered.subColumn(j) = mat.subColumn(i);
175 ++j;
176 }
177
178 return reordered;
179}
180
181
182// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
183
184template<class Type>
186{
187public:
188
190};
191
192
193// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
195} // End namespace Foam
196
197// ************************************************************************* //
label n
Templated identity and dual space identity tensors derived from SphericalTensor.
Definition: Identity.H:52
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:77
A templated (m x n) matrix of objects of <T>. The layout is (mRows x nCols) - row-major order:
Definition: Matrix.H:81
labelPair sizes() const
Return row/column sizes.
Definition: MatrixI.H:117
label m() const noexcept
The number of rows.
Definition: MatrixI.H:96
ConstMatrixBlock< mType > subColumn(const label colIndex, const label rowIndex=0, label len=-1) const
Return const column or column's subset of Matrix.
Definition: MatrixI.H:289
A templated (N x N) square matrix of objects of <Type>, containing N*N elements, derived from Matrix.
Definition: SquareMatrix.H:66
List< label > sortPermutation(CompOp &compare) const
void applyPermutation(const List< label > &p)
Definition: SquareMatrix.C:56
SquareMatrix & operator=(const SquareMatrix &)=default
Copy assignment.
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
A class representing the concept of 1 (one) that can be used to avoid manipulating objects known to b...
Definition: one.H:62
volScalarField & p
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Namespace for OpenFOAM.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
dimensionedScalar sign(const dimensionedScalar &ds)
void LUDecompose(scalarSquareMatrix &matrix, labelList &pivotIndices)
LU decompose the matrix with pivoting.
errorManip< error > abort(error &err)
Definition: errorManip.H:144
DiagonalMatrix< Type > applyPermutation(const DiagonalMatrix< Type > &mat, const List< label > &p)
error FatalError
scalar detDecomposed(const SquareMatrix< Type > &matrix, const label sign)
Return the determinant of the LU decomposed SquareMatrix.
Definition: SquareMatrix.C:107
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53