src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.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 "scalarMatrices.H"
  28 #include "SVD.H"
  29
  30 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  31
  32 void Foam::LUDecompose
  33 (
  34     scalarSquareMatrix& matrix,
  35     labelList& pivotIndices
  36 )
  37 {
  38     label n = matrix.n();
  39     scalar vv[n];
  40
  41     for (register label i=0; i<n; i++)
  42     {
  43         scalar largestCoeff = 0.0;
  44         scalar temp;
  45         const scalar* __restrict__ matrixi = matrix[i];
  46
  47         for (register label j=0; j<n; j++)
  48         {
  49             if ((temp = mag(matrixi[j])) > largestCoeff)
  50             {
  51                 largestCoeff = temp;
  52             }
  53         }
  54
  55         if (largestCoeff == 0.0)
  56         {
  57             FatalErrorIn
  58             (
  59                 "LUdecompose"
  60                 "(scalarSquareMatrix& matrix, labelList& rowIndices)"
  61             )   << "Singular matrix" << exit(FatalError);
  62         }
  63
  64         vv[i] = 1.0/largestCoeff;
  65     }
  66
  67     for (register label j=0; j<n; j++)
  68     {
  69         scalar* __restrict__ matrixj = matrix[j];
  70
  71         for (register label i=0; i<j; i++)
  72         {
  73             scalar* __restrict__ matrixi = matrix[i];
  74
  75             scalar sum = matrixi[j];
  76             for (register label k=0; k<i; k++)
  77             {
  78                 sum -= matrixi[k]*matrix[k][j];
  79             }
  80             matrixi[j] = sum;
  81         }
  82
  83         label iMax = 0;
  84
  85         scalar largestCoeff = 0.0;
  86         for (register label i=j; i<n; i++)
  87         {
  88             scalar* __restrict__ matrixi = matrix[i];
  89             scalar sum = matrixi[j];
  90
  91             for (register label k=0; k<j; k++)
  92             {
  93                 sum -= matrixi[k]*matrix[k][j];
  94             }
  95
  96             matrixi[j] = sum;
  97
  98             scalar temp;
  99             if ((temp = vv[i]*mag(sum)) >= largestCoeff)
 100             {
 101                 largestCoeff = temp;
 102                 iMax = i;
 103             }
 104         }
 105
 106         pivotIndices[j] = iMax;
 107
 108         if (j != iMax)
 109         {
 110             scalar* __restrict__ matrixiMax = matrix[iMax];
 111
 112             for (register label k=0; k<n; k++)
 113             {
 114                 Swap(matrixj[k], matrixiMax[k]);
 115             }
 116
 117             vv[iMax] = vv[j];
 118         }
 119
 120         if (matrixj[j] == 0.0)
 121         {
 122             matrixj[j] = SMALL;
 123         }
 124
 125         if (j != n-1)
 126         {
 127             scalar rDiag = 1.0/matrixj[j];
 128
 129             for (register label i=j+1; i<n; i++)
 130             {
 131                 matrix[i][j] *= rDiag;
 132             }
 133         }
 134     }
 135 }
 136
 137
 138 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 139
 140 void Foam::multiply
 141 (
 142     scalarRectangularMatrix& ans,         // value changed in return
 143     const scalarRectangularMatrix& A,
 144     const scalarRectangularMatrix& B
 145 )
 146 {
 147     if (A.m() != B.n())
 148     {
 149         FatalErrorIn
 150         (
 151             "multiply("
 152             "scalarRectangularMatrix& answer "
 153             "const scalarRectangularMatrix& A, "
 154             "const scalarRectangularMatrix& B)"
 155         )   << "A and B must have identical inner dimensions but A.m = "
 156             << A.m() << " and B.n = " << B.n()
 157             << abort(FatalError);
 158     }
 159
 160     ans = scalarRectangularMatrix(A.n(), B.m(), scalar(0));
 161
 162     for(register label i = 0; i < A.n(); i++)
 163     {
 164         for(register label j = 0; j < B.m(); j++)
 165         {
 166             for(register label l = 0; l < B.n(); l++)
 167             {
 168                 ans[i][j] += A[i][l]*B[l][j];
 169             }
 170         }
 171     }
 172 }
 173
 174
 175 void Foam::multiply
 176 (
 177     scalarRectangularMatrix& ans,         // value changed in return
 178     const scalarRectangularMatrix& A,
 179     const scalarRectangularMatrix& B,
 180     const scalarRectangularMatrix& C
 181 )
 182 {
 183     if (A.m() != B.n())
 184     {
 185         FatalErrorIn
 186         (
 187             "multiply("
 188             "const scalarRectangularMatrix& A, "
 189             "const scalarRectangularMatrix& B, "
 190             "const scalarRectangularMatrix& C, "
 191             "scalarRectangularMatrix& answer)"
 192         )   << "A and B must have identical inner dimensions but A.m = "
 193             << A.m() << " and B.n = " << B.n()
 194             << abort(FatalError);
 195     }
 196
 197     if (B.m() != C.n())
 198     {
 199         FatalErrorIn
 200         (
 201             "multiply("
 202             "const scalarRectangularMatrix& A, "
 203             "const scalarRectangularMatrix& B, "
 204             "const scalarRectangularMatrix& C, "
 205             "scalarRectangularMatrix& answer)"
 206         )   << "B and C must have identical inner dimensions but B.m = "
 207             << B.m() << " and C.n = " << C.n()
 208             << abort(FatalError);
 209     }
 210
 211     ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0));
 212
 213     for(register label i = 0; i < A.n(); i++)
 214     {
 215         for(register label g = 0; g < C.m(); g++)
 216         {
 217             for(register label l = 0; l < C.n(); l++)
 218             {
 219                 scalar ab = 0;
 220                 for(register label j = 0; j < A.m(); j++)
 221                 {
 222                     ab += A[i][j]*B[j][l];
 223                 }
 224                 ans[i][g] += C[l][g] * ab;
 225             }
 226         }
 227     }
 228 }
 229
 230
 231 void Foam::multiply
 232 (
 233     scalarRectangularMatrix& ans,         // value changed in return
 234     const scalarRectangularMatrix& A,
 235     const DiagonalMatrix<scalar>& B,
 236     const scalarRectangularMatrix& C
 237 )
 238 {
 239     if (A.m() != B.size())
 240     {
 241         FatalErrorIn
 242         (
 243             "multiply("
 244             "const scalarRectangularMatrix& A, "
 245             "const DiagonalMatrix<scalar>& B, "
 246             "const scalarRectangularMatrix& C, "
 247             "scalarRectangularMatrix& answer)"
 248         )   << "A and B must have identical inner dimensions but A.m = "
 249             << A.m() << " and B.n = " << B.size()
 250             << abort(FatalError);
 251     }
 252
 253     if (B.size() != C.n())
 254     {
 255         FatalErrorIn
 256         (
 257             "multiply("
 258             "const scalarRectangularMatrix& A, "
 259             "const DiagonalMatrix<scalar>& B, "
 260             "const scalarRectangularMatrix& C, "
 261             "scalarRectangularMatrix& answer)"
 262         )   << "B and C must have identical inner dimensions but B.m = "
 263             << B.size() << " and C.n = " << C.n()
 264             << abort(FatalError);
 265     }
 266
 267     ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0));
 268
 269     for(register label i = 0; i < A.n(); i++)
 270     {
 271         for(register label g = 0; g < C.m(); g++)
 272         {
 273             for(register label l = 0; l < C.n(); l++)
 274             {
 275                 ans[i][g] += C[l][g] * A[i][l]*B[l];
 276             }
 277         }
 278     }
 279 }
 280
 281
 282 Foam::RectangularMatrix<Foam::scalar> Foam::SVDinv
 283 (
 284     const scalarRectangularMatrix& A,
 285     scalar minCondition
 286 )
 287 {
 288     SVD svd(A, minCondition);
 289     return svd.VSinvUt();
 290 }
 291
 292
 293 // ************************************************************************* //