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