src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 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.m(), A.m()),
  37     S_(A.m()),
  38     VSinvUt_(A.m(), A.n()),
  39     nZeros_(0)
  40 {
  41     // SVDcomp to find U_, V_ and S_ - the singular values
  42
  43     const label Um = U_.m();
  44     const label Un = U_.n();
  45
  46     scalarList rv1(Um);
  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 < Um; i++)
  54     {
  55         l = i+2;
  56         rv1[i] = scale*g;
  57         g = s = scale = 0;
  58
  59         if (i < Un)
  60         {
  61             for (label k = i; k < Un; k++)
  62             {
  63                 scale += mag(U_[k][i]);
  64             }
  65
  66             if (scale != 0)
  67             {
  68                 for (label k = i; k < Un; 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(Foam::sqrt(s), f);
  76                 scalar h = f*g - s;
  77                 U_[i][i] = f - g;
  78
  79                 for (label j = l-1; j < Um; j++)
  80                 {
  81                     s = 0;
  82                     for (label k = i; k < Un; k++)
  83                     {
  84                         s += U_[k][i]*U_[k][j];
  85                     }
  86
  87                     f = s/h;
  88                     for (label k = i; k < A.n(); k++)
  89                     {
  90                         U_[k][j] += f*U_[k][i];
  91                     }
  92                 }
  93
  94                 for (label k = i; k < Un; 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 <= Un && i != Um)
 106         {
 107             for (label k = l-1; k < Um; k++)
 108             {
 109                 scale += mag(U_[i][k]);
 110             }
 111
 112             if (scale != 0)
 113             {
 114                 for (label k=l-1; k < Um; 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(Foam::sqrt(s),f);
 122                 scalar h = f*g - s;
 123                 U_[i][l-1] = f - g;
 124
 125                 for (label k = l-1; k < Um; k++)
 126                 {
 127                     rv1[k] = U_[i][k]/h;
 128                 }
 129
 130                 for (label j = l-1; j < Un; j++)
 131                 {
 132                     s = 0;
 133                     for (label k = l-1; k < Um; k++)
 134                     {
 135                         s += U_[j][k]*U_[i][k];
 136                     }
 137
 138                     for (label k = l-1; k < Um; k++)
 139                     {
 140                         U_[j][k] += s*rv1[k];
 141                     }
 142                 }
 143                 for (label k = l-1; k < Um; k++)
 144                 {
 145                     U_[i][k] *= scale;
 146                 }
 147             }
 148         }
 149
 150         anorm = max(anorm, mag(S_[i]) + mag(rv1[i]));
 151     }
 152
 153     for (label i = Um-1; i >= 0; i--)
 154     {
 155         if (i < Um-1)
 156         {
 157             if (g != 0)
 158             {
 159                 for (label j = l; j < Um; j++)
 160                 {
 161                     V_[j][i] = (U_[i][j]/U_[i][l])/g;
 162                 }
 163
 164                 for (label j=l; j < Um; j++)
 165                 {
 166                     s = 0;
 167                     for (label k = l; k < Um; k++)
 168                     {
 169                         s += U_[i][k]*V_[k][j];
 170                     }
 171
 172                     for (label k = l; k < Um; k++)
 173                     {
 174                         V_[k][j] += s*V_[k][i];
 175                     }
 176                 }
 177             }
 178
 179             for (label j = l; j < Um;j++)
 180             {
 181                 V_[i][j] = V_[j][i] = 0.0;
 182             }
 183         }
 184
 185         V_[i][i] = 1;
 186         g = rv1[i];
 187         l = i;
 188     }
 189
 190     for (label i = min(Um, Un) - 1; i >= 0; i--)
 191     {
 192         l = i+1;
 193         g = S_[i];
 194
 195         for (label j = l; j < Um; j++)
 196         {
 197             U_[i][j] = 0.0;
 198         }
 199
 200         if (g != 0)
 201         {
 202             g = 1.0/g;
 203
 204             for (label j = l; j < Um; j++)
 205             {
 206                 s = 0;
 207                 for (label k = l; k < Un; k++)
 208                 {
 209                     s += U_[k][i]*U_[k][j];
 210                 }
 211
 212                 scalar f = (s/U_[i][i])*g;
 213
 214                 for (label k = i; k < Un; k++)
 215                 {
 216                     U_[k][j] += f*U_[k][i];
 217                 }
 218             }
 219
 220             for (label j = i; j < Un; j++)
 221             {
 222                 U_[j][i] *= g;
 223             }
 224         }
 225         else
 226         {
 227             for (label j = i; j < Un; j++)
 228             {
 229                 U_[j][i] = 0.0;
 230             }
 231         }
 232
 233         ++U_[i][i];
 234     }
 235
 236     for (label k = Um-1; k >= 0; k--)
 237     {
 238         for (label its = 0; its < 35; its++)
 239         {
 240             bool flag = true;
 241
 242             label nm;
 243             for (l = k; l >= 0; l--)
 244             {
 245                 nm = l-1;
 246                 if (mag(rv1[l]) + anorm == anorm)
 247                 {
 248                     flag = false;
 249                     break;
 250                 }
 251                 if (mag(S_[nm]) + anorm == anorm) break;
 252             }
 253
 254             if (flag)
 255             {
 256                 scalar c = 0.0;
 257                 s = 1.0;
 258                 for (label i = l-1; i < k+1; i++)
 259                 {
 260                     scalar f = s*rv1[i];
 261                     rv1[i] = c*rv1[i];
 262
 263                     if (mag(f) + anorm == anorm) break;
 264
 265                     g = S_[i];
 266                     scalar h = sqrtSumSqr(f, g);
 267                     S_[i] = h;
 268                     h = 1.0/h;
 269                     c = g*h;
 270                     s = -f*h;
 271
 272                     for (label j = 0; j < Un; j++)
 273                     {
 274                         scalar y = U_[j][nm];
 275                         scalar z = U_[j][i];
 276                         U_[j][nm] = y*c + z*s;
 277                         U_[j][i] = z*c - y*s;
 278                     }
 279                 }
 280             }
 281
 282             scalar z = S_[k];
 283
 284             if (l == k)
 285             {
 286                 if (z < 0.0)
 287                 {
 288                     S_[k] = -z;
 289
 290                     for (label j = 0; j < Um; j++) V_[j][k] = -V_[j][k];
 291                 }
 292                 break;
 293             }
 294             if (its == 34)
 295             {
 296                 WarningIn
 297                 (
 298                     "SVD::SVD"
 299                     "(scalarRectangularMatrix& A, const scalar minCondition)"
 300                 )   << "no convergence in 35 SVD iterations"
 301                     << endl;
 302             }
 303
 304             scalar x = S_[l];
 305             nm = k-1;
 306             scalar y = S_[nm];
 307             g = rv1[nm];
 308             scalar h = rv1[k];
 309             scalar f = ((y - z)*(y + z) + (g - h)*(g + h))/(2.0*h*y);
 310             g = sqrtSumSqr(f, scalar(1));
 311             f = ((x - z)*(x + z) + h*((y/(f + sign(g, f))) - h))/x;
 312             scalar c = 1.0;
 313             s = 1.0;
 314
 315             for (label j = l; j <= nm; j++)
 316             {
 317                 label i = j + 1;
 318                 g = rv1[i];
 319                 y = S_[i];
 320                 h = s*g;
 321                 g = c*g;
 322                 scalar z = sqrtSumSqr(f, h);
 323                 rv1[j] = z;
 324                 c = f/z;
 325                 s = h/z;
 326                 f = x*c + g*s;
 327                 g = g*c - x*s;
 328                 h = y*s;
 329                 y *= c;
 330
 331                 for (label jj = 0; jj < Um; jj++)
 332                 {
 333                     x = V_[jj][j];
 334                     z = V_[jj][i];
 335                     V_[jj][j] = x*c + z*s;
 336                     V_[jj][i] = z*c - x*s;
 337                 }
 338
 339                 z = sqrtSumSqr(f, h);
 340                 S_[j] = z;
 341                 if (z)
 342                 {
 343                     z = 1.0/z;
 344                     c = f*z;
 345                     s = h*z;
 346                 }
 347                 f = c*g + s*y;
 348                 x = c*y - s*g;
 349
 350                 for (label jj=0; jj < Un; jj++)
 351                 {
 352                     y = U_[jj][j];
 353                     z = U_[jj][i];
 354                     U_[jj][j] = y*c + z*s;
 355                     U_[jj][i] = z*c - y*s;
 356                 }
 357             }
 358             rv1[l] = 0.0;
 359             rv1[k] = f;
 360             S_[k] = x;
 361         }
 362     }
 363
 364     // zero singular values that are less than minCondition*maxS
 365     const scalar minS = minCondition*S_[findMax(S_)];
 366     forAll(S_, i)
 367     {
 368         if (S_[i] <= minS)
 369         {
 370             //Info<< "Removing " << S_[i] << " < " << minS << endl;
 371             S_[i] = 0;
 372             nZeros_++;
 373         }
 374     }
 375
 376     // now multiply out to find the pseudo inverse of A, VSinvUt_
 377     multiply(VSinvUt_, V_, inv(S_), U_.T());
 378
 379     // test SVD
 380     /*scalarRectangularMatrix SVDA(A.n(), A.m());
 381     multiply(SVDA, U_, S_, transpose(V_));
 382     scalar maxDiff = 0;
 383     scalar diff = 0;
 384     for (label i = 0; i < A.n(); i++)
 385     {
 386         for (label j = 0; j < A.m(); j++)
 387         {
 388             diff = mag(A[i][j] - SVDA[i][j]);
 389             if (diff > maxDiff) maxDiff = diff;
 390         }
 391     }
 392     Info<< "Maximum discrepancy between A and svd(A) = " << maxDiff << endl;
 393
 394     if (maxDiff > 4)
 395     {
 396         Info<< "singular values " << S_ << endl;
 397     }
 398     */
 399 }
 400
 401
 402 // ************************************************************************* //