Herezh_dev/Resolin/Matrices/Mat_pleine.h

// FICHIER : Mat_pleine.h
// CLASSE : Mat_pleine

// 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_pleine derive de la classe Mat_abstraite et permet de declarer
 * une matrice pleine dont les composantes sont de type double. Toute instance
 * de cette classe n'a pas de restrictions particulieres et est capable de representer
 * n'importe quel type de matrices.
 *
 * Une telle matrice est stockee a l'aide d'un tableau de vecteurs à l'aide des classes
 * Tableau<T> et Vecteur de facon a beneficier des methodes introduites au niveau de
 * ces dernieres.
 *                                                                $     *
 *     ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''     *                                                                      *
 *     VERIFICATION:                                                    *
 *                                                                      *
 *     !  date  !   auteur   !       but                          !     *
 *     ------------------------------------------------------------     *
 *     !        !            !                                    !     *
 *                                                                $     *
 *     ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''     *
 *     MODIFICATIONS:                                                   *
 *     !  date  !   auteur   !       but                          !     *
 *     ------------------------------------------------------------     *
 *                                                                $     *
 ************************************************************************/

// La classe Mat_pleine derive de la classe Mat_abstraite et permet de declarer 
// une matrice pleine dont les composantes sont de type double. Toute instance
// de cette classe n'a pas de restrictions particulieres et est capable de representer
// n'importe quel type de matrices.
//
// Une telle matrice est stockee a l'aide d'un tableau de vecteurs à l'aide des classes
// Tableau<T> et Vecteur de facon a beneficier des methodes introduites au niveau de 
// ces dernieres.


#ifndef MAT_PLEINE_H
#define MAT_PLEINE_H


#include "Mat_abstraite.h"
#include "Tableau_T.h"
#include "Vecteur.h"
#include "Coordonnee.h"
#include "Coordonnee3.h"
#ifdef UTILISATION_MPI
  #include "mpi.h"
  #include <boost/mpi/environment.hpp>
  #include <boost/mpi/communicator.hpp>
  #include <boost/serialization/string.hpp>
  #include <boost/mpi.hpp>
  namespace mpi = boost::mpi;
#endif
# include "ParaGlob.h"

class Mat_abstraite;

/// @addtogroup Les_classes_Matrices
///  @{
///


class Mat_pleine : public Mat_abstraite
{      // surcharge de l'operator de lecture typée
       friend istream & operator >> (istream &, Mat_pleine &);
       // surcharge de l'operator d'ecriture typée
       friend ostream & operator << (ostream &, const Mat_pleine &);
       // ResRaid_MPI permet de faire un stockage globale de tranche de matrice
       friend class ResRaid_MPI;
	public : 
		// CONSTRUCTEURS :
		
		// Constructeur par defaut
		Mat_pleine ();
		
		// Constructeur se servant d'un nombre de lignes et de colonnes, et eventuellement
		// d'une valeur d'initialisation
		Mat_pleine (int nb_ligne,int nb_colonne,double val_init=0.0);
		
		// Constructeur de copie
		Mat_pleine (const Mat_pleine& mat_pl);
		
		
		// DESTRUCTEUR :
		
		virtual ~Mat_pleine ();
		
		
		// METHODES :
		
  // fonction permettant de creer une nouvelle instance d'element
  Mat_abstraite * NouvelElement() const ;
		
  // surcharge de l'opérateur d'affectation : cas de matrices abstraites
  // il y a ajustement de la taille de la matrice en fonction de mat_pl
  Mat_abstraite & operator = ( const Mat_abstraite & mat_pl);
 
  // transfert des informations de *this dans la matrice passée en paramètre
  // la matrice paramètre est au préalable, mise à 0.
  void Transfert_vers_mat( Mat_abstraite & A );
 
  //-----------------------------------------------------------------------------------
  // --- 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)
		
		// Surcharge de l'operateur += : addition d'une matrice a la matrice courante
		void operator+= (const Mat_abstraite& mat_pl);

