FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

SVD.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "SVD.H"
00027 #include <OpenFOAM/scalarList.H>
00028 #include <OpenFOAM/scalarMatrices.H>
00029 #include <OpenFOAM/ListOps.H>
00030 
00031 // * * * * * * * * * * * * * * * * Constructor  * * * * * * * * * * * * * * //
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     // SVDcomp to find U_, V_ and S_ - the singular values
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     // zero singular values that are less than minCondition*maxS
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             //Info << "Removing " << S_[i] << " < " << minS << endl;
00371             S_[i] = 0;
00372             nZeros_++;
00373         }
00374     }
00375 
00376     // now multiply out to find the pseudo inverse of A, VSinvUt_
00377     multiply(VSinvUt_, V_, inv(S_), U_.T());
00378 
00379     // test SVD
00380     /*scalarRectangularMatrix SVDA(A.n(), A.m());
00381     multiply(SVDA, U_, S_, transpose(V_));
00382     scalar maxDiff = 0;
00383     scalar diff = 0;
00384     for(label i = 0; i < A.n(); i++)
00385     {
00386         for(label j = 0; j < A.m(); j++)
00387         {
00388             diff = mag(A[i][j] - SVDA[i][j]);
00389             if (diff > maxDiff) maxDiff = diff;
00390         }
00391     }
00392     Info << "Maximum discrepancy between A and svd(A) = " << maxDiff << endl;
00393 
00394     if (maxDiff > 4)
00395     {
00396         Info << "singular values " << S_ << endl;
00397     }
00398     */
00399 }
00400 
00401 
00402 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines