Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #include "SVD.H"
00027 #include <OpenFOAM/scalarList.H>
00028 #include <OpenFOAM/scalarMatrices.H>
00029 #include <OpenFOAM/ListOps.H>
00030
00031
00032
00033 Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
00034 :
00035 U_(A),
00036 V_(A.m(), A.m()),
00037 S_(A.m()),
00038 VSinvUt_(A.m(), A.n()),
00039 nZeros_(0)
00040 {
00041
00042
00043 const label Um = U_.m();
00044 const label Un = U_.n();
00045
00046 scalarList rv1(Um);
00047 scalar g = 0;
00048 scalar scale = 0;
00049 scalar s = 0;
00050 scalar anorm = 0;
00051 label l = 0;
00052
00053 for (label i = 0; i < Um; i++)
00054 {
00055 l = i+2;
00056 rv1[i] = scale*g;
00057 g = s = scale = 0;
00058
00059 if (i < Un)
00060 {
00061 for (label k = i; k < Un; k++)
00062 {
00063 scale += mag(U_[k][i]);
00064 }
00065
00066 if (scale != 0)
00067 {
00068 for (label k = i; k < Un; k++)
00069 {
00070 U_[k][i] /= scale;
00071 s += U_[k][i]*U_[k][i];
00072 }
00073
00074 scalar f = U_[i][i];
00075 g = -sign(Foam::sqrt(s), f);
00076 scalar h = f*g - s;
00077 U_[i][i] = f - g;
00078
00079 for (label j = l-1; j < Um; j++)
00080 {
00081 s = 0;
00082 for (label k = i; k < Un; k++)
00083 {
00084 s += U_[k][i]*U_[k][j];
00085 }
00086
00087 f = s/h;
00088 for (label k = i; k < A.n(); k++)
00089 {
00090 U_[k][j] += f*U_[k][i];
00091 }
00092 }
00093
00094 for (label k = i; k < Un; k++)
00095 {
00096 U_[k][i] *= scale;
00097 }
00098 }
00099 }
00100
00101 S_[i] = scale*g;
00102
00103 g = s = scale = 0;
00104
00105 if (i+1 <= Un && i != Um)
00106 {
00107 for (label k = l-1; k < Um; k++)
00108 {
00109 scale += mag(U_[i][k]);
00110 }
00111
00112 if (scale != 0)
00113 {
00114 for (label k=l-1; k < Um; k++)
00115 {
00116 U_[i][k] /= scale;
00117 s += U_[i][k]*U_[i][k];
00118 }
00119
00120 scalar f = U_[i][l-1];
00121 g = -sign(Foam::sqrt(s),f);
00122 scalar h = f*g - s;
00123 U_[i][l-1] = f - g;
00124
00125 for (label k = l-1; k < Um; k++)
00126 {
00127 rv1[k] = U_[i][k]/h;
00128 }
00129
00130 for (label j = l-1; j < Un; j++)
00131 {
00132 s = 0;
00133 for (label k = l-1; k < Um; k++)
00134 {
00135 s += U_[j][k]*U_[i][k];
00136 }
00137
00138 for (label k = l-1; k < Um; k++)
00139 {
00140 U_[j][k] += s*rv1[k];
00141 }
00142 }
00143 for (label k = l-1; k < Um; k++)
00144 {
00145 U_[i][k] *= scale;
00146 }
00147 }
00148 }
00149
00150 anorm = max(anorm, mag(S_[i]) + mag(rv1[i]));
00151 }
00152
00153 for (label i = Um-1; i >= 0; i--)
00154 {
00155 if (i < Um-1)
00156 {
00157 if (g != 0)
00158 {
00159 for (label j = l; j < Um; j++)
00160 {
00161 V_[j][i] = (U_[i][j]/U_[i][l])/g;
00162 }
00163
00164 for (label j=l; j < Um; j++)
00165 {
00166 s = 0;
00167 for (label k = l; k < Um; k++)
00168 {
00169 s += U_[i][k]*V_[k][j];
00170 }
00171
00172 for (label k = l; k < Um; k++)
00173 {
00174 V_[k][j] += s*V_[k][i];
00175 }
00176 }
00177 }
00178
00179 for (label j = l; j < Um;j++)
00180 {
00181 V_[i][j] = V_[j][i] = 0.0;
00182 }
00183 }
00184
00185 V_[i][i] = 1;
00186 g = rv1[i];
00187 l = i;
00188 }
00189
00190 for (label i = min(Um, Un) - 1; i >= 0; i--)
00191 {
00192 l = i+1;
00193 g = S_[i];
00194
00195 for (label j = l; j < Um; j++)
00196 {
00197 U_[i][j] = 0.0;
00198 }
00199
00200 if (g != 0)
00201 {
00202 g = 1.0/g;
00203
00204 for (label j = l; j < Um; j++)
00205 {
00206 s = 0;
00207 for (label k = l; k < Un; k++)
00208 {
00209 s += U_[k][i]*U_[k][j];
00210 }
00211
00212 scalar f = (s/U_[i][i])*g;
00213
00214 for (label k = i; k < Un; k++)
00215 {
00216 U_[k][j] += f*U_[k][i];
00217 }
00218 }
00219
00220 for (label j = i; j < Un; j++)
00221 {
00222 U_[j][i] *= g;
00223 }
00224 }
00225 else
00226 {
00227 for (label j = i; j < Un; j++)
00228 {
00229 U_[j][i] = 0.0;
00230 }
00231 }
00232
00233 ++U_[i][i];
00234 }
00235
00236 for (label k = Um-1; k >= 0; k--)
00237 {
00238 for (label its = 0; its < 35; its++)
00239 {
00240 bool flag = true;
00241
00242 label nm;
00243 for (l = k; l >= 0; l--)
00244 {
00245 nm = l-1;
00246 if (mag(rv1[l]) + anorm == anorm)
00247 {
00248 flag = false;
00249 break;
00250 }
00251 if (mag(S_[nm]) + anorm == anorm) break;
00252 }
00253
00254 if (flag)
00255 {
00256 scalar c = 0.0;
00257 s = 1.0;
00258 for (label i = l-1; i < k+1; i++)
00259 {
00260 scalar f = s*rv1[i];
00261 rv1[i] = c*rv1[i];
00262
00263 if (mag(f) + anorm == anorm) break;
00264
00265 g = S_[i];
00266 scalar h = sqrtSumSqr(f, g);
00267 S_[i] = h;
00268 h = 1.0/h;
00269 c = g*h;
00270 s = -f*h;
00271
00272 for (label j = 0; j < Un; j++)
00273 {
00274 scalar y = U_[j][nm];
00275 scalar z = U_[j][i];
00276 U_[j][nm] = y*c + z*s;
00277 U_[j][i] = z*c - y*s;
00278 }
00279 }
00280 }
00281
00282 scalar z = S_[k];
00283
00284 if (l == k)
00285 {
00286 if (z < 0.0)
00287 {
00288 S_[k] = -z;
00289
00290 for (label j = 0; j < Um; j++) V_[j][k] = -V_[j][k];
00291 }
00292 break;
00293 }
00294 if (its == 34)
00295 {
00296 WarningIn
00297 (
00298 "SVD::SVD"
00299 "(scalarRectangularMatrix& A, const scalar minCondition)"
00300 ) << "no convergence in 35 SVD iterations"
00301 << endl;
00302 }
00303
00304 scalar x = S_[l];
00305 nm = k-1;
00306 scalar y = S_[nm];
00307 g = rv1[nm];
00308 scalar h = rv1[k];
00309 scalar f = ((y - z)*(y + z) + (g - h)*(g + h))/(2.0*h*y);
00310 g = sqrtSumSqr(f, scalar(1));
00311 f = ((x - z)*(x + z) + h*((y/(f + sign(g, f))) - h))/x;
00312 scalar c = 1.0;
00313 s = 1.0;
00314
00315 for (label j = l; j <= nm; j++)
00316 {
00317 label i = j + 1;
00318 g = rv1[i];
00319 y = S_[i];
00320 h = s*g;
00321 g = c*g;
00322 scalar z = sqrtSumSqr(f, h);
00323 rv1[j] = z;
00324 c = f/z;
00325 s = h/z;
00326 f = x*c + g*s;
00327 g = g*c - x*s;
00328 h = y*s;
00329 y *= c;
00330
00331 for (label jj = 0; jj < Um; jj++)
00332 {
00333 x = V_[jj][j];
00334 z = V_[jj][i];
00335 V_[jj][j] = x*c + z*s;
00336 V_[jj][i] = z*c - x*s;
00337 }
00338
00339 z = sqrtSumSqr(f, h);
00340 S_[j] = z;
00341 if (z)
00342 {
00343 z = 1.0/z;
00344 c = f*z;
00345 s = h*z;
00346 }
00347 f = c*g + s*y;
00348 x = c*y - s*g;
00349
00350 for (label jj=0; jj < Un; jj++)
00351 {
00352 y = U_[jj][j];
00353 z = U_[jj][i];
00354 U_[jj][j] = y*c + z*s;
00355 U_[jj][i] = z*c - y*s;
00356 }
00357 }
00358 rv1[l] = 0.0;
00359 rv1[k] = f;
00360 S_[k] = x;
00361 }
00362 }
00363
00364
00365 const scalar minS = minCondition*S_[findMax(S_)];
00366 for (label i = 0; i < S_.size(); i++)
00367 {
00368 if (S_[i] <= minS)
00369 {
00370
00371 S_[i] = 0;
00372 nZeros_++;
00373 }
00374 }
00375
00376
00377 multiply(VSinvUt_, V_, inv(S_), U_.T());
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399 }
00400
00401
00402