		// Surcharge de l'operateur -= : soustraction d'une matrice a la matrice courante
		void operator-= (const Mat_abstraite& mat_pl);
		
        //------------------------------------------------------------------------------------
        // --- fin de plusieurs fonctions virtuelles qui agissent en général sur la matrice --
        //------------------------------------------------------------------------------------
        
		// Surcharge de l'operateur == : test d'egalite entre deux matrices
		int operator== (const Mat_abstraite& mat_pl) const ;
		
		// Surcharge de l'operateur = : affectation d'une matrice a une autre
        // il y a ajustement de la taille de la matrice en fonction de mat_pl
		Mat_pleine& operator= (const Mat_pleine& mat_pl);
		
		// Affichage des donnees de la matrice
		void Affiche () const ;
		
		// Initialisation des composantes de la matrice (par defaut a 0, sinon a val_init)		
		void Initialise (double val_init=0.0)
      {  val_de_base.Inita(val_init);};
		
		// Initialisation d'une matrice a un nombre de lignes et de colonnes ainsi que
		// ces composantes a 0 par defaut ou a val_init sinon
		void Initialise (int nb_ligne,int nb_colonne,double val_init=0.0);
		
		inline void Libere ()
		// A la suite de l'appel de cette methode, la matrice est identique a ce qu'elle
		// serait a la suite d'un appel du constructeur par defaut
		{	val_de_base.Libere();	};
		
		// Retour du nombre de lignes de la matrice vue comme une matrice_rectangulaire
		inline int Nb_ligne () const
		{	return val.Taille();	};
				
		// Retour du nombre de colonnes de la matrice vue comme une matrice_rectangulaire
		// N.B. : La taille des vecteurs du tableau val est supposee etre la meme pour tous
		// les vecteurs du tableau (aucun test n'est realise pour le verifier)
		inline int Nb_colonne () const
		{	if ( Nb_ligne()!=0 )
				return val(1)->Taille();
			else 
				return 0;
		};
		
		// Retourne la ieme ligne de la matrice
		inline Vecteur Ligne (int i) const 
		{	return *(val(i));	};
		
		// 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
		inline Vecteur LigneSpe (int i) const 
		{	return *(val(i));	};
		// remplace la ligne de la matrice par la ligne fournie
		void RemplaceLigneSpe(int i,const Vecteur & v);
		//met une valeur identique sur toute la ligne
		void MetValLigne(int i,double x);
		
		// Retourne la jieme colonne de la matrice		
		inline Vecteur Colonne (int j) const 
		{	Vecteur result(Nb_ligne());
			for (int i=1;i<=Nb_ligne();i++)
				result(i)=(*this)(i,j);
			return result;
		};
		
		// 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
		inline Vecteur ColonneSpe (int j) const 
		// Retourne la jieme colonne de la matrice
		{	return Colonne(j);
		};
		// remplace la Colonne de la matrice par la colonne fournie
		void RemplaceColonneSpe(int j,const Vecteur & v);
		//met une valeur identique sur toute la colonne
		void MetValColonne(int j,double y);
		
		// Test sur la symetrie de la matrice
		int Symetrie () const ;
		// affiche les termes non symétriques de la matrice s'il y en a
		void AfficheNonSymetries() const;
		
