src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
  31
  32 void Foam::multiply
  33 (
  34     scalarRectangularMatrix& ans,         // value changed in return
  35     const scalarRectangularMatrix& A,
  36     const scalarRectangularMatrix& B
  37 )
  38 {
  39     if (A.m() != B.n())
  40     {
  41         FatalErrorIn
  42         (
  43             "multiply("
  44             "scalarRectangularMatrix& answer "
  45             "const scalarRectangularMatrix& A, "
  46             "const scalarRectangularMatrix& B)"
  47         )   << "A and B must have identical inner dimensions but A.m = "
  48             << A.m() << " and B.n = " << B.n()
  49             << abort(FatalError);
  50     }
  51
  52     ans = scalarRectangularMatrix(A.n(), B.m(), scalar(0));
  53
  54     for(register label i = 0; i < A.n(); i++)
  55     {
  56         for(register label j = 0; j < B.m(); j++)
  57         {
  58             for(register label l = 0; l < B.n(); l++)
  59             {
  60                 ans[i][j] += A[i][l]*B[l][j];
  61             }
  62         }
  63     }
  64 }
  65
  66
  67 void Foam::multiply
  68 (
  69     scalarRectangularMatrix& ans,         // value changed in return
  70     const scalarRectangularMatrix& A,
  71     const scalarRectangularMatrix& B,
  72     const scalarRectangularMatrix& C
  73 )
  74 {
  75     if (A.m() != B.n())
  76     {
  77         FatalErrorIn
  78         (
  79             "multiply("
  80             "const scalarRectangularMatrix& A, "
  81             "const scalarRectangularMatrix& B, "
  82             "const scalarRectangularMatrix& C, "
  83             "scalarRectangularMatrix& answer)"
  84         )   << "A and B must have identical inner dimensions but A.m = "
  85             << A.m() << " and B.n = " << B.n()
  86             << abort(FatalError);
  87     }
  88
  89     if (B.m() != C.n())
  90     {
  91         FatalErrorIn
  92         (
  93             "multiply("
  94             "const scalarRectangularMatrix& A, "
  95             "const scalarRectangularMatrix& B, "
  96             "const scalarRectangularMatrix& C, "
  97             "scalarRectangularMatrix& answer)"
  98         )   << "B and C must have identical inner dimensions but B.m = "
  99             << B.m() << " and C.n = " << C.n()
 100             << abort(FatalError);
 101     }
 102
 103     ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0));
 104
 105     for(register label i = 0; i < A.n(); i++)
 106     {
 107         for(register label g = 0; g < C.m(); g++)
 108         {
 109             for(register label l = 0; l < C.n(); l++)
 110             {
 111                 scalar ab = 0;
 112                 for(register label j = 0; j < A.m(); j++)
 113                 {
 114                     ab += A[i][j]*B[j][l];
 115                 }
 116                 ans[i][g] += C[l][g] * ab;
 117             }
 118         }
 119     }
 120 }
 121
 122
 123 void Foam::multiply
 124 (
 125     scalarRectangularMatrix& ans,         // value changed in return
 126     const scalarRectangularMatrix& A,
 127     const DiagonalMatrix<scalar>& B,
 128     const scalarRectangularMatrix& C
 129 )
 130 {
 131     if (A.m() != B.size())
 132     {
 133         FatalErrorIn
 134         (
 135             "multiply("
 136             "const scalarRectangularMatrix& A, "
 137             "const DiagonalMatrix<scalar>& B, "
 138             "const scalarRectangularMatrix& C, "
 139             "scalarRectangularMatrix& answer)"
 140         )   << "A and B must have identical inner dimensions but A.m = "
 141             << A.m() << " and B.n = " << B.size()
 142             << abort(FatalError);
 143     }
 144
 145     if (B.size() != C.n())
 146     {
 147         FatalErrorIn
 148         (
 149             "multiply("
 150             "const scalarRectangularMatrix& A, "
 151             "const DiagonalMatrix<scalar>& B, "
 152             "const scalarRectangularMatrix& C, "
 153             "scalarRectangularMatrix& answer)"
 154         )   << "B and C must have identical inner dimensions but B.m = "
 155             << B.size() << " and C.n = " << C.n()
 156             << abort(FatalError);
 157     }
 158
 159     ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0));
 160
 161     for(register label i = 0; i < A.n(); i++)
 162     {
 163         for(register label g = 0; g < C.m(); g++)
 164         {
 165             for(register label l = 0; l < C.n(); l++)
 166             {
 167                 ans[i][g] += C[l][g] * A[i][l]*B[l];
 168             }
 169         }
 170     }
 171 }
 172
 173
 174 Foam::RectangularMatrix<Foam::scalar> Foam::SVDinv
 175 (
 176     const scalarRectangularMatrix& A,
 177     scalar minCondition
 178 )
 179 {
 180     SVD svd(A, minCondition);
 181     return svd.VSinvUt();
 182 }
 183
 184
 185 // ************************************************************************* //