src/OpenFOAM/matrices/scalarMatrix/scalarMatrix.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2008 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 "scalarMatrix.H"
  28
  29 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  30
  31 Foam::scalarMatrix::scalarMatrix()
  32 {}
  33
  34
  35 Foam::scalarMatrix::scalarMatrix(const label mSize)
  36 :
  37     Matrix<scalar>(mSize, mSize, 0.0)
  38 {}
  39
  40
  41 Foam::scalarMatrix::scalarMatrix(const Matrix<scalar>& matrix)
  42 :
  43     Matrix<scalar>(matrix)
  44 {}
  45
  46
  47 Foam::scalarMatrix::scalarMatrix(Istream& is)
  48 :
  49     Matrix<scalar>(is)
  50 {}
  51
  52
  53 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  54
  55 void Foam::scalarMatrix::LUDecompose
  56 (
  57     Matrix<scalar>& matrix,
  58     labelList& pivotIndices
  59 )
  60 {
  61     label n = matrix.n();
  62     scalar vv[n];
  63
  64     for (register label i=0; i<n; i++)
  65     {
  66         scalar largestCoeff = 0.0;
  67         scalar temp;
  68         const scalar* __restrict__ matrixi = matrix[i];
  69
  70         for (register label j=0; j<n; j++)
  71         {
  72             if ((temp = mag(matrixi[j])) > largestCoeff)
  73             {
  74                 largestCoeff = temp;
  75             }
  76         }
  77
  78         if (largestCoeff == 0.0)
  79         {
  80             FatalErrorIn
  81             (
  82                 "scalarMatrix::LUdecompose"
  83                 "(Matrix<scalar>& matrix, labelList& rowIndices)"
  84             )   << "Singular matrix" << exit(FatalError);
  85         }
  86
  87         vv[i] = 1.0/largestCoeff;
  88     }
  89
  90     for (register label j=0; j<n; j++)
  91     {
  92         scalar* __restrict__ matrixj = matrix[j];
  93
  94         for (register label i=0; i<j; i++)
  95         {
  96             scalar* __restrict__ matrixi = matrix[i];
  97
  98             scalar sum = matrixi[j];
  99             for (register label k=0; k<i; k++)
 100             {
 101                 sum -= matrixi[k]*matrix[k][j];
 102             }
 103             matrixi[j] = sum;
 104         }
 105
 106         label iMax = 0;
 107
 108         scalar largestCoeff = 0.0;
 109         for (register label i=j; i<n; i++)
 110         {
 111             scalar* __restrict__ matrixi = matrix[i];
 112             scalar sum = matrixi[j];
 113
 114             for (register label k=0; k<j; k++)
 115             {
 116                 sum -= matrixi[k]*matrix[k][j];
 117             }
 118
 119             matrixi[j] = sum;
 120
 121             scalar temp;
 122             if ((temp = vv[i]*mag(sum)) >= largestCoeff)
 123             {
 124                 largestCoeff = temp;
 125                 iMax = i;
 126             }
 127         }
 128
 129         pivotIndices[j] = iMax;
 130
 131         if (j != iMax)
 132         {
 133             scalar* __restrict__ matrixiMax = matrix[iMax];
 134
 135             for (register label k=0; k<n; k++)
 136             {
 137                 Swap(matrixj[k], matrixiMax[k]);
 138             }
 139
 140             vv[iMax] = vv[j];
 141         }
 142
 143         if (matrixj[j] == 0.0)
 144         {
 145             matrixj[j] = SMALL;
 146         }
 147
 148         if (j != n-1)
 149         {
 150             scalar rDiag = 1.0/matrixj[j];
 151
 152             for (register label i=j+1; i<n; i++)
 153             {
 154                 matrix[i][j] *= rDiag;
 155             }
 156         }
 157     }
 158 }
 159
 160
 161 // ************************************************************************* //