		// Resolution du systeme Ax=b par la methode de Cholesky, de cramer ou appel de la résolution générale 
		// par exemple gradient conjugué
		    //1) avec en sortie un new vecteur
		Vecteur Resol_syst (const Vecteur& b,const double &tol = tol_defaut
		       ,const int maxit = maxit_defaut,const int restart = restart_defaut);
		    //2) avec en sortie le vecteur d'entree 
		Vecteur& Resol_systID (Vecteur& b,const double &tol = tol_defaut
		       ,const int maxit = maxit_defaut,const int restart = restart_defaut);
		    //3) avec en entrée un tableau de vecteur second membre et
		    //    en sortie un nouveau tableau de vecteurs  
		Tableau <Vecteur> Resol_syst 
		        (const Tableau <Vecteur>& b,const double &tol = tol_defaut
		       ,const int maxit = maxit_defaut,const int restart = restart_defaut) ;
		    //4) avec en entrée un tableau de vecteur second membre et
		    //    en sortie le tableau de vecteurs d'entree 
		Tableau <Vecteur>& Resol_systID 
		       (Tableau <Vecteur>& b,const double &tol = tol_defaut
		       ,const int maxit = maxit_defaut,const int restart = restart_defaut) ;
		//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
  Vecteur& Resol_systID_2 (const Vecteur& b,Vecteur& vortie
		                     , const double &tol = tol_defaut,const int maxit = maxit_defaut
		                     ,const int restart = restart_defaut);
		                     
		// ===== 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 
		void Preparation_resol();
		//    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
		Vecteur Simple_Resol_syst (const Vecteur& b,const double &tol = tol_defaut
		       ,const int maxit = maxit_defaut,const int restart = restart_defaut) const ;
		      // b) avec en sortie le vecteur d'entree 
		Vecteur& Simple_Resol_systID (Vecteur& b,const double &tol = tol_defaut
		       ,const int maxit = maxit_defaut,const int restart = restart_defaut) const ;
		      // c) avec en entrée un tableau de vecteur second membre et
		      //    en sortie un nouveau tableau de vecteurs  
		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 ;
		      // d) avec en entrée un tableau de vecteur second membre et
		      //    en sortie le tableau de vecteurs d'entree 
		Tableau <Vecteur>& Simple_Resol_systID 
		                        (Tableau <Vecteur>& b,const double &tol = tol_defaut
		       ,const int maxit = maxit_defaut,const int restart = restart_defaut) const ;
		      // 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
		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 ;
		// ===== FIN RÉSOLUTION EN DEUX TEMPS ================ :
		
		// Produit d'un vecteur par une matrice
		Vecteur Prod_vec_mat (  const Vecteur& vec) const;
		// idem mais on utilise la place du second vecteur pour le résultat
		Vecteur& Prod_vec_mat ( const Vecteur& vec, Vecteur & resul) const ;		
		// Produit d'une matrice par un vecteur
		Vecteur Prod_mat_vec ( const Vecteur& vec) const;
		// idem mais on utilise la place du second vecteur pour le résultat
		Vecteur& Prod_mat_vec ( const Vecteur& vec, Vecteur & resul) const;
		// Determine la transpose d'une matrice
		Mat_pleine Transpose()  const;
		
  // determine l'inverse d'une matrice carre
  // actuellement uniquement implemente pour dim = 1,2,3
  Mat_pleine Inverse()  const;
  // idem mais en retournant l'inverse dans la matrice passée en paramètre
  // qui doit être de même type et de même dimension que this
  Mat_pleine& Inverse(Mat_pleine& res) const;
 
  // détermine le déterminant d'une matrice
  // actuellement uniquement implemente pour dim = 1,2,3
  double Determinant() const;
 
		// Surcharge de l'operateur + : addition de deux matrices
		Mat_pleine operator+ (const Mat_pleine& mat_pl) const;
		
		// Surcharge de l'operateur - : soustraction entre deux matrices
		Mat_pleine operator- (const Mat_pleine& mat_pl)  const;
		
		// Surcharge de l'operateur - : oppose d'une matrice 
		Mat_pleine operator- () const;
		
		// Surcharge de l'operateur += : addition d'une matrice a la matrice courante
		void operator+= (const Mat_pleine& mat_pl);
		
		// Surcharge de l'operateur -= : soustraction d'une matrice a la matrice courante
		void operator-= (const Mat_pleine& mat_pl);
		
		// Surcharge de l'operateur * : multiplication de deux matrices
		Mat_pleine operator* (const Mat_pleine& mat_pl) const;
		
		// Surcharge de l'operateur *= : multiplication de la  matrice courante par un scalaire
		void operator*= (const double r);
		
