FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

scalarMatrices.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "scalarMatrices.H"
00027 #include <OpenFOAM/SVD.H>
00028 
00029 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00030 
00031 void Foam::LUDecompose
00032 (
00033     scalarSquareMatrix& matrix,
00034     labelList& pivotIndices
00035 )
00036 {
00037     label n = matrix.n();
00038     scalar vv[n];
00039 
00040     for (register label i=0; i<n; i++)
00041     {
00042         scalar largestCoeff = 0.0;
00043         scalar temp;
00044         const scalar* __restrict__ matrixi = matrix[i];
00045 
00046         for (register label j=0; j<n; j++)
00047         {
00048             if ((temp = mag(matrixi[j])) > largestCoeff)
00049             {
00050                 largestCoeff = temp;
00051             }
00052         }
00053 
00054         if (largestCoeff == 0.0)
00055         {
00056             FatalErrorIn
00057             (
00058                 "LUdecompose"
00059                 "(scalarSquareMatrix& matrix, labelList& rowIndices)"
00060             )   << "Singular matrix" << exit(FatalError);
00061         }
00062 
00063         vv[i] = 1.0/largestCoeff;
00064     }
00065 
00066     for (register label j=0; j<n; j++)
00067     {
00068         scalar* __restrict__ matrixj = matrix[j];
00069 
00070         for (register label i=0; i<j; i++)
00071         {
00072             scalar* __restrict__ matrixi = matrix[i];
00073 
00074             scalar sum = matrixi[j];
00075             for (register label k=0; k<i; k++)
00076             {
00077                 sum -= matrixi[k]*matrix[k][j];
00078             }
00079             matrixi[j] = sum;
00080         }
00081 
00082         label iMax = 0;
00083 
00084         scalar largestCoeff = 0.0;
00085         for (register label i=j; i<n; i++)
00086         {
00087             scalar* __restrict__ matrixi = matrix[i];
00088             scalar sum = matrixi[j];
00089 
00090             for (register label k=0; k<j; k++)
00091             {
00092                 sum -= matrixi[k]*matrix[k][j];
00093             }
00094 
00095             matrixi[j] = sum;
00096 
00097             scalar temp;
00098             if ((temp = vv[i]*mag(sum)) >= largestCoeff)
00099             {
00100                 largestCoeff = temp;
00101                 iMax = i;
00102             }
00103         }
00104 
00105         pivotIndices[j] = iMax;
00106 
00107         if (j != iMax)
00108         {
00109             scalar* __restrict__ matrixiMax = matrix[iMax];
00110 
00111             for (register label k=0; k<n; k++)
00112             {
00113                 Swap(matrixj[k], matrixiMax[k]);
00114             }
00115 
00116             vv[iMax] = vv[j];
00117         }
00118 
00119         if (matrixj[j] == 0.0)
00120         {
00121             matrixj[j] = SMALL;
00122         }
00123 
00124         if (j != n-1)
00125         {
00126             scalar rDiag = 1.0/matrixj[j];
00127 
00128             for (register label i=j+1; i<n; i++)
00129             {
00130                 matrix[i][j] *= rDiag;
00131             }
00132         }
00133     }
00134 }
00135 
00136 
00137 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
00138 
00139 void Foam::multiply
00140 (
00141     scalarRectangularMatrix& ans,         // value changed in return
00142     const scalarRectangularMatrix& A,
00143     const scalarRectangularMatrix& B
00144 )
00145 {
00146     if (A.m() != B.n())
00147     {
00148         FatalErrorIn
00149         (
00150             "multiply("
00151             "scalarRectangularMatrix& answer "
00152             "const scalarRectangularMatrix& A, "
00153             "const scalarRectangularMatrix& B)"
00154         )   << "A and B must have identical inner dimensions but A.m = "
00155             << A.m() << " and B.n = " << B.n()
00156             << abort(FatalError);
00157     }
00158 
00159     ans = scalarRectangularMatrix(A.n(), B.m(), scalar(0));
00160 
00161     for(register label i = 0; i < A.n(); i++)
00162     {
00163         for(register label j = 0; j < B.m(); j++)
00164         {
00165             for(register label l = 0; l < B.n(); l++)
00166             {
00167                 ans[i][j] += A[i][l]*B[l][j];
00168             }
00169         }
00170     }
00171 }
00172 
00173 
00174 void Foam::multiply
00175 (
00176     scalarRectangularMatrix& ans,         // value changed in return
00177     const scalarRectangularMatrix& A,
00178     const scalarRectangularMatrix& B,
00179     const scalarRectangularMatrix& C
00180 )
00181 {
00182     if (A.m() != B.n())
00183     {
00184         FatalErrorIn
00185         (
00186             "multiply("
00187             "const scalarRectangularMatrix& A, "
00188             "const scalarRectangularMatrix& B, "
00189             "const scalarRectangularMatrix& C, "
00190             "scalarRectangularMatrix& answer)"
00191         )   << "A and B must have identical inner dimensions but A.m = "
00192             << A.m() << " and B.n = " << B.n()
00193             << abort(FatalError);
00194     }
00195 
00196     if (B.m() != C.n())
00197     {
00198         FatalErrorIn
00199         (
00200             "multiply("
00201             "const scalarRectangularMatrix& A, "
00202             "const scalarRectangularMatrix& B, "
00203             "const scalarRectangularMatrix& C, "
00204             "scalarRectangularMatrix& answer)"
00205         )   << "B and C must have identical inner dimensions but B.m = "
00206             << B.m() << " and C.n = " << C.n()
00207             << abort(FatalError);
00208     }
00209 
00210     ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0));
00211 
00212     for(register label i = 0; i < A.n(); i++)
00213     {
00214         for(register label g = 0; g < C.m(); g++)
00215         {
00216             for(register label l = 0; l < C.n(); l++)
00217             {
00218                 scalar ab = 0;
00219                 for(register label j = 0; j < A.m(); j++)
00220                 {
00221                     ab += A[i][j]*B[j][l];
00222                 }
00223                 ans[i][g] += C[l][g] * ab;
00224             }
00225         }
00226     }
00227 }
00228 
00229 
00230 void Foam::multiply
00231 (
00232     scalarRectangularMatrix& ans,         // value changed in return
00233     const scalarRectangularMatrix& A,
00234     const DiagonalMatrix<scalar>& B,
00235     const scalarRectangularMatrix& C
00236 )
00237 {
00238     if (A.m() != B.size())
00239     {
00240         FatalErrorIn
00241         (
00242             "multiply("
00243             "const scalarRectangularMatrix& A, "
00244             "const DiagonalMatrix<scalar>& B, "
00245             "const scalarRectangularMatrix& C, "
00246             "scalarRectangularMatrix& answer)"
00247         )   << "A and B must have identical inner dimensions but A.m = "
00248             << A.m() << " and B.n = " << B.size()
00249             << abort(FatalError);
00250     }
00251 
00252     if (B.size() != C.n())
00253     {
00254         FatalErrorIn
00255         (
00256             "multiply("
00257             "const scalarRectangularMatrix& A, "
00258             "const DiagonalMatrix<scalar>& B, "
00259             "const scalarRectangularMatrix& C, "
00260             "scalarRectangularMatrix& answer)"
00261         )   << "B and C must have identical inner dimensions but B.m = "
00262             << B.size() << " and C.n = " << C.n()
00263             << abort(FatalError);
00264     }
00265 
00266     ans = scalarRectangularMatrix(A.n(), C.m(), scalar(0));
00267 
00268     for(register label i = 0; i < A.n(); i++)
00269     {
00270         for(register label g = 0; g < C.m(); g++)
00271         {
00272             for(register label l = 0; l < C.n(); l++)
00273             {
00274                 ans[i][g] += C[l][g] * A[i][l]*B[l];
00275             }
00276         }
00277     }
00278 }
00279 
00280 
00281 Foam::RectangularMatrix<Foam::scalar> Foam::SVDinv
00282 (
00283     const scalarRectangularMatrix& A,
00284     scalar minCondition
00285 )
00286 {
00287     SVD svd(A, minCondition);
00288     return svd.VSinvUt();
00289 }
00290 
00291 
00292 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines