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