		inline Vecteur operator* (const Vecteur& vec) const
		// Realise la multiplication d'une matrice par un vecteur
		{	return Prod_mat_vec((Vecteur&) vec);	};
		
		// Surcharge de l'operateur * : multiplication d'une matrice par un scalaire
		Mat_pleine operator* (const double coeff) const;
		
		// Surcharge de l'operateur == : test d'egalite entre deux matrices
		int operator== (const Mat_pleine& mat_pl) const ;
		
		inline int operator!= (const Mat_pleine& 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;
		};
		
		inline Vecteur& Ligne_set(int i)
		// Retourne la ieme ligne de la matrice
		{	return *(val(i));	};
		
		inline Vecteur& operator() (int i) 
		// Retourne également  de la ieme ligne de la matrice
		// (on a changé: un moment on ne ramenait qu'une copie, mais c'est trop dangereux)
		// ( car naturellement on s'attend à avoir la ligne et cela engendre des erreurs) 
		// ( très difficile à trouver, car tant que la copie persiste c'est ok, et quand le destructeur)
		// ( est appelé -> plus rien !! ) 
		{	return *(val(i));	};

		inline Vecteur operator() (int i) const
		// Retourne une copie  de la ieme ligne de la matrice
		// utile quand on travaille avec des const
		{	return *(val(i));	};
		
		inline double& operator () (int i,int j)
		// Retourne la ieme jieme composante de la matrice
		{	return (*(val(i)))(j);	};
		
		// 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
		inline bool Existe(int , int ) const {return true;};

		// Acces aux valeurs de la matrices
		// cas ou l'on ne veut pas modifier les valeurs
		// i : indice de ligne, j : indice de colonne
		inline double operator () ( int i, int j) const 
		  {	return (*(val(i)))(j);	};
		
  // Multiplication d'une ligne iligne de la matrice avec un vecteur de
  // dimension = le nombre de colonne de la matrice
  double Prod_Ligne_vec ( int iligne, const Vecteur& vec) const;

  // Multiplication d'un vecteur avec une colonne icol de la matrice 
  // dimension = le nombre de ligne de la matrice
  double Prod_vec_col( int icol, const Vecteur& vec) const;
		
		// calcul du produit : (vec_1)^T * A * (vect_2)
		double vectT_mat_vec(const Vecteur& vec1, const Vecteur& vec2) const ;

  // mise a zero de tous les composantes de la matrice
  void Zero()
     {val_de_base.Inita(0.);};
		  
  // retourne la place que prend la matrice en entier
  int Place() const { return (2*(Nb_ligne() * Nb_colonne()));};
 
  // --------- méthodes particulières (en particulier non virtuelles) ----------
 
  // ramène un tableau de coordonnées (avec 3 variances possibles) correspondant à la matrice
  // tab(i) = la colonne i de la matrice
  Tableau <Coordonnee > Coordonnee_Base_associee() const;
  Tableau <CoordonneeH > CoordonneeH_Base_associee() const;
  Tableau <CoordonneeB > CoordonneeB_Base_associee() const;
  // idem mais une seule colonne i
  Coordonnee  Coordonnee_Base_associee(int i) const;
  CoordonneeH  CoordonneeH_Base_associee(int i) const;
  CoordonneeB  CoordonneeB_Base_associee(int i) const;
 
 
  // ramène le maxi en valeur absolue des valeurs de la matrice, et les indices associées
  double MaxiValAbs(int & i, int & j) const;
  // ramène le maxi en valeur absolue de la somme des valeurs absolues de chaque ligne,
  // et l'indice de ligne associé
  // permet d'avoir une approximation de la valeur propre maximale de la matrice via le
  // théorème de Gerschgorin : |lambda| < Max_i (|lambda_i|)
  // avec:  |lambda_i| < Max_j somme_j |k_ij|)
  double Maxi_ligne_ValAbs(int & i) const;

   
#ifdef UTILISATION_MPI // spécifique au calcul parallèle
  // envoi et récup brut du tableau de composantes : nécessite pour la réception d'avoir
  // un tableau de même dimension
  void Envoi_MPI(int dest, int tag) const
      {val_de_base.Envoi_MPI(dest,tag);};
  void Recup_MPI(int source, int tag) const
      {val_de_base.Recup_MPI(source,tag);};
  mpi::request Ienvoi_MPI(int dest, int tag) const
      {return val_de_base.Ienvoi_MPI(dest,tag);};
  mpi::request Irecup_MPI(int source, int tag) const
      {return val_de_base.Irecup_MPI(source,tag);};
#endif


	protected : 
	
	 Vecteur val_de_base; // le vecteur qui contient tous les vecteurs de val
		Tableau<Vecteur * > val;	// Valeurs des composantes: il s'agit de pointeurs
             // sur des vecteurs liés: ils sont à la même place que val_de_base
		// val(i) = une ligne
  
  // fonction protégé pour l'accès directe au stockage global, utilisé
  // par ResRaid_MPI par exemple: pour un stockage par tranche de matrice
  Vecteur & Val_de_base() {return val_de_base;}
  
  /// --------- cas très particulier d'un vecteur val_de_base lié
  /// --------- utilisé uniquement par des classes friend
  // c-a-d qui ne gère pas sa mémoire -> memoire = false
  /// il s'agit d'un tuillage, val_de_base est positionné a "vec"
  /// aucun test n'est fait concernant la taille disponible à partir de "vec"
  /// routine dangereuse !!!!!
  /// avec_recopie == true : les infos de mat_pl sont recopiées, sinon non
  Mat_pleine(bool avec_recopie,double* vec, const Mat_pleine& mat_pl);

  // méthode interne pour créer val en fonction des dimensions
  // ce qui crée une liaison entre les deux stockages
  void Creation_tab_val(int nb_lig,int nb_col);
  // idem pour modifier éventuellement les tailles de val uniquement
  void Change_tailles_tab_val(int nb_lig,int nb_col);
                                  
#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)
  {   int nb_ligne = val.Taille();
      int nb_colonne = Nb_colonne();
      ar & nb_ligne;
      ar & nb_colonne;
      val_de_base.Change_taille(nb_ligne*nb_colonne);
      Change_tailles_tab_val(nb_ligne,nb_colonne); // en écriture, ne fait rien,
                                 // en lecture permet d'avoir la bonne taille
      // puis le vecteur global
      ar & val_de_base;
  }
#endif

  // calcul de la matrice triangulée dans le cadre de la méthode de cholesky
  // utilisation particulière : la matrice résultat B est défini dans le programme
  // appelant.
  void Triangulation (int N,Mat_pleine& B);
  // résolution du problème triangularisé
  // second membre b, et résultat res
  void Resolution (int N, const Mat_pleine& B, const Vecteur& b,Vecteur& res) const ; 
  // résolution directe par la méthode de cramer (pour les petits systèmes dim 1,2,3)
  void Cramer( const Vecteur& b,Vecteur& res) const;
  // idem pour un tableau de vecteur
  void Cramer( const Tableau <Vecteur>& b,Tableau <Vecteur>& res) const; 
	 // symétrisation de la matrice (il faut qu'elle soit carrée
	 void Symetrisation();    
		
};
/// @}  // end of group


inline Vecteur operator* (const Vecteur& vec,const Mat_pleine& mat_pl)
// Permet de realiser la multiplication entre un vecteur et une matrice
  {	return mat_pl.Prod_vec_mat((Vecteur&) vec);
  };

inline Mat_pleine operator* (const double coeff,const Mat_pleine& mat_pl)
// Permet de realiser la multiplication entre un scalaire et une matrice
  {	return (mat_pl*coeff);	
  };

#include "Mat_abstraite.h"

#ifndef MISE_AU_POINT
  #include "Mat_pleine.cc"
  #define  Mat_pleine_H_deja_inclus
#endif


#endif