2 // class : Resolution d'un systeme de Cramer
4 // Copyright (C) 2009-2023 CEA/DEN, EDF R&D
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
25 #include "hexa_base.hxx"
30 class HexaExport Cramer
36 void defMatrix (double matuser[]);
37 int resoudre (double matr[], double second[], double sol[]);
38 void multiply (double vect[], double result[]);
41 int resoudre (double msecond[], double solution[]);
42 double determinant ();
43 void defCofacteur (Cramer* orig, int lig, int col);
44 bool cestNul (double x) { return x > -1e-10 && x < 1e-10 ; }
51 // ========================================================= Constructeur
52 inline Cramer::Cramer (int size)
55 matrix = new double* [mat_size];
56 col_saved = new double [mat_size];
59 co_facteur = new Cramer (mat_size-1);
61 for (int nl=0 ; nl<mat_size ; nl++)
63 matrix[nl] = new double [mat_size];
64 for (int nc=0 ; nc<mat_size ; nc++)
68 // ========================================================= Destructeur
69 inline Cramer::~Cramer ()
71 for (int nl=0 ; nl<mat_size ; nl++)
78 // ========================================================= defMatrix
79 inline void Cramer::defMatrix (double matuser[])
81 for (int nl=0 ; nl<mat_size ; nl++)
82 for (int nc=0 ; nc<mat_size ; nc++)
83 matrix [nl][nc] = matuser [nl*mat_size + nc];
85 // ========================================================= determinant
86 inline double Cramer::determinant ()
92 else if (mat_size == 2)
94 det =matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0];
98 for (int ni=0 ; ni<mat_size ; ni++)
100 if (NOT cestNul (matrix[ni][0]))
102 co_facteur->defCofacteur (this, ni, 0);
103 double signe = (ni MODULO 2) ? -1.0 : 1.0;
104 det += signe*matrix[ni][0] * co_facteur->determinant ();
109 // ========================================================= defCofacteur
110 inline void Cramer::defCofacteur (Cramer* orig, int lig, int col)
112 for (int nl=0 ; nl<mat_size ; nl++)
114 int ni = nl < lig ? nl : nl+1;
115 for (int nc=0 ; nc<mat_size ; nc++)
117 int nj = nc < col ? nc : nc+1;
118 matrix[nl][nc] = orig->matrix[ni][nj];
122 // ========================================================= resoudre
123 inline int Cramer::resoudre (double mat[], double second[], double sol[])
126 int ier = resoudre (second, sol);
129 // ========================================================= multiply
130 inline void Cramer::multiply (double facteur[], double produit[])
132 for (int ni=0 ; ni<mat_size ; ni++)
135 for (int nj=0 ; nj<mat_size ; nj++)
136 produit [ni] += matrix [ni][nj] * facteur [nj];
139 // ========================================================= resoudre
140 inline int Cramer::resoudre (double msecond[], double solution[])
142 double det = determinant();
146 for (int nj=0 ; nj<mat_size ; nj++)
148 for (int ni=0 ; ni<mat_size ; ni++)
150 col_saved [ni] = matrix [ni][nj];
151 matrix [ni][nj] = msecond [ni];
153 solution [nj] = determinant() / det;
155 for (int ni=0 ; ni<mat_size ; ni++)
156 matrix [ni][nj] = col_saved [ni];