rio/Herezh_dev
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity
Herezh_dev/Resolin/Matrices/Mat_abstraite.cc

568 lines
23 KiB
C++
Raw Normal View History

ajout du répertoire Resolin
2021-09-28 17:28:23 +02:00
// 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-2021 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/>.
#include "Mat_abstraite.h"
#include "ConstMath.h"
#include <math.h>
#include "MathUtil.h"
#include <iomanip>
#include "CharUtil.h"
// partie relative à la résolution
#include "vecdefs_GR.h" // modif GR
//#include "mvblasd.h" // MV_Vector level 1 BLAS
#include "diagpre_double_GR.h" // Diagonal preconditioner
#include "icpre_double_GR.h" // preconditioner cholesky incomplet
#include "ilupre_double_GR.h" // preconditioner ILU
#include "cgs.h" // IML++ CGS template
#include "cg.h" // IML++ CG template
#include "bicg.h" // IML++ BiCG template
#include "cheby.h" // IML++ Cheby template
#include "bicgstab.h" // IML++ BiCGSTAB template
#include "gmres.h" // IML++ GMRES template
#include "ir.h" // IML++ IR template
#include "qmr.h" // IML++ QMR template
// grandeur par défaut déclarées avant car on s'en sert dans les passages de paramètres
const double Mat_abstraite::tol_defaut = 1.E-7; // tolérance sur le résidu dans les méthodes itératives
const int Mat_abstraite::maxit_defaut = 150 ; // maximum d'itération dans les méthodes itératives
const int Mat_abstraite::restart_defaut = 32;// maximum de restart itération (nb de vecteur sauvegardé)
//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> * Mat_abstraite::V_Propres(Mat_abstraite& KG,Tableau <VeurPropre> & VP )
{ // la méthode est décrite dans Touzot/Dhatt seconde édition p 393-395
#ifdef MISE_AU_POINT
// vérification des entrées
// 1) (*this) et KG sont carrées et identiques en taille
if ( (Nb_ligne() != KG.Nb_ligne()) || (Nb_colonne() != KG.Nb_colonne()) ||
(Nb_ligne() != Nb_colonne()))
{ cout << "\n ***** Erreur : la taille des matrices a traiter n\'est pas correcte !\n";
cout << "\n les matrices doivent être carrés !! \n";
cout << "Mat_abstraite::V_Propres(... \n";
cout << " KG CxL " << KG.Nb_colonne() << " x " << KG.Nb_ligne();
cout << " (*this) CxL : " << (Nb_colonne()) << " x " << (Nb_ligne()) << "\n";
Sortie(1);
};
#endif
// stockage inter
Vecteur Vinter(Nb_colonne());
bool sorIter ; // gestion de la sortie de la boucle sur les valeurs propres
// on boucle sur le nombre de valeur propre demandé
for (int il=1; il<= VP.Taille(); il++)
{ sorIter = true; // tout va bien à priori
// choix d'un vecteur initial
Vecteur V(Nb_colonne(),1.);
// calcul d'un décalage éventuelle du aux vecteurs déjà calculés
Vecteur Ybarre = KG.Prod_mat_vec(V);
for (int ill=1; ill< il; ill++)
V = V - (VP(ill).Pv() * Ybarre) * VP(ill).Pv();
// calcul initial de la sollicitation
Vecteur F = KG.Prod_mat_vec(V);
// initialisation en taille (uniquement) de la sollicitation intermédiaire
Vecteur Fbarre = F;
// des scalaires intermédiaires
double d,lambda = 0.;
double lambdaPrec = ConstMath::grand; // initialisation à une valeur grande
int maxiBoucle = 500; // maxi de boucle d'itération inverse
int nbBoucle = 0; // indice courant de nombre de boucle
while (Abs(lambda - lambdaPrec) > ConstMath::petit*MaX(Abs(lambda),Abs(lambdaPrec)))
{ // a chaque itération
nbBoucle++; // incrémentation du nombre de boucle
if (nbBoucle > maxiBoucle)
{cout << " \n PB dans le calcul de la valeur propre nb = "<< il;
cout << " Mat_abstraite::V_Propres(.. \n ";
cout << " lambda = "<< lambda;
sorIter = false; // pour l'arrêt de la boucle globale
break; // sortie des itérations
};
// sauvegarde de la raideur
Mat_abstraite * KSS = NouvelElement();
V = KSS->Resol_syst(F); // résolution
// calcul d'un décalage éventuelle du aux vecteurs déjà calculés
Ybarre = KG.Prod_mat_vec(V);
for (int ill=1; ill< il; ill++)
V = V - (VP(ill).Pv() * Ybarre) * VP(ill).Pv();
Fbarre = KG.Prod_mat_vec(V); // sollicitation intermédiaire
d = V * Fbarre; // facteur d
lambdaPrec = lambda; // sauvegarde
lambda = (V * F) / d; // approximation de lambda
F = Fbarre / sqrt(d); // nouvelle sollicitation
};
if (!sorIter) break; // sortie au cas ou l'on n'a pas convergé
// calcul du vecteur propre normalisé
Vinter = V / sqrt(d);
VP(il) = VeurPropre(lambda,Vinter); // sauvegarde
};
// retour du tableau
return &VP;
};
// définie la surcharge de multiplication d'une matrice par un MV_Vector
MV_Vector<double> Mat_abstraite::operator * (const MV_Vector<double> & vec) const
{ // tout d'abord on construit un vecteur en place du MV_Vector
const MV_Vecteur vecte(vec);
const Vecteur & vecv = (Vecteur &) vecte;
// multiplication et retour
// MV_Vecteur resulta = (Prod_mat_vec(vecv)).MV_vecteur_double();
//// MV_Vector<double> resulta = (Prod_mat_vec(vecv)).MV_vecteur_double();
// autre solution
MV_Vecteur resulta(vec);
Vecteur & vecretour = (Vecteur &) resulta;
Prod_mat_vec(vecv,vecretour);
return resulta;
};
// multiplication matrice transposée fois vecteur ce qui est équivalent à
// la transposée du résultat de : multiplication vecteur transposé fois matrice
MV_Vector<double> Mat_abstraite::trans_mult(const MV_Vector<double> &vec) const
{ // tout d'abord on construit un vecteur en place du MV_Vector
const MV_Vecteur vecte(vec);
const Vecteur & vecv = (Vecteur &) vecte;
// multiplication et retour
// MV_Vecteur resulta = (Prod_vec_mat (vecv)).MV_vecteur_double();
MV_Vector<double> resulta = (Prod_vec_mat (vecv)).MV_vecteur_double();
return resulta;
};
// -------------------------- 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
// en sortie : sol de même dimension que sol
void Mat_abstraite::Resolution_syst
(const Vecteur &b, Vecteur &sol,const double &tole, const int maxi, const int rest) const
{ // récup de l'homologue en vecteur_double
// const VECTOR_double bb_double(b.MV_vecteur_double());
const VECTOR_double * bb_double_pt = b.Nouveau_MV_Vector_double_const();
const VECTOR_double& bb_double = (*bb_double_pt);
// const VECTOR_double bb_double(b.Pointeur_vect(),b.Taille(),MV_Vector_::ref_type(1));
// MV_Vector<TYPE>::MV_Vector(TYPE* d, int n, MV_Vector_::ref_type i) : // modif GR
// idem pour la solution
// VECTOR_double sol_double (sol.MV_vecteur_double());
VECTOR_double* sol_double_pt = sol.Nouveau_MV_Vector_double();
VECTOR_double& sol_double = (*sol_double_pt);
int result,maxit = maxi,restart = rest; // Maximum iterations
double tol = tole;
// choix du préconditionnement
Pre_cond_double * D = NULL; // pointeur non affecté avec le choix
switch (type_preconditionnement)
{ case DIAGONAL :
// cas diagonal
D = new DiagPreconditioner_double (*this);
break;
case ICP :
// cas de la factorisation de cholesky incomplete
D = new ICPreconditioner_double (*this);
break;
case ILU :
// cas de la factorisation LU
D = new ILUPreconditioner_double (*this);
break;
default :
cout << "\nErreur : valeur incorrecte du type de preconditionnement !\n"
<< " = " << type_preconditionnement ;
cout << "\n Mat_abstraite::Resolution_syst (... ) \n";
Sortie(1);
};
if (ParaGlob::NiveauImpression() > 2) cout << "\n Resolution ";
// choix du type de résolution
switch (type_resolution)
{ case BI_CONJUG :
// biconjugué
{result = BiCG(*this, sol_double,bb_double, *D, maxit, tol);
if (ParaGlob::NiveauImpression() > 2)
cout << "BiCG flag = " << result << endl;
break;
}
case BI_CONJUG_STAB:
// biconjugué stabilisé
{result = BiCGSTAB( *this, sol_double, bb_double, *D, maxit, tol);
if (ParaGlob::NiveauImpression() > 2)
cout << "BiCGSTAB flag = " << result << endl;
break;
}
case CONJUG_GRAD :
// conjugué gradient
{result = CG(*this, sol_double, bb_double, *D, maxit, tol);
if (ParaGlob::NiveauImpression() > 2)
cout << "CG flag = " << result << endl;
break;
}
case CONJUG_GRAD_SQUARE :
// conjugué square method
{result = CGS(*this, sol_double, bb_double, *D, maxit, tol);
if (ParaGlob::NiveauImpression() > 2)
cout << "CGS flag = " << result << endl;
break;
}
case CHEBYSHEV :
// itération de chebyshev
{if (ParaGlob::NiveauImpression() > 2)
cout << "\n attention utilisation de bornes bidon pour le rayon spectrale"
<< " dans les itération de chebyshev (résolution de système matricielle)";
double mineig = .01; // n'importe quoi
double maxeig = 3; // devra être modifié pour avoir un résultat correcte
result = CHEBY(*this, sol_double, bb_double, *D, maxit, tol, mineig, maxeig);
if (ParaGlob::NiveauImpression() > 2)
cout << "CHEBY flag = " << result << endl;
break;
}
case GENE_MINI_RESIDUAL :
// résidu minimal généralisé
{MATRIX_double Htt(restart+1, restart, 0.0); // storage for upper Hessenberg H
result = GMRES(*this, sol_double, bb_double, *D, Htt,restart, maxit, tol);
if (ParaGlob::NiveauImpression() > 2)
cout << "GMRES flag = " << result << endl;
break;
}
case ITERATION_RICHARSON :
// itération de richarson
{result = IR( *this, sol_double, bb_double, *D, maxit, tol);
if (ParaGlob::NiveauImpression() > 2)
cout << "IR flag = " << result << endl;
break;
}
case QUASI_MINI_RESIDUAL :
// résidu quasi minimal
{result = QMR( *this, sol_double, bb_double, *D, *D, maxit, tol);
if (ParaGlob::NiveauImpression() > 2)
cout << "QMR flag = " << result << endl;
break;
}
default :
cout << "\nErreur : valeur incorrecte du type de résolution !\n"
<< " = " << type_resolution ;
cout << "\n Mat_abstraite::Resolution_syst (... ) \n";
Sortie(1);
};
if (ParaGlob::NiveauImpression() > 2)
{ cout << "iterations performed: " << maxit << endl;
cout << "tolerance achieved : " << tol << endl;
};
// on recopie le résultat
sol = sol_double;
// on efface le préconditionnement
if (D != NULL) delete D;
// on efface les vecteurs de travail
delete bb_double_pt; delete sol_double_pt;
};
// idem avec en entrée un tableau de vecteur second membre et
// en sortie un tableau de vecteurs solution
void Mat_abstraite::Resolution_syst
(const Tableau <Vecteur>& b, Tableau <Vecteur> &sol,const double &tole
, const int maxi, const int rest) const
{
// choix du préconditionnement
Pre_cond_double * D = NULL; // pointeur non affecté avec le choix
switch (type_preconditionnement)
{ case DIAGONAL :
// cas diagonal
D = new DiagPreconditioner_double (*this);
break;
case ICP :
// cas de la factorisation de cholesky incomplete
D = new ICPreconditioner_double (*this);
break;
case ILU :
// cas de la factorisation LU
D = new ILUPreconditioner_double (*this);
break;
default :
cout << "\nErreur : valeur incorrecte du type de preconditionnement !\n"
<< " = " << type_preconditionnement ;
cout << "\n Mat_abstraite::Resolution_syst (... ) \n";
Sortie(1);
};
// on boucle sur le nombre de second membre
int nbsecondmembre = b.Taille();
for (int ise=1; ise<= nbsecondmembre;ise++)
{
// récup de l'homologue en vecteur_double
const VECTOR_double bb_double ((b(ise)).MV_vecteur_double());
// idem pour la solution
VECTOR_double sol_double ((sol(ise)).MV_vecteur_double());
int result,maxit = maxi,restart = rest; // Maximum iterations
double tol = tole;
if (ParaGlob::NiveauImpression() > 2)
cout << "\n resolution du second membre nb = " << ise << " ";
// choix du type de résolution
switch (type_resolution)
{ case BI_CONJUG :
// biconjugué
{result = BiCG(*this, sol_double,bb_double, *D, maxit, tol);
if (ParaGlob::NiveauImpression() > 2)
cout << "BiCG flag = " << result << endl;
break;
}
case BI_CONJUG_STAB:
// biconjugué stabilisé
{result = BiCGSTAB( *this, sol_double, bb_double, *D, maxit, tol);
cout << "BiCGSTAB flag = " << result << endl;
break;
}
case CONJUG_GRAD :
// conjugué gradient
{result = CG(*this, sol_double, bb_double, *D, maxit, tol);
if (ParaGlob::NiveauImpression() > 2)
cout << "CG flag = " << result << endl;
break;
}
case CONJUG_GRAD_SQUARE :
// conjugué square method
{result = CGS(*this, sol_double, bb_double, *D, maxit, tol);
if (ParaGlob::NiveauImpression() > 2)
cout << "CGS flag = " << result << endl;
break;
}
case CHEBYSHEV :
// itération de chebyshev
{if (ParaGlob::NiveauImpression() > 2)
cout << "\n attention utilisation de bornes bidon pour le rayon spectrale"
<< " dans les iteration de chebyshev (resolution de systeme matricielle)";
double mineig = .01; // n'importe quoi
double maxeig = 3; // devra être modifié pour avoir un résultat correcte
result = CHEBY(*this, sol_double, bb_double, *D, maxit, tol, mineig, maxeig);
if (ParaGlob::NiveauImpression() > 2)
cout << "CHEBY flag = " << result << endl;
break;
}
case GENE_MINI_RESIDUAL :
// résidu minimal généralisé
{MATRIX_double Htt(restart+1, restart, 0.0); // storage for upper Hessenberg H
result = GMRES(*this, sol_double, bb_double, *D, Htt,restart, maxit, tol);
cout << "GMRES flag = " << result << endl;
break;
}
case ITERATION_RICHARSON :
// itération de richarson
{result = IR( *this, sol_double, bb_double, *D, maxit, tol);
if (ParaGlob::NiveauImpression() > 2)
cout << "IR flag = " << result << endl;
break;
}
case QUASI_MINI_RESIDUAL :
// résidu quasi minimal
{result = QMR( *this, sol_double, bb_double, *D, *D, maxit, tol);
if (ParaGlob::NiveauImpression() > 2)
cout << "QMR flag = " << result << endl;
break;
}
default :
cout << "\nErreur : valeur incorrecte du type de résolution !\n"
<< " = " << type_resolution ;
cout << "\n Mat_abstraite::Resolution_syst (... ) \n";
Sortie(1);
};
if (ParaGlob::NiveauImpression() > 2)
{cout << "iterations performed: " << maxit << endl;
cout << "tolerance achieved : " << tol << endl;
};
// on recopie le résultat
sol(ise) = sol_double;
};
// on efface le préconditionnement
if (D != NULL) delete D;
};
// 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 Mat_abstraite::Affiche1(int min_i,int max_i,int pas_i,int min_j,int max_j,int pas_j)
const
{ cout << "\n affichage d une matrice " << '\n';
int i,j;
cout << " ";
for (j=min_j;j<=max_j;j+=pas_j) cout << "col " << setw(5) << j << " " ;
cout << '\n';
for (i=min_i;i<=max_i;i+=pas_i)
{ cout << "lig " << setw(4) << i << " * ";
for (j=min_j;j<=max_j;j+=pas_j)
cout << setw(11) << setprecision(7) << (*this)(i,j) << " ";
cout << '\n';
};
cout << endl;
};
// 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
// avec en plus le nb de digit = nd
void Mat_abstraite::Affiche2
(int min_i,int max_i,int pas_i,int min_j,int max_j,int pas_j,int nd)
const
{ cout << "\n affichage d une matrice " << '\n';
if (nd < 7)
{ cout << " dimension des colonnes trop faible (<7 ! ) " << endl;
return;
};
int i,j;
SortBlanc(16);
for (j=min_j;j<=max_j;j+=pas_j)
{ cout << "col " << setw(2) << j ;
SortBlanc(nd-6);
};
cout << '\n';
for (i=min_i;i<=max_i;i+=pas_i)
{ cout << "li " << setw(3) << i << "* ";
for (j=min_j;j<=max_j;j+=pas_j)
cout << setw(nd) << setprecision(7) << (*this)(i,j);
cout << '\n';
};
cout << endl;
};
// affichage à l'écran après une demande interactive de la matrice
// entete: une chaine explicative, a afficher en entête
void Mat_abstraite::Affichage_ecran(string entete) const
{// affichage éventuelle de la matrice
cout << "\n --- " << entete << " -----";
string rep;
cout << "\n voulez-vous afficher la matrice ? (o ou n) ";
rep = lect_return_defaut(false,"n");
if (rep == "o")
{cout << "\n parametre d'affichage par defaut (toute la matrice) ? (o ou n) " ;
rep = lect_return_defaut(false,"n");
int min_i,max_i,pas_i,min_j,max_j,pas_j,nd;
if (rep =="o")
{ min_i=1; max_i=this->Nb_ligne(); min_j=1; max_j=this->Nb_ligne();
pas_i=1;pas_j=1; nd=14;
}
else
{try
{ cout << "\n a partir de quel numero de ligne ? (1) "; min_i=(int)lect_double();
cout << " jusqu'a quel numero de ligne ? ("<<this->Nb_ligne()<<") "; max_i=(int)lect_double();
cout << " par pas de ? (1) "; pas_i=(int)lect_double();
cout << "\n a partir de quel numero de colonne ? (1) "; min_j=(int)lect_double();
cout << " jusqu'a quel numero de colonne ? ("<<this->Nb_ligne()<<") "; max_j=(int)lect_double();
cout << " par pas de ? (1) "; pas_j=(int)lect_double();
cout << "\n nombre de diggit ? (14) "; nd=(int)lect_double();
}
catch (ErrSortieFinale)
// cas d'une direction voulue vers la sortie
// on relance l'interuption pour le niveau supérieur
{ ErrSortieFinale toto;
throw (toto);
}
catch (...)
{ cout << "\n Erreur sur l'entree des donnee !!, on prend les valeurs par defaut ";
min_i=1; max_i=this->Nb_ligne(); min_j=1; max_j=this->Nb_ligne();
pas_i=1;pas_j=1; nd=14;
}
};
this->Affiche2(min_i,max_i,pas_i,min_j,max_j,pas_j,nd);
cout << "\n entrez une lettre pour continuer ? ";
rep = lect_chaine();
}; //-- fin du choix visu ou pas de la matrice
};
// 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 Mat_abstraite::MinMaxMoy(bool affiche) const
{ Coordonnee3 retour;
// calcul des grandeurs
double min=ConstMath::tresgrand,max=-ConstMath::tresgrand,moy=0.;
int nbl = this->Nb_ligne();
if (affiche)
{for (int i=1;i<=nbl;i++)
{double dia = (*this)(i,i);
min = MiN(min, dia);
max = MaX(max, dia);
moy += Dabs(dia);
};
moy /= nbl;
// affichage
cout << "\n matrice masse: valeurs diagonales: min= "<<min <<", max= "<<max<<", moy des |(i,i)| "<<moy<<flush;
}
else
{for (int i=1;i<=nbl;i++)
{double dia = (*this)(i,i);
min = MiN(min, dia);
max = MaX(max, dia);
moy += Dabs(dia);
};
moy /= nbl;
};
retour(1) = min; retour(2) = max; retour(3)=moy;
return retour;
};
// 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 Mat_abstraite::Limitation_min_diag(double seuil_bas, double val_a_imposer)
{bool modif = false;
int nbl = this->Nb_ligne();
for (int i=1;i<=nbl;i++)
{double dia = (*this)(i,i);
if (dia < seuil_bas)
{ modif = true;
(*this)(i,i) = val_a_imposer;
};
};
return modif;
};
Reference in a new issue Copy permalink