src/foam/matrices/scalarMatrices/scalarMatrices.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend 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 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "scalarMatrices.H"
  27 #include "SVD.H"
  28
  29 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
  30
  31 void Foam::multiply
  32 (
  33     scalarRectangularMatrix& ans,         // value changed in return
  34     const scalarRectangularMatrix& A,
  35     const scalarRectangularMatrix& B
  36 )
  37 {
  38     if (A.m() != B.n())
  39     {
  40         FatalErrorIn
  41         (
  42             "multiply("
  43             "scalarRectangularMatrix& answer "
  44             "const scalarRectangularMatrix& A, "
  45             "const scalarRectangularMatrix& B)"
  46         )   << "A and B must have identical inner dimensions but A.m = "
  47             << A.m() << " and B.n = " << B.n()
  48             << abort(FatalError);
  49     }
  50
  51     ans = scalarRectangularMatrix(A.n(), B.m(), scalar(0));
  52
  53     for(register label i = 0; i < A.n(); i++)
  54     {
  55         for(register label j = 0; j < B.m(); j++)
  56         {
  57             for(register label l = 0; l < B.n(); l++)
  58             {
  59                 ans[i][j] += A[i][l]*B[l][j];
  60             }
  61         }
  62     }
  63 }
  64
  65
  66 void Foam::multiply
  67 (
  68     scalarRectangularMatrix& ans,         // value changed in return
  69     const scalarRectangularMatrix& A,
  70     const scalarRectangularMatrix& B,
  71     const scalarRectangularMatrix& C
  72 )
  73 {
  74     if (A.m() != B.n())
  75     {
  76         FatalErrorIn
  77         (
  78             "multiply("
  79             "const scalarRectangularMatrix& A, "
  80             "const scalarRectangularMatrix& B, "
  81             "const scalarRectangularMatrix& C, "
  82             "scalarRectangularMatrix& answer)"
  83         )   << "A and B must have identical inner dimensions but A.m = "
  84             << A.m() << " and B.n = " << B.n()
  85             << abort(FatalError);
  86     }
  87
  88     if (B.m() != C.n())
  89     {
  90         FatalErrorIn
  91         (
  92             "multiply("
  93             "const scalarRectangularMatrix& A, "
  94             "const scalarRectangularMatrix& B, "
  95             "const scalarRectangularMatrix& C, "
  96             "scalarRectangularMatrix& answer)"
  97         )   << "B and C must have identical inner dimensions but B.m = "
  98             << B.m() << " and C.n = " << C.n()
  99             << abort(FatalError);
 100     }
 101
 102     ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0));
 103
 104     for(register label i = 0; i < A.n(); i++)
 105     {
 106         for(register label g = 0; g < C.m(); g++)
 107         {
 108             for(register label l = 0; l < C.n(); l++)
 109             {
 110                 scalar ab = 0;
 111                 for(register label j = 0; j < A.m(); j++)
 112                 {
 113                     ab += A[i][j]*B[j][l];
 114                 }
 115                 ans[i][g] += C[l][g] * ab;
 116             }
 117         }
 118     }
 119 }
 120
 121
 122 void Foam::multiply
 123 (
 124     scalarRectangularMatrix& ans,         // value changed in return
 125     const scalarRectangularMatrix& A,
 126     const DiagonalMatrix<scalar>& B,
 127     const scalarRectangularMatrix& C
 128 )
 129 {
 130     if (A.m() != B.size())
 131     {
 132         FatalErrorIn
 133         (
 134             "multiply("
 135             "const scalarRectangularMatrix& A, "
 136             "const DiagonalMatrix<scalar>& B, "
 137             "const scalarRectangularMatrix& C, "
 138             "scalarRectangularMatrix& answer)"
 139         )   << "A and B must have identical inner dimensions but A.m = "
 140             << A.m() << " and B.n = " << B.size()
 141             << abort(FatalError);
 142     }
 143
 144     if (B.size() != C.n())
 145     {
 146         FatalErrorIn
 147         (
 148             "multiply("
 149             "const scalarRectangularMatrix& A, "
 150             "const DiagonalMatrix<scalar>& B, "
 151             "const scalarRectangularMatrix& C, "
 152             "scalarRectangularMatrix& answer)"
 153         )   << "B and C must have identical inner dimensions but B.m = "
 154             << B.size() << " and C.n = " << C.n()
 155             << abort(FatalError);
 156     }
 157
 158     ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0));
 159
 160     for(register label i = 0; i < A.n(); i++)
 161     {
 162         for(register label g = 0; g < C.m(); g++)
 163         {
 164             for(register label l = 0; l < C.n(); l++)
 165             {
 166                 ans[i][g] += C[l][g] * A[i][l]*B[l];
 167             }
 168         }
 169     }
 170 }
 171
 172
 173 Foam::RectangularMatrix<Foam::scalar> Foam::SVDinv
 174 (
 175     const scalarRectangularMatrix& A,
 176     scalar minCondition
 177 )
 178 {
 179     SVD svd(A, minCondition);
 180     return svd.VSinvUt();
 181 }
 182
 183
 184 // ************************************************************************* //