SVD.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) 2011-2017 OpenFOAM Foundation
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26\*---------------------------------------------------------------------------*/
27
28#include "SVD.H"
29#include "scalarList.H"
30#include "scalarMatrices.H"
31#include "ListOps.H"
32
33// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * //
34
35Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
36:
37 U_(A),
38 V_(A.n(), A.n()),
39 S_(A.n()),
40 converged_(true),
41 nZeros_(0)
42{
43 // SVDcomp to find U_, V_ and S_ - the singular values
44
45 const label Un = U_.n();
46 const label Um = U_.m();
47
48 scalarList rv1(Un);
49 scalar g = 0;
50 scalar scale = 0;
51 scalar s = 0;
52 scalar anorm = 0;
53 label l = 0;
54
55 for (label i=0; i<Un; i++)
56 {
57 l = i + 2;
58 rv1[i] = scale*g;
59 g = s = scale = 0;
60
61 if (i < Um)
62 {
63 for (label k=i; k<Um; k++)
64 {
65 scale += mag(U_(k, i));
66 }
67
68 if (scale != 0)
69 {
70 for (label k=i; k<Um; k++)
71 {
72 U_(k, i) /= scale;
73 s += U_(k, i)*U_(k, i);
74 }
75
76 scalar f = U_(i, i);
77 g = -sign(sqrt(s), f);
78 scalar h = f*g - s;
79 U_(i, i) = f - g;
80
81 for (label j=l-1; j<Un; j++)
82 {
83 s = 0;
84 for (label k=i; k<Um; k++)
85 {
86 s += U_(k, i)*U_(k, j);
87 }
88
89 f = s/h;
90 for (label k=i; k<A.m(); k++)
91 {
92 U_(k, j) += f*U_(k, i);
93 }
94 }
95
96 for (label k=i; k<Um; k++)
97 {
98 U_(k, i) *= scale;
99 }
100 }
101 }
102
103 S_[i] = scale*g;
104
105 g = s = scale = 0;
106
107 if (i+1 <= Um && i != Un)
108 {
109 for (label k=l-1; k<Un; k++)
110 {
111 scale += mag(U_(i, k));
112 }
113
114 if (scale != 0)
115 {
116 for (label k=l-1; k<Un; k++)
117 {
118 U_(i, k) /= scale;
119 s += U_(i, k)*U_(i, k);
120 }
121
122 scalar f = U_(i, l-1);
123 g = -sign(sqrt(s), f);
124 scalar h = f*g - s;
125 U_(i, l-1) = f - g;
126
127 for (label k=l-1; k<Un; k++)
128 {
129 rv1[k] = U_(i, k)/h;
130 }
131
132 for (label j=l-1; j<Um; j++)
133 {
134 s = 0;
135 for (label k=l-1; k<Un; k++)
136 {
137 s += U_(j, k)*U_(i, k);
138 }
139
140 for (label k=l-1; k<Un; k++)
141 {
142 U_(j, k) += s*rv1[k];
143 }
144 }
145 for (label k=l-1; k<Un; k++)
146 {
147 U_(i, k) *= scale;
148 }
149 }
150 }
151
152 anorm = max(anorm, mag(S_[i]) + mag(rv1[i]));
153 }
154
155 anorm *= SMALL;
156
157 for (label i=Un-1; i >= 0; i--)
158 {
159 if (i < Un-1)
160 {
161 if (g*U_(i, l) != 0)
162 {
163 for (label j=l; j<Un; j++)
164 {
165 V_(j, i) = (U_(i, j)/U_(i, l))/g;
166 }
167
168 for (label j=l; j<Un; j++)
169 {
170 s = 0;
171 for (label k=l; k<Un; k++)
172 {
173 s += U_(i, k)*V_(k, j);
174 }
175
176 for (label k=l; k<Un; k++)
177 {
178 V_(k, j) += s*V_(k, i);
179 }
180 }
181 }
182
183 for (label j=l; j<Un;j++)
184 {
185 V_(i, j) = V_(j, i) = 0;
186 }
187 }
188
189 V_(i, i) = 1;
190 g = rv1[i];
191 l = i;
192 }
193
194 for (label i=min(Un, Um) - 1; i>=0; i--)
195 {
196 l = i+1;
197 g = S_[i];
198
199 for (label j=l; j<Un; j++)
200 {
201 U_(i, j) = 0;
202 }
203
204 if (g*U_(i, i) != 0)
205 {
206 g = 1.0/g;
207
208 for (label j=l; j<Un; j++)
209 {
210 s = 0;
211 for (label k=l; k<Um; k++)
212 {
213 s += U_(k, i)*U_(k, j);
214 }
215
216 scalar f = (s/U_(i, i))*g;
217
218 for (label k=i; k<Um; k++)
219 {
220 U_(k, j) += f*U_(k, i);
221 }
222 }
223
224 for (label j=i; j<Um; j++)
225 {
226 U_(j, i) *= g;
227 }
228 }
229 else
230 {
231 for (label j=i; j<Um; j++)
232 {
233 U_(j, i) = 0;
234 }
235 }
236
237 ++U_(i, i);
238 }
239
240 for (label k=Un-1; k >= 0; k--)
241 {
242 for (label its = 0; its < 30; its++)
243 {
244 bool flag = true;
245
246 label mn;
247 for (l = k; l >= 0; l--)
248 {
249 mn = l - 1;
250
251 if (l == 0 || mag(rv1[l]) <= anorm)
252 {
253 flag = false;
254 break;
255 }
256
257 if (mag(S_[mn]) <= anorm)
258 {
259 break;
260 }
261 }
262
263 if (flag)
264 {
265 scalar c = 0;
266 s = 1;
267 for (label i=l; i<k+1; i++)
268 {
269 scalar f = s*rv1[i];
270 rv1[i] = c*rv1[i];
271
272 if (mag(f) <= anorm)
273 {
274 break;
275 }
276
277 g = S_[i];
278 scalar h = sqrtSumSqr(f, g);
279 S_[i] = h;
280 h = 1.0/h;
281 c = g*h;
282 s = -f*h;
283
284 for (label j=0; j<Um; j++)
285 {
286 scalar y = U_(j, mn);
287 scalar z = U_(j, i);
288 U_(j, mn) = y*c + z*s;
289 U_(j, i) = z*c - y*s;
290 }
291 }
292 }
293
294 scalar z = S_[k];
295
296 if (l == k)
297 {
298 if (z < 0)
299 {
300 S_[k] = -z;
301 for (label j=0; j<Un; j++)
302 {
303 V_(j, k) = -V_(j, k);
304 }
305 }
306 break;
307 }
308
309 if (its == 29)
310 {
311 converged_ = false;
312 }
313
314 scalar x = S_[l];
315 mn = k-1;
316 scalar y = S_[mn];
317 g = rv1[mn];
318 scalar h = rv1[k];
319 scalar f = ((y - z)*(y + z) + (g - h)*(g + h))/(2*h*y);
320 g = sqrtSumSqr(f, scalar(1));
321 f = ((x - z)*(x + z) + h*((y/(f + sign(g, f))) - h))/x;
322 scalar c = 1;
323 s = 1;
324
325 for (label j=l; j <= mn; j++)
326 {
327 label i = j + 1;
328 g = rv1[i];
329 y = S_[i];
330 h = s*g;
331 g = c*g;
332 scalar z = sqrtSumSqr(f, h);
333 rv1[j] = z;
334 c = f/z;
335 s = h/z;
336 f = x*c + g*s;
337 g = g*c - x*s;
338 h = y*s;
339 y *= c;
340
341 for (label jj = 0; jj < Un; jj++)
342 {
343 x = V_(jj, j);
344 z = V_(jj, i);
345 V_(jj, j) = x*c + z*s;
346 V_(jj, i) = z*c - x*s;
347 }
348
349 z = sqrtSumSqr(f, h);
350 S_[j] = z;
351 if (z)
352 {
353 z = 1.0/z;
354 c = f*z;
355 s = h*z;
356 }
357 f = c*g + s*y;
358 x = c*y - s*g;
359
360 for (label jj=0; jj < Um; jj++)
361 {
362 y = U_(jj, j);
363 z = U_(jj, i);
364 U_(jj, j) = y*c + z*s;
365 U_(jj, i) = z*c - y*s;
366 }
367 }
368 rv1[l] = 0;
369 rv1[k] = f;
370 S_[k] = x;
371 }
372 }
373
374 // Zero singular values that are less than minCondition*maxS
375 const scalar minS = minCondition*S_[findMax(S_)];
376 forAll(S_, i)
377 {
378 if (S_[i] <= minS)
379 {
380 S_[i] = 0;
381 nZeros_++;
382 }
383 }
384}
385
386
388{
390 multiply(VSinvUt, V_, inv(S_), U_.T());
391 return VSinvUt;
392}
393
394
395// ************************************************************************* //
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
scalar y
label k
Various functions to operate on Lists.
label n
const uniformDimensionedVectorField & g
label n() const noexcept
The number of columns.
Definition: MatrixI.H:103
label m() const noexcept
The number of rows.
Definition: MatrixI.H:96
Singular value decomposition of a rectangular matrix.
Definition: SVD.H:54
scalarRectangularMatrix VSinvUt() const
Return the matrix product V S^(-1) U^T (the pseudo inverse)
Definition: SVD.C:387
dimensioned< Type > T() const
Return transpose.
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
void multiply(FieldField< Field, Type > &f, const FieldField< Field, Type > &f1, const FieldField< Field, scalar > &f2)
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
label findMax(const ListType &input, label start=0)
Scalar sqrtSumSqr(const Scalar a, const Scalar b)
Definition: Scalar.H:417
labelList f(nPoints)
volScalarField & h
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333