// FICHIER : Mat_abstraite.h
// CLASSE : Mat_abstraite
// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
// AUTHOR : Gérard Rio
// E-MAIL : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.
/************************************************************************
* DATE: 23/01/97 *
* $ *
* AUTEUR: G RIO (mailto:gerardrio56@free.fr) *
* $ *
* PROJET: Herezh++ *
* $ *
************************************************************************
* La classe Mat_abstraite est une classe abstraite qui permet de cacher les
* implementations specifiques a un type particulier de matrices.
* Cette classe fournit un ensemble de methodes qui seront implementees dans
* les classes derivees caracterisant un type particulier de matrice.
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' * *
* VERIFICATION: *
* *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* ! ! ! ! *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' *
* MODIFICATIONS: *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* $ *
************************************************************************/
// La classe Mat_abstraite est une classe abstraite qui permet de cacher les
// implementations specifiques a un type particulier de matrices.
// Cette classe fournit un ensemble de methodes qui seront implementees dans
// les classes derivees caracterisant un type particulier de matrice.
#ifndef MAT_ABSTRAITE_H
#define MAT_ABSTRAITE_H
#ifdef UTILISATION_MPI
#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>
#endif
#include "mvvtp_GR.h" // classe template MV++
#include "Vecteur.h"
#include "Tableau_T.h"
#include "Enum_matrice.h"
#include "Enum_type_resolution_matri.h"
#include "ExceptionsMatrices.h"
#include "Coordonnee3.h"
class VeurPropre;
/** @defgroup Les_classes_Matrices
*
* BUT: les différentes classes de matrices.
*
*
* \author Gérard Rio
* \version 1.0
* \date 23/01/97
* \brief les différentes classes de matrices.
*
*/
/// @addtogroup Les_classes_Matrices
/// @{
///
class Mat_abstraite
{
protected :
// grandeur par défaut déclarées avant car on s'en sert dans les passages de paramètres
static const double tol_defaut; // tolérance sur le résidu dans les méthodes itératives
static const int maxit_defaut ; // maximum d'itération dans les méthodes itératives
static const int restart_defaut ;// maximum de restart itération (nb de vecteur sauvegardé)
public :
// Constructeur par défaut
Mat_abstraite (Enum_matrice type_mat = RIEN_MATRICE,
Enum_type_resolution_matri type_resol=RIEN_TYPE_RESOLUTION_MATRI,
Enum_preconditionnement type_precondi= RIEN_PRECONDITIONNEMENT ):
type_matrice(type_mat),type_resolution(type_resol),
type_preconditionnement(type_precondi)
{}
// de copie
Mat_abstraite (const Mat_abstraite& a) :
type_matrice(a.type_matrice),type_resolution(a.type_resolution),
type_preconditionnement(a.type_preconditionnement)
{};
// DESTRUCTEUR VIRTUEL :
virtual ~Mat_abstraite ()
{};
// METHODES VIRTUELLES PURES :
// surcharge de l'opérateur d'affectation
// IMPORTANT : le fonctionnement de l'implantation dépend de la classe
// dérivée : en particulier il peut y avoir redimentionnement automatique
// de la matrice en fonction de l'attribut, cas par exemple des matrices
// pleines, ou message d'erreur car les tailles ne sont pas identiques
// exemple de la classe mat_band
virtual Mat_abstraite & operator = ( const Mat_abstraite &) = 0;
// transfert des informations de *this dans la matrice passée en paramètre
// la matrice paramètre est au préalable, mise à 0.
virtual void Transfert_vers_mat( Mat_abstraite & A ) = 0;
//-----------------------------------------------------------------------------------
// --- plusieurs fonctions virtuelles qui agissent en général sur la matrice ---------
//-----------------------------------------------------------------------------------
// ici l'idée est d'éviter de construire une nouvelle matrice pour des questions
// d'encombrement, c'est surtout des méthodes utiles pour le stockage de matrice
// de raideur. On met le nombre mini de fonction, pour pouvoir faire des expressions.
// Cependant aucune de ces fonctions n'est en désaccord avec les fonctions membres
// qui crée une matrice en résultat (cf. Mat_pleine)
// pour les 2 premières méthodes : += et -= :
// les matrices arguments sont en général du mêmes type que celui de la matrice this
// mais cela fonctionne également avec comme argument une matrice diagonale (pour tous les types)
// Surcharge de l'operateur += : addition d'une matrice a la matrice courante
virtual void operator+= (const Mat_abstraite& mat_pl) = 0;
// Surcharge de l'operateur -= : soustraction d'une matrice a la matrice courante
virtual void operator-= (const Mat_abstraite& mat_pl) = 0;
// Surcharge de l'operateur *= : multiplication de la matrice courante par un scalaire
virtual void operator*= (const double r) = 0;
//------------------------------------------------------------------------------------
// --- fin de plusieurs fonctions virtuelles qui agissent en général sur la matrice --
//------------------------------------------------------------------------------------
// Surcharge de l'operateur == : test d'egalite entre deux matrices
virtual int operator== (const Mat_abstraite& mat_pl) const =0;
inline int operator!= (const Mat_abstraite& mat_pl) const
// Surcharge de l'operateur != : test de non egalite entre deux matrices
// Renvoie 1 si les deux matrices ne sont pas egales
// Renvoie 0 sinon
{ if ( (*this)==mat_pl )
return 0;
else
return 1;
};
// fonction permettant de creer une nouvelle instance d'element
// dérivé. la nouvelle instance = *this, il y a donc utilisation du constructeur de copie
virtual Mat_abstraite * NouvelElement() const = 0;
// Affichage des valeurs de la matrice
virtual void Affiche () const =0;
// Initialisation des valeurs des composantes
virtual void Initialise (double val_init) =0;
// Liberation de la place memoire
virtual void Libere () =0;
// Retour du nombre de lignes de la matrice vue comme une matrice_rectangulaire
virtual int Nb_ligne() const =0;
// Retour du nombre de colonnes de la matrice vue comme une matrice_rectangulaire
virtual int Nb_colonne() const =0;
// Symetrie de la matrice
virtual int Symetrie () const =0;
// Acces aux valeurs de la matrices en écriture
// i : indice de ligne, j : indice de colonne
virtual double& operator () (int i, int j) =0;
// ramène true si la place (i,j) existe, false sinon
// ici on ne test pas le fait que i et j puissent être négatif ou pas
virtual bool Existe(int i, int j) const = 0;
// Acces aux valeurs de la matrices en lecture
// cas ou l'on ne veut pas modifier les valeurs
// i : indice de ligne, j : indice de colonne
// !!! important, si l'élément n'existe pas du au type de stockage
// utilisé, (mais que les indices sont possible) le retour est 0
virtual double operator () ( int i, int j) const =0;
// Retourne la ieme ligne de la matrice
// lorsque implemente ( a verifier )
virtual Vecteur& Ligne_set(int i) = 0;
// Retourne la ieme ligne de la matrice
virtual Vecteur Ligne(int i) const = 0;
// Retourne la ieme ligne de la matrice
// sous le format de stokage propre a la matrice
// donc a n'utiliser que comme sauvegarde en parralele
// avec la fonction RemplaceLigne
virtual Vecteur LigneSpe(int i) const = 0;
// remplace la ligne de la matrice par la ligne fournie
virtual void RemplaceLigneSpe(int i,const Vecteur & v) = 0;
//met une valeur identique sur toute la ligne
virtual void MetValLigne(int i,double val) = 0;
// Retourne la jeme colonne de la matrice
virtual Vecteur Colonne(int j) const = 0;
// Retourne la jeme colonne de la matrice
// sous le format de stokage propre a la matrice
// donc a n'utiliser que comme sauvegarde en parralele
// avec la fonction RemplaceColonne
virtual Vecteur ColonneSpe(int j) const = 0;
// remplace la Colonne de la matrice par la colonne fournie
virtual void RemplaceColonneSpe(int j,const Vecteur & v) = 0;
//met une valeur identique sur toute la colonne
virtual void MetValColonne(int j,double val) = 0;
// Resolution du systeme Ax=b
//1) avec en sortie un new vecteur
virtual Vecteur Resol_syst (const Vecteur& b,const double &tol = tol_defaut
,const int maxit = maxit_defaut,const int restart = restart_defaut) =0;
//2) avec en sortie le vecteur d'entree
virtual Vecteur& Resol_systID (Vecteur& b,const double &tol = tol_defaut
,const int maxit = maxit_defaut,const int restart = restart_defaut) =0;
//3) avec en entrée un tableau de vecteur second membre et
// en sortie un nouveau tableau de vecteurs
virtual Tableau <Vecteur> Resol_syst
(const Tableau <Vecteur>& b,const double &tol = tol_defaut
,const int maxit = maxit_defaut,const int restart = restart_defaut) =0;
//4) avec en entrée un tableau de vecteur second membre et
// en sortie le tableau de vecteurs d'entree
virtual Tableau <Vecteur>& Resol_systID
(Tableau <Vecteur>& b,const double &tol = tol_defaut
,const int maxit = maxit_defaut,const int restart = restart_defaut) =0;
//5) avec en sortie le dernier vecteur d'entree, le premier étant le second membre
// et restant inchangé, en sortie c'est donc soit le retour ou soit vortie, les
// deux étant identiques
virtual Vecteur& Resol_systID_2 (const Vecteur& b,Vecteur& vortie
, const double &tol = tol_defaut
,const int maxit = maxit_defaut
,const int restart = restart_defaut) =0;
// ===== RÉSOLUTION EN DEUX TEMPS ================ :
// 1) préparation de la matrice donc modification de la matrice éventuellement
// par exemple pour les matrices bandes avec cholesky : triangulation
virtual void Preparation_resol() =0 ;
// 2) *** résolution sans modification de la matrice DOIT ÊTRE PRÉCÉDÉ DE L'APPEL DE
// Preparation_resol
// a) avec en sortie un new vecteur
virtual Vecteur Simple_Resol_syst (const Vecteur& b,const double &tol = tol_defaut
,const int maxit = maxit_defaut,const int restart = restart_defaut) const =0 ;
// b) avec en sortie le vecteur d'entree
virtual Vecteur& Simple_Resol_systID (Vecteur& b,const double &tol = tol_defaut
,const int maxit = maxit_defaut,const int restart = restart_defaut) const =0 ;
// c) avec en entrée un tableau de vecteur second membre et
// en sortie un nouveau tableau de vecteurs
virtual Tableau <Vecteur> Simple_Resol_syst
(const Tableau <Vecteur>& b,const double &tol = tol_defaut
,const int maxit = maxit_defaut,const int restart = restart_defaut) const =0 ;
// d) avec en entrée un tableau de vecteur second membre et
// en sortie le tableau de vecteurs d'entree
virtual Tableau <Vecteur>& Simple_Resol_systID
(Tableau <Vecteur>& b,const double &tol = tol_defaut
,const int maxit = maxit_defaut,const int restart = restart_defaut) const =0 ;
// e) avec en sortie le dernier vecteur d'entree, le premier étant le second membre
// et restant inchangé, en sortie c'est donc soit le retour ou soit vortie, les
// deux étant identiques
virtual Vecteur& Simple_Resol_systID_2 (const Vecteur& b,Vecteur& vortie
, const double &tol = tol_defaut
,const int maxit = maxit_defaut
,const int restart = restart_defaut) const =0 ;
// ===== FIN RÉSOLUTION EN DEUX TEMPS ================ :
// Multiplication d'un vecteur par une matrice ( (vec)t * A )
virtual Vecteur Prod_vec_mat ( const Vecteur& vec) const =0;
// idem mais on utilise la place du second vecteur pour le résultat
virtual Vecteur& Prod_vec_mat ( const Vecteur& vec, Vecteur & resul) const =0;
// Multiplication d'une matrice par un vecteur ( A * vec )
virtual Vecteur Prod_mat_vec ( const Vecteur& vec) const =0;
// idem mais on utilise la place du second vecteur pour le résultat
virtual Vecteur& Prod_mat_vec ( const Vecteur& vec, Vecteur & resul) const =0;
// Multiplication d'une ligne iligne de la matrice avec un vecteur de
// dimension = le nombre de colonne de la matrice
virtual double Prod_Ligne_vec ( int iligne,const Vecteur& vec) const =0;
// Multiplication d'un vecteur avec une colonne icol de la matrice
// dimension = le nombre de ligne de la matrice
virtual double Prod_vec_col( int icol,const Vecteur& vec) const =0;
// calcul du produit : (vec_1)^T * A * (vect_2)
virtual double vectT_mat_vec(const Vecteur& vec1, const Vecteur& vec2) const =0;
// retourne la place que prend la matrice en entier
virtual int Place() const = 0;
// ------------ méthode non virtuelle pur -> générales ----------------
// Affiche une partie de la matrice (util pour le debug)
// min_ et max_ sont les bornes de la sous_matrice
// pas_ indique le pas en i et j pour les indices
virtual void Affiche1(int min_i,int max_i,int pas_i,int min_j,int max_j,int pas_j) const ;
// Affiche une partie de la matrice idem si dessus
// mais avec un nombre de digit (>7) = nd
// si < 7 ne fait rien
virtual void Affiche2
(int min_i,int max_i,int pas_i,int min_j,int max_j,int pas_j,int nd) const ;
// changement du choix de la méthode de résolution si c'est possible
// peut être surchargé par les classes dérivées en fonction des spécificités
virtual void Change_Choix_resolution(Enum_type_resolution_matri type_resol,
Enum_preconditionnement type_precondi)
{ type_resolution = type_resol; type_preconditionnement = type_precondi;};
// retourne le type de la matrice
inline const Enum_matrice Type_matrice() const {return type_matrice;};
// affichage à l'écran après une demande interactive de la matrice
// entete: une chaine explicative, a afficher en entête
void Affichage_ecran(string entete) const;
// calcul, récupération et affichage éventuelle
// des mini, maxi, et en valeur absolue la moyenne de la diagonale de la matrice
// en retour: le min, le max et la moyenne en valeur absolue
Coordonnee3 MinMaxMoy(bool affiche) const;
// limite la valeur mini de la diagonale: si un terme de la diagonale
// est inférieure à la limite passée en argument (seuil_bas), d'une manière arbitraire
// le terme à mis à une valeur égale au paramètre passé en paramètre val_a_imposer
// s'il y a eu changement retour de true
bool Limitation_min_diag(double seuil_bas, double val_a_imposer) ;
//====== quelques méthodes spécifiques à la classe MV_Vector<double>==============
// définie la surcharge de multiplication d'une matrice par un MV_Vector
// ici de type constant, ramène un nouveau MV_Vector<double>
MV_Vector <double> operator * (const MV_Vector <double> & vec) const ;
// multiplication matrice transposée fois vecteur ce qui est équivalent à
// la transposée du résultat de : multiplication vecteur transposé fois matrice
virtual MV_Vector <double> trans_mult(const MV_Vector <double> &x) const ;
//====== fin des quelques méthodes spécifiques à la classe MV_Vector <double>=======
//calcul des valeurs propres par la méthode des puissances itérées
// inverses.
// --> Utilisable que pour des matrices carrées
// Intéressant lorsque l'on veut les plus petites
// valeurs propres, en nombre restreind car seule les premières valeurs
// sont calculées avec une précision convenable. On obtiend également
// les vecteurs propres correspondant.
// résolution de : ((*this) - lambda KG) X = 0
// en entrée : KG, VP dont la dimension est défini et fourni le nombre
// de valeurs propres que l'on veut calculer
// On considère que les conditions limites sont déjà appliquées à (*this) et KG.
// c'est à dire des 1 sur la diagonale pour (*this) et des 0 sur la diagonale pour KG
// en sortie : VP contiend les valeurs propres et vecteurs propres
Tableau <VeurPropre>* V_Propres(Mat_abstraite& KG,Tableau <VeurPropre> & VP );
protected :
// VARIABLES PROTEGEES :
Enum_matrice type_matrice; // le type de la matrice, défini dans les classes dérivée
Enum_type_resolution_matri type_resolution; // le type de résolution envisagée
Enum_preconditionnement type_preconditionnement; // le type éventuelle de préconditionnement
// pour la méthode "résidu minimal généralisé"
#ifdef UTILISATION_MPI
friend class boost::serialization::access;
// When the class Archive corresponds to an output archive, the
// & operator is defined similar to <<. Likewise, when the class Archive
// is a type of input archive the & operator is defined similar to >>.
template<class Archive>
void serialize(Archive & ar, const unsigned int version)
{
ar & type_matrice;
ar & type_resolution;
ar & type_preconditionnement;
}
#endif
// Méthodes protégées :
// Resolution du systeme Ax=b , méthodes générales pouvant être appelées par les classes dérivée
// en entrée : b : comme second membre
// tol : tolérance sur le résidu
// maxit : le maximum d'itération
// restart : le Maximum de restart iterations (vecteur sauvegardé) dans le cas
// de la méthode résidu minimal généralisé (GENE_MINI_RESIDUAL)
// en sortie : sol de même dimension que sol
void Resolution_syst
(const Vecteur &b, Vecteur &sol,const double &tol, const int maxit = 300
, const int restart = 32) const ;
// idem avec en entrée un tableau de vecteur second membre et
// en sortie un tableau de vecteurs solution
void Resolution_syst
(const Tableau <Vecteur>& b, Tableau <Vecteur> &sol,const double &tol
, const int maxit = 300
, const int restart = 32) const ;
};
/// @} // end of group
#include "VeurPropre.h"
#endif