src/HEXABLOCK/HexCramer.hxx

   1
   2 // class : Resolution d'un systeme de Cramer
   3
   4 // Copyright (C) 2009-2023  CEA, EDF
   5 //
   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.
  10 //
  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.
  15 //
  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
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 #ifndef __CRAMER_H_
  23 #define __CRAMER_H_
  24
  25 #include "hexa_base.hxx"
  26 // #include <cmath>
  27
  28 BEGIN_NAMESPACE_HEXA
  29
  30 class HexaExport Cramer
  31 {
  32 public:
  33     Cramer (int size);
  34    ~Cramer ();
  35
  36     void defMatrix (double matuser[]);
  37     int  resoudre  (double matr[], double second[],  double sol[]);
  38     void multiply  (double vect[], double result[]);
  39
  40 private:
  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 ; }
  45 private:
  46     const int mat_size;
  47     double** matrix;
  48     double*  col_saved;
  49     Cramer*  co_facteur;
  50 };
  51 // ========================================================= Constructeur
  52 inline Cramer::Cramer (int size)
  53              : mat_size(size)
  54 {
  55    matrix     = new double* [mat_size];
  56    col_saved  = new double  [mat_size];
  57    co_facteur = NULL;
  58    if (mat_size > 2)
  59       co_facteur =  new Cramer (mat_size-1);
  60
  61    for (int nl=0 ; nl<mat_size ; nl++)
  62        {
  63        matrix[nl] = new double [mat_size];
  64        for (int nc=0 ; nc<mat_size ; nc++)
  65            matrix [nl][nc] = 0;
  66        }
  67 }
  68 // ========================================================= Destructeur
  69 inline Cramer::~Cramer ()
  70 {
  71    for (int nl=0 ; nl<mat_size ; nl++)
  72        delete [] matrix[nl];
  73
  74    delete [] matrix;
  75    delete [] col_saved;
  76    delete    co_facteur;
  77 }
  78 // ========================================================= defMatrix
  79 inline void Cramer::defMatrix (double matuser[])
  80 {
  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];
  84 }
  85 // ========================================================= determinant
  86 inline double Cramer::determinant ()
  87 {
  88    double det = 0;
  89    if (mat_size <= 1)
  90       return matrix[0][0];
  91
  92    else if (mat_size == 2)
  93       {
  94       det =matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0];
  95       return det;
  96       }
  97
  98    for (int ni=0 ; ni<mat_size ; ni++)
  99        {
 100        if (NOT cestNul (matrix[ni][0]))
 101           {
 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 ();
 105           }
 106        }
 107    return det;
 108 }
 109 // ========================================================= defCofacteur
 110 inline void Cramer::defCofacteur (Cramer* orig, int lig, int col)
 111 {
 112    for (int nl=0 ; nl<mat_size ; nl++)
 113        {
 114        int ni = nl < lig ? nl : nl+1;
 115        for (int nc=0 ; nc<mat_size ; nc++)
 116            {
 117            int nj = nc < col ? nc : nc+1;
 118            matrix[nl][nc]  = orig->matrix[ni][nj];
 119            }
 120        }
 121 }
 122 // ========================================================= resoudre
 123 inline int Cramer::resoudre (double mat[], double second[],  double sol[])
 124 {
 125    defMatrix (mat);
 126    int ier = resoudre (second, sol);
 127    return ier;
 128 }
 129 // ========================================================= multiply
 130 inline void Cramer::multiply (double facteur[],  double produit[])
 131 {
 132    for (int ni=0 ; ni<mat_size ; ni++)
 133        {
 134        produit [ni] = 0;
 135        for (int nj=0 ; nj<mat_size ; nj++)
 136            produit [ni] += matrix  [ni][nj] * facteur [nj];
 137        }
 138 }
 139 // ========================================================= resoudre
 140 inline int Cramer::resoudre (double msecond[],  double solution[])
 141 {
 142    double det = determinant();
 143    if (cestNul (det))
 144       return HERR;
 145
 146    for (int nj=0 ; nj<mat_size ; nj++)
 147        {
 148        for (int ni=0 ; ni<mat_size ; ni++)
 149            {
 150            col_saved [ni]  = matrix  [ni][nj];
 151            matrix [ni][nj] = msecond [ni];
 152            }
 153        solution [nj] = determinant() / det;
 154
 155        for (int ni=0 ; ni<mat_size ; ni++)
 156            matrix [ni][nj] = col_saved [ni];
 157        }
 158    return HOK;
 159 }
 160 END_NAMESPACE_HEXA
 161 #endif