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