// 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++ *
* $ *
************************************************************************
* BUT: definition de matrices bandes symetriques a valeurs reelles.*
* Seule la partie inferieure est stockee. *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' * *
* VERIFICATION: *
* *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* ! ! ! ! *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' *
* MODIFICATIONS: *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* $ *
************************************************************************/
#ifndef MATBAND_H
#define MATBAND_H
//#include "Debug.h"
#include <iostream>
#include "Mat_abstraite.h"
#include "MathUtil.h"
/// @addtogroup Les_classes_Matrices
/// @{
///
class MatBand : public Mat_abstraite
{ // surcharge de l'operator de lecture typée
friend istream & operator >> (istream &, MatBand &);
// surcharge de l'operator d'ecriture typée
friend ostream & operator << (ostream &, const MatBand &);
public :
// CONSTRUCTEURS :
MatBand (); // par defaut
MatBand (Enum_matrice type_mat , int lb, int dim ); // def d'une matrice bande
MatBand (Enum_matrice type_mat, int lb, int dim , double a); // def d'une matrice bande avec
// initialisation des composantes a la valeur a
MatBand (const MatBand& m) ; // de copie, il y a creation d'une deuxieme matrice
// DESTRUCTEUR :
~MatBand ();
// fonction permettant de creer une nouvelle instance d'element
Mat_abstraite * NouvelElement() const;
// METHODES PUBLIQUES :
// surcharge de l'opérateur d'affectation : cas de matrices abstraites
Mat_abstraite & operator = ( const Mat_abstraite &);
// 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 );
// surcharge de l'opérateur d'affectation : cas de matrices bandes
MatBand & operator = ( const MatBand &);
//-----------------------------------------------------------------------------------
// --- 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
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);
// Surcharge de l'operateur *= : multiplication de la matrice courante par un scalaire
void operator*= (const double r);
//------------------------------------------------------------------------------------
// --- 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 ;
// acces aux composantes comme pour une matrice carree
// l'utilisateur doit gerer les depassements de taille
// il y a un message en debuggage
double& operator ()(int i, int 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 i, int j) const {return ((abs(i-j)-lb) < 0);};
// acces en lecture seul, aux composantes comme pour une matrice carree
// l'utilisateur doit gerer les depassements de taille
// il y a un message en debuggage
// si l'indice est possible (c'est à dire >0 et < à la dimension) mais
// hors la bande -> retour 0
double operator () (int i, int j ) const ;
// Retourne la ieme ligne de la matrice
// pas util dans le cas des bande donc
// non implemente
Vecteur& Ligne_set(int i)
{ cout <<" fonction non implementee : "
<< " Vecteur& MatBand::operator() (int i) " << endl;
Sortie(1);
i = i; // pour ne pas avoir de message a la compilation !!
Vecteur* toto; toto = new Vecteur();
return *toto;
};
// Retourne la ieme ligne de la matrice uniquement en lecture
// pas utile dans le cas des bandes donc
// non implementee
Vecteur operator() (int i) const
{ cout <<" fonction non implementee : "
<< " Vecteur& MatBand::operator() (int i) " << endl;
i = i; // pour ne pas avoir de message a la compilation !!
return Vecteur(0);
};
// Retourne la ieme ligne de la matrice
Vecteur Ligne(int i) const ;
// 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
Vecteur LigneSpe(int i) const ;
// 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 val);
// Retourne la jeme colonne de la matrice
Vecteur Colonne(int j) const ;
// 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
Vecteur ColonneSpe(int j) const ;
// remplace la Colonne de la matrice par la ligne fournie
void RemplaceColonneSpe(int j, const Vecteur & v);
//met une valeur identique sur toute la colonne
void MetValColonne(int j,double val);
// Affichage des valeurs de la matrice
// uniquement le valeurs de la bande inferieur
void Affiche() const ;
// 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
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
void Affiche2(int min_i,int max_i,int pas_i,int min_j,int max_j,int pas_j,int nd) const ;
void Change_taille(int lb, int dim ); // changement de la taille de la matrice
inline int Nb_colonne () const{ return lb;} ; // retour de la largeur de bande
inline int Nb_ligne() const { return dim;}; // retour de la dimension de la matrice
void Initialise (double a); // initialisation de la matrice a la valeur "a"
void Libere () ; // Liberation de la place memoire
int Symetrie () const { return 1; } // la matrice est Symetrique
// Resolution du systeme Ax=b
// la verification de taille n'est faite que pour le debug
//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 ================ :
// Multiplication d'un vecteur par une matrice ( (vec)t * A )
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 ;
// Multiplication d'une matrice par un vecteur ( A * vec )
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 ;
// 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 ;
// retourne la place que prend la matrice en entier
int Place() const { return (dim * lb );};
// ++++ méthodes spécifique à la classe
// Resolution du systeme Ax=b, sans triangulation, c'est-à-dire que l'on
// considère que la matrice stocke actuellement la triangulation
Vecteur& Resol_systID_sans_triangul ( Vecteur& b);
protected :
// VARIABLES PROTEGEES :
double * pt; // pointeur de la matrice
int lb; // largeur de la bande
int dim; // dimension de la matrice
// METHODES PROTEGEES :
//reduction de la matrice de raideur avant resolution
void REDUCT (MatBand & TRAID,int LB,int NDDL,Vecteur& CC);
// resolution
void RESOLT(const MatBand & TRAID,Vecteur & TSM,int NDDL,int LB) const;
// enchainement de la resolution avec controle des pivots
void ResoBand (MatBand & TRAID,int LB,int NDDL,Vecteur& CC);
// idem précédent mais avec plusieurs seconds membres
void ResoBand (MatBand & TRAID,int LB,int NDDL,Tableau <Vecteur>& CC);
// acces aux coordonnées bandes : utilisés uniquement pour la résolution:
// de reduct et rosolt
// l'adressage est différent de celui de la surcharge de l'opérateur ()
// donc à ne pas utiliser autre part !!
inline double& c(int ii, int jj )
{
#ifdef MISE_AU_POINT
if (ii>lb || jj>dim )
{cout << "erreur d acces aux composantes bande";
cout << " ii = " << ii <<" jj = " << jj << '\n';
cout << "double& c(int ii, int jj )" << endl;
}
#endif
return pt[(jj-1)*lb+ii-1];
};
// idem comme acces aux coordonnées bandes mais en constante
inline const double& cConst(int ii, int jj ) const
{
#ifdef MISE_AU_POINT
if (ii>lb || jj>dim )
{cout << "erreur d acces aux composantes bande";
cout << " ii = " << ii <<" jj = " << jj << '\n';
cout << "double& cConst(int ii, int jj )" << endl;
}
#endif
return pt[(jj-1)*lb+ii-1];
};
};
/// @} // end of group
#endif