FreeFOAM The Cross-Platform CFD Toolkit
scalarMatrices.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 | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "scalarMatrices.H"
27 #include <OpenFOAM/SVD.H>
28 
29 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
30 
32 (
33  scalarSquareMatrix& matrix,
34  labelList& pivotIndices
35 )
36 {
37  label n = matrix.n();
38  scalar vv[n];
39 
40  for (register label i=0; i<n; i++)
41  {
42  scalar largestCoeff = 0.0;
43  scalar temp;
44  const scalar* __restrict__ matrixi = matrix[i];
45 
46  for (register label j=0; j<n; j++)
47  {
48  if ((temp = mag(matrixi[j])) > largestCoeff)
49  {
50  largestCoeff = temp;
51  }
52  }
53 
54  if (largestCoeff == 0.0)
55  {
57  (
58  "LUdecompose"
59  "(scalarSquareMatrix& matrix, labelList& rowIndices)"
60  ) << "Singular matrix" << exit(FatalError);
61  }
62 
63  vv[i] = 1.0/largestCoeff;
64  }
65 
66  for (register label j=0; j<n; j++)
67  {
68  scalar* __restrict__ matrixj = matrix[j];
69 
70  for (register label i=0; i<j; i++)
71  {
72  scalar* __restrict__ matrixi = matrix[i];
73 
74  scalar sum = matrixi[j];
75  for (register label k=0; k<i; k++)
76  {
77  sum -= matrixi[k]*matrix[k][j];
78  }
79  matrixi[j] = sum;
80  }
81 
82  label iMax = 0;
83 
84  scalar largestCoeff = 0.0;
85  for (register label i=j; i<n; i++)
86  {
87  scalar* __restrict__ matrixi = matrix[i];
88  scalar sum = matrixi[j];
89 
90  for (register label k=0; k<j; k++)
91  {
92  sum -= matrixi[k]*matrix[k][j];
93  }
94 
95  matrixi[j] = sum;
96 
97  scalar temp;
98  if ((temp = vv[i]*mag(sum)) >= largestCoeff)
99  {
100  largestCoeff = temp;
101  iMax = i;
102  }
103  }
104 
105  pivotIndices[j] = iMax;
106 
107  if (j != iMax)
108  {
109  scalar* __restrict__ matrixiMax = matrix[iMax];
110 
111  for (register label k=0; k<n; k++)
112  {
113  Swap(matrixj[k], matrixiMax[k]);
114  }
115 
116  vv[iMax] = vv[j];
117  }
118 
119  if (matrixj[j] == 0.0)
120  {
121  matrixj[j] = SMALL;
122  }
123 
124  if (j != n-1)
125  {
126  scalar rDiag = 1.0/matrixj[j];
127 
128  for (register label i=j+1; i<n; i++)
129  {
130  matrix[i][j] *= rDiag;
131  }
132  }
133  }
134 }
135 
136 
137 // * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
138 
139 void Foam::multiply
140 (
141  scalarRectangularMatrix& ans, // value changed in return
142  const scalarRectangularMatrix& A,
143  const scalarRectangularMatrix& B
144 )
145 {
146  if (A.m() != B.n())
147  {
149  (
150  "multiply("
151  "scalarRectangularMatrix& answer "
152  "const scalarRectangularMatrix& A, "
153  "const scalarRectangularMatrix& B)"
154  ) << "A and B must have identical inner dimensions but A.m = "
155  << A.m() << " and B.n = " << B.n()
156  << abort(FatalError);
157  }
158 
159  ans = scalarRectangularMatrix(A.n(), B.m(), scalar(0));
160 
161  for(register label i = 0; i < A.n(); i++)
162  {
163  for(register label j = 0; j < B.m(); j++)
164  {
165  for(register label l = 0; l < B.n(); l++)
166  {
167  ans[i][j] += A[i][l]*B[l][j];
168  }
169  }
170  }
171 }
172 
173 
174 void Foam::multiply
175 (
176  scalarRectangularMatrix& ans, // value changed in return
177  const scalarRectangularMatrix& A,
178  const scalarRectangularMatrix& B,
179  const scalarRectangularMatrix& C
180 )
181 {
182  if (A.m() != B.n())
183  {
185  (
186  "multiply("
187  "const scalarRectangularMatrix& A, "
188  "const scalarRectangularMatrix& B, "
189  "const scalarRectangularMatrix& C, "
190  "scalarRectangularMatrix& answer)"
191  ) << "A and B must have identical inner dimensions but A.m = "
192  << A.m() << " and B.n = " << B.n()
193  << abort(FatalError);
194  }
195 
196  if (B.m() != C.n())
197  {
199  (
200  "multiply("
201  "const scalarRectangularMatrix& A, "
202  "const scalarRectangularMatrix& B, "
203  "const scalarRectangularMatrix& C, "
204  "scalarRectangularMatrix& answer)"
205  ) << "B and C must have identical inner dimensions but B.m = "
206  << B.m() << " and C.n = " << C.n()
207  << abort(FatalError);
208  }
209 
210  ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0));
211 
212  for(register label i = 0; i < A.n(); i++)
213  {
214  for(register label g = 0; g < C.m(); g++)
215  {
216  for(register label l = 0; l < C.n(); l++)
217  {
218  scalar ab = 0;
219  for(register label j = 0; j < A.m(); j++)
220  {
221  ab += A[i][j]*B[j][l];
222  }
223  ans[i][g] += C[l][g] * ab;
224  }
225  }
226  }
227 }
228 
229 
230 void Foam::multiply
231 (
232  scalarRectangularMatrix& ans, // value changed in return
233  const scalarRectangularMatrix& A,
234  const DiagonalMatrix<scalar>& B,
235  const scalarRectangularMatrix& C
236 )
237 {
238  if (A.m() != B.size())
239  {
241  (
242  "multiply("
243  "const scalarRectangularMatrix& A, "
244  "const DiagonalMatrix<scalar>& B, "
245  "const scalarRectangularMatrix& C, "
246  "scalarRectangularMatrix& answer)"
247  ) << "A and B must have identical inner dimensions but A.m = "
248  << A.m() << " and B.n = " << B.size()
249  << abort(FatalError);
250  }
251 
252  if (B.size() != C.n())
253  {
255  (
256  "multiply("
257  "const scalarRectangularMatrix& A, "
258  "const DiagonalMatrix<scalar>& B, "
259  "const scalarRectangularMatrix& C, "
260  "scalarRectangularMatrix& answer)"
261  ) << "B and C must have identical inner dimensions but B.m = "
262  << B.size() << " and C.n = " << C.n()
263  << abort(FatalError);
264  }
265 
266  ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0));
267 
268  for(register label i = 0; i < A.n(); i++)
269  {
270  for(register label g = 0; g < C.m(); g++)
271  {
272  for(register label l = 0; l < C.n(); l++)
273  {
274  ans[i][g] += C[l][g] * A[i][l]*B[l];
275  }
276  }
277  }
278 }
279 
280 
282 (
283  const scalarRectangularMatrix& A,
284  scalar minCondition
285 )
286 {
287  SVD svd(A, minCondition);
288  return svd.VSinvUt();
289 }
290 
291 
292 // ************************ vim: set sw=4 sts=4 et: ************************ //