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 -------------------------------------------------------------------------------
10 License
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 
35 Foam::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 
387 Foam::scalarRectangularMatrix Foam::SVD::VSinvUt() const
388 {
389  scalarRectangularMatrix VSinvUt;
390  multiply(VSinvUt, V_, inv(S_), U_.T());
391  return VSinvUt;
392 }
393 
394 
395 // ************************************************************************* //
Foam::sqrtSumSqr
Scalar sqrtSumSqr(const Scalar a, const Scalar b)
Definition: Scalar.H:412
s
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))
Definition: gmvOutputSpray.H:25
SVD.H
Foam::Matrix::n
label n() const noexcept
The number of columns.
Definition: MatrixI.H:103
A
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
Foam::dimensioned::T
dimensioned< Type > T() const
Return transpose.
Definition: dimensionedSphericalTensor.C:38
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
scalarList.H
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::SVD::VSinvUt
scalarRectangularMatrix VSinvUt() const
Return the matrix product V S^(-1) U^T (the pseudo inverse)
Definition: SVD.C:387
Foam::constant::universal::h
const dimensionedScalar h
Planck constant.
Definition: setRegionSolidFields.H:33
Foam::inv
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
Definition: dimensionedSphericalTensor.C:73
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Foam::RectangularMatrix< scalar >
g
const uniformDimensionedVectorField & g
Definition: createFluidFields.H:26
Foam::Matrix::m
label m() const noexcept
The number of rows.
Definition: MatrixI.H:96
f
labelList f(nPoints)
Foam::findMax
label findMax(const ListType &input, label start=0)
Foam::List< scalar >
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
k
label k
Boltzmann constant.
Definition: LISASMDCalcMethod2.H:41
x
x
Definition: LISASMDCalcMethod2.H:52
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
scalarMatrices.H
ListOps.H
Various functions to operate on Lists.
Foam::multiply
void multiply(FieldField< Field, Type > &f, const FieldField< Field, Type > &f1, const FieldField< Field, scalar > &f2)
y
scalar y
Definition: LISASMDCalcMethod1.H:14