src/HEXABLOCK/HexCramer.hxx

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