src/OpenFOAM/matrices/scalarMatrices/scalarMatricesTemplates.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2010 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
  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 "Swap.H"
  28
  29 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  30
  31 template<class Type>
  32 void Foam::solve
  33 (
  34     scalarSquareMatrix& tmpMatrix,
  35     Field<Type>& sourceSol
  36 )
  37 {
  38     label n = tmpMatrix.n();
  39
  40     // Elimination
  41     for (register label i=0; i<n; i++)
  42     {
  43         label iMax = i;
  44         scalar largestCoeff = mag(tmpMatrix[iMax][i]);
  45
  46         // Swap entries around to find a good pivot
  47         for (register label j=i+1; j<n; j++)
  48         {
  49             if (mag(tmpMatrix[j][i]) > largestCoeff)
  50             {
  51                 iMax = j;
  52                 largestCoeff = mag(tmpMatrix[iMax][i]);
  53             }
  54         }
  55
  56         if (i != iMax)
  57         {
  58             //Info<< "Pivoted on " << i << " " << iMax << endl;
  59
  60             for (register label k=i; k<n; k++)
  61             {
  62                 Swap(tmpMatrix[i][k], tmpMatrix[iMax][k]);
  63             }
  64             Swap(sourceSol[i], sourceSol[iMax]);
  65         }
  66
  67         // Check that the system of equations isn't singular
  68         if (mag(tmpMatrix[i][i]) < 1e-20)
  69         {
  70             FatalErrorIn("solve(scalarSquareMatrix&, Field<Type>& sourceSol)")
  71                 << "Singular Matrix"
  72                 << exit(FatalError);
  73         }
  74
  75         // Reduce to upper triangular form
  76         for (register label j=i+1; j<n; j++)
  77         {
  78             sourceSol[j] -= sourceSol[i]*(tmpMatrix[j][i]/tmpMatrix[i][i]);
  79
  80             for (register label k=n-1; k>=i; k--)
  81             {
  82                 tmpMatrix[j][k] -=
  83                     tmpMatrix[i][k]*tmpMatrix[j][i]/tmpMatrix[i][i];
  84             }
  85         }
  86     }
  87
  88     // Back-substitution
  89     for (register label j=n-1; j>=0; j--)
  90     {
  91         Type ntempvec = pTraits<Type>::zero;
  92
  93         for (register label k=j+1; k<n; k++)
  94         {
  95             ntempvec += tmpMatrix[j][k]*sourceSol[k];
  96         }
  97
  98         sourceSol[j] = (sourceSol[j] - ntempvec)/tmpMatrix[j][j];
  99     }
 100 }
 101
 102
 103 template<class Type>
 104 void Foam::solve
 105 (
 106     Field<Type>& psi,
 107     const scalarSquareMatrix& matrix,
 108     const Field<Type>& source
 109 )
 110 {
 111     scalarSquareMatrix tmpMatrix = matrix;
 112     psi = source;
 113     solve(tmpMatrix, psi);
 114 }
 115
 116
 117 template<class Type>
 118 void Foam::LUBacksubstitute
 119 (
 120     const scalarSquareMatrix& luMatrix,
 121     const labelList& pivotIndices,
 122     Field<Type>& sourceSol
 123 )
 124 {
 125     label n = luMatrix.n();
 126
 127     label ii = 0;
 128
 129     for (register label i=0; i<n; i++)
 130     {
 131         label ip = pivotIndices[i];
 132         Type sum = sourceSol[ip];
 133         sourceSol[ip] = sourceSol[i];
 134         const scalar* __restrict__ luMatrixi = luMatrix[i];
 135
 136         if (ii != 0)
 137         {
 138             for (label j=ii-1; j<i; j++)
 139             {
 140                 sum -= luMatrixi[j]*sourceSol[j];
 141             }
 142         }
 143         else if (sum != pTraits<Type>::zero)
 144         {
 145             ii = i+1;
 146         }
 147
 148         sourceSol[i] = sum;
 149     }
 150
 151     for (register label i=n-1; i>=0; i--)
 152     {
 153         Type sum = sourceSol[i];
 154         const scalar* __restrict__ luMatrixi = luMatrix[i];
 155
 156         for (register label j=i+1; j<n; j++)
 157         {
 158             sum -= luMatrixi[j]*sourceSol[j];
 159         }
 160
 161         sourceSol[i] = sum/luMatrixi[i];
 162     }
 163 }
 164
 165
 166 template<class Type>
 167 void Foam::LUsolve
 168 (
 169     scalarSquareMatrix& matrix,
 170     Field<Type>& sourceSol
 171 )
 172 {
 173     labelList pivotIndices(matrix.n());
 174     LUDecompose(matrix, pivotIndices);
 175     LUBacksubstitute(matrix, pivotIndices, sourceSol);
 176 }
 177
 178
 179 // ************************************************************************* //