SVD.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
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 "SVD.H"
27 #include "scalarList.H"
28 #include "scalarMatrices.H"
29 #include "ListOps.H"
30 
31 // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * //
32 
33 Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
34 :
35  U_(A),
36  V_(A.n(), A.n()),
37  S_(A.n()),
38  converged_(true),
39  nZeros_(0)
40 {
41  // SVDcomp to find U_, V_ and S_ - the singular values
42 
43  const label Un = U_.n();
44  const label Um = U_.m();
45 
46  scalarList rv1(Un);
47  scalar g = 0;
48  scalar scale = 0;
49  scalar s = 0;
50  scalar anorm = 0;
51  label l = 0;
52 
53  for (label i=0; i<Un; i++)
54  {
55  l = i + 2;
56  rv1[i] = scale*g;
57  g = s = scale = 0;
58 
59  if (i < Um)
60  {
61  for (label k=i; k<Um; k++)
62  {
63  scale += mag(U_(k, i));
64  }
65 
66  if (scale != 0)
67  {
68  for (label k=i; k<Um; k++)
69  {
70  U_(k, i) /= scale;
71  s += U_(k, i)*U_(k, i);
72  }
73 
74  scalar f = U_(i, i);
75  g = -sign(sqrt(s), f);
76  scalar h = f*g - s;
77  U_(i, i) = f - g;
78 
79  for (label j=l-1; j<Un; j++)
80  {
81  s = 0;
82  for (label k=i; k<Um; k++)
83  {
84  s += U_(k, i)*U_(k, j);
85  }
86 
87  f = s/h;
88  for (label k=i; k<A.m(); k++)
89  {
90  U_(k, j) += f*U_(k, i);
91  }
92  }
93 
94  for (label k=i; k<Um; k++)
95  {
96  U_(k, i) *= scale;
97  }
98  }
99  }
100 
101  S_[i] = scale*g;
102 
103  g = s = scale = 0;
104 
105  if (i+1 <= Um && i != Un)
106  {
107  for (label k=l-1; k<Un; k++)
108  {
109  scale += mag(U_(i, k));
110  }
111 
112  if (scale != 0)
113  {
114  for (label k=l-1; k<Un; k++)
115  {
116  U_(i, k) /= scale;
117  s += U_(i, k)*U_(i, k);
118  }
119 
120  scalar f = U_(i, l-1);
121  g = -sign(sqrt(s), f);
122  scalar h = f*g - s;
123  U_(i, l-1) = f - g;
124 
125  for (label k=l-1; k<Un; k++)
126  {
127  rv1[k] = U_(i, k)/h;
128  }
129 
130  for (label j=l-1; j<Um; j++)
131  {
132  s = 0;
133  for (label k=l-1; k<Un; k++)
134  {
135  s += U_(j, k)*U_(i, k);
136  }
137 
138  for (label k=l-1; k<Un; k++)
139  {
140  U_(j, k) += s*rv1[k];
141  }
142  }
143  for (label k=l-1; k<Un; k++)
144  {
145  U_(i, k) *= scale;
146  }
147  }
148  }
149 
150  anorm = max(anorm, mag(S_[i]) + mag(rv1[i]));
151  }
152 
153  anorm *= small;
154 
155  for (label i=Un-1; i >= 0; i--)
156  {
157  if (i < Un-1)
158  {
159  if (g*U_(i, l) != 0)
160  {
161  for (label j=l; j<Un; j++)
162  {
163  V_(j, i) = (U_(i, j)/U_(i, l))/g;
164  }
165 
166  for (label j=l; j<Un; j++)
167  {
168  s = 0;
169  for (label k=l; k<Un; k++)
170  {
171  s += U_(i, k)*V_(k, j);
172  }
173 
174  for (label k=l; k<Un; k++)
175  {
176  V_(k, j) += s*V_(k, i);
177  }
178  }
179  }
180 
181  for (label j=l; j<Un;j++)
182  {
183  V_(i, j) = V_(j, i) = 0;
184  }
185  }
186 
187  V_(i, i) = 1;
188  g = rv1[i];
189  l = i;
190  }
191 
192  for (label i=min(Un, Um) - 1; i>=0; i--)
193  {
194  l = i+1;
195  g = S_[i];
196 
197  for (label j=l; j<Un; j++)
198  {
199  U_(i, j) = 0;
200  }
201 
202  if (g*U_(i, i) != 0)
203  {
204  g = 1.0/g;
205 
206  for (label j=l; j<Un; j++)
207  {
208  s = 0;
209  for (label k=l; k<Um; k++)
210  {
211  s += U_(k, i)*U_(k, j);
212  }
213 
214  scalar f = (s/U_(i, i))*g;
215 
216  for (label k=i; k<Um; k++)
217  {
218  U_(k, j) += f*U_(k, i);
219  }
220  }
221 
222  for (label j=i; j<Um; j++)
223  {
224  U_(j, i) *= g;
225  }
226  }
227  else
228  {
229  for (label j=i; j<Um; j++)
230  {
231  U_(j, i) = 0;
232  }
233  }
234 
235  ++U_(i, i);
236  }
237 
238  for (label k=Un-1; k >= 0; k--)
239  {
240  for (label its = 0; its < 30; its++)
241  {
242  bool flag = true;
243 
244  label mn;
245  for (l = k; l >= 0; l--)
246  {
247  mn = l - 1;
248 
249  if (l == 0 || mag(rv1[l]) <= anorm)
250  {
251  flag = false;
252  break;
253  }
254 
255  if (mag(S_[mn]) <= anorm)
256  {
257  break;
258  }
259  }
260 
261  if (flag)
262  {
263  scalar c = 0;
264  s = 1;
265  for (label i=l; i<k+1; i++)
266  {
267  scalar f = s*rv1[i];
268  rv1[i] = c*rv1[i];
269 
270  if (mag(f) <= anorm)
271  {
272  break;
273  }
274 
275  g = S_[i];
276  scalar h = sqrtSumSqr(f, g);
277  S_[i] = h;
278  h = 1.0/h;
279  c = g*h;
280  s = -f*h;
281 
282  for (label j=0; j<Um; j++)
283  {
284  scalar y = U_(j, mn);
285  scalar z = U_(j, i);
286  U_(j, mn) = y*c + z*s;
287  U_(j, i) = z*c - y*s;
288  }
289  }
290  }
291 
292  scalar z = S_[k];
293 
294  if (l == k)
295  {
296  if (z < 0)
297  {
298  S_[k] = -z;
299  for (label j=0; j<Un; j++)
300  {
301  V_(j, k) = -V_(j, k);
302  }
303  }
304  break;
305  }
306 
307  if (its == 29)
308  {
309  converged_ = false;
310  }
311 
312  scalar x = S_[l];
313  mn = k-1;
314  scalar y = S_[mn];
315  g = rv1[mn];
316  scalar h = rv1[k];
317  scalar f = ((y - z)*(y + z) + (g - h)*(g + h))/(2*h*y);
318  g = sqrtSumSqr(f, scalar(1));
319  f = ((x - z)*(x + z) + h*((y/(f + sign(g, f))) - h))/x;
320  scalar c = 1;
321  s = 1;
322 
323  for (label j=l; j <= mn; j++)
324  {
325  label i = j + 1;
326  g = rv1[i];
327  y = S_[i];
328  h = s*g;
329  g = c*g;
330  scalar z = sqrtSumSqr(f, h);
331  rv1[j] = z;
332  c = f/z;
333  s = h/z;
334  f = x*c + g*s;
335  g = g*c - x*s;
336  h = y*s;
337  y *= c;
338 
339  for (label jj = 0; jj < Un; jj++)
340  {
341  x = V_(jj, j);
342  z = V_(jj, i);
343  V_(jj, j) = x*c + z*s;
344  V_(jj, i) = z*c - x*s;
345  }
346 
347  z = sqrtSumSqr(f, h);
348  S_[j] = z;
349  if (z)
350  {
351  z = 1.0/z;
352  c = f*z;
353  s = h*z;
354  }
355  f = c*g + s*y;
356  x = c*y - s*g;
357 
358  for (label jj=0; jj < Um; jj++)
359  {
360  y = U_(jj, j);
361  z = U_(jj, i);
362  U_(jj, j) = y*c + z*s;
363  U_(jj, i) = z*c - y*s;
364  }
365  }
366  rv1[l] = 0;
367  rv1[k] = f;
368  S_[k] = x;
369  }
370  }
371 
372  // Zero singular values that are less than minCondition*maxS
373  const scalar minS = minCondition*S_[findMax(S_)];
374  forAll(S_, i)
375  {
376  if (S_[i] <= minS)
377  {
378  S_[i] = 0;
379  nZeros_++;
380  }
381  }
382 }
383 
384 
385 Foam::scalarRectangularMatrix Foam::SVD::VSinvUt() const
386 {
387  scalarRectangularMatrix VSinvUt;
388  multiply(VSinvUt, V_, inv(S_), U_.T());
389  return VSinvUt;
390 }
391 
392 
393 // ************************************************************************* //
Foam::Matrix::n
label n() const
Return the number of columns.
Definition: MatrixI.H:64
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Foam::max
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
Definition: fvMeshStitchersMoving.C:102
Foam::multiply
void multiply(FieldField< Field, Type > &f, const FieldField< Field, Type > &f1, const FieldField< Field, scalar > &f2)
Foam::SVD::SVD
SVD(const scalarRectangularMatrix &A, const scalar minCondition=0)
Construct from a rectangular Matrix.
Definition: SVD.C:33
Foam::Matrix::T
Form T() const
Return the transpose of the matrix.
Definition: Matrix.C:264
Foam::inv
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
Definition: dimensionedSphericalTensor.C:71
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:142
k
label k
Boltzmann constant.
Definition: LISASMDCalcMethod2.H:41
Foam::constant::universal::h
const dimensionedScalar h
Planck constant.
Foam::constant::universal::c
const dimensionedScalar c
Speed of light in a vacuum.
ListOps.H
Various functions to operate on Lists.
y
scalar y
Definition: LISASMDCalcMethod1.H:14
s
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().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))
Foam::sqrtSumSqr
Scalar sqrtSumSqr(const Scalar a, const Scalar b)
Definition: Scalar.H:335
Foam::min
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
Definition: fvMeshStitchersMoving.C:108
Foam::findMax
label findMax(const ListType &, const label start=0)
Find index of max element (and larger than given element).
f
labelList f(nPoints)
Foam::SVD::VSinvUt
scalarRectangularMatrix VSinvUt() const
Return the matrix product V S^(-1) U^T (the pseudo inverse)
Definition: SVD.C:385
Foam::Matrix::m
label m() const
Return the number of rows.
Definition: MatrixI.H:57
Foam::mag
dimensioned< scalar > mag(const dimensioned< Type > &)
n
label n
Definition: TABSMDCalcMethod2.H:31
g
const dimensionedVector & g
Definition: setRegionFluidFields.H:27