// 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) .
//
// 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 .
//
// For more information, please consult: .
/************************************************************************
* DATE: 28/03/2005 *
* $ *
* AUTEUR: G RIO (mailto:gerardrio56@free.fr) *
* $ *
* PROJET: Herezh++ *
* $ *
************************************************************************
* BUT: Utilisation de routines lapack dans le cas de matrices *
* rectangulairesn, carrées ou bandes, symétriques ou pas. *
* Les coefficients sont des doubles. *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' * *
* VERIFICATION: *
* *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* ! ! ! ! *
* $ *
* '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' *
* MODIFICATIONS: *
* ! date ! auteur ! but ! *
* ------------------------------------------------------------ *
* $ *
************************************************************************/
#ifndef MatLapack_H
#define MatLapack_H
// au niveau de la bibliothèque apple
// ./System/Library/Frameworks/vecLib.framework/Versions/A/Headers/clapack.h
// ./System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/Headers/clapack.h
//#include "Debug.h"
#ifdef SYSTEM_MAC_OS_CARBON
#include "LAPACK.h"
#endif
#ifdef ENLINUX_2009
// #include // on supprime l'inclusion, et on la remplace par une introduction dans le debut de clapack.h !!
// ceci pour la version linux !!
// #include "INCLUDE/clapack.h"
#include "Include_lapack/clapack.h"
#include "cblas.h"
typedef __CLPK_integer ENTIER_POUR_CLAPACK;
// typedef long int ENTIER_POUR_CLAPACK;
// on n'arrive pas à faire un include en linux ??
// mais l'option de mise au point est supprimée à la fin du matlapack.cc
//#define MISE_AU_POINT
#else
#ifdef SYSTEM_MAC_OS_X_unix
// #include
// #include "Headers/clapack.h"
#include
// #include
typedef __CLPK_integer ENTIER_POUR_CLAPACK;
#endif
#endif
#ifdef SYSTEM_MAC_OS_X
#include
typedef __CLPK_integer ENTIER_POUR_CLAPACK;
#endif
#include "Mat_abstraite.h"
/// @addtogroup Les_classes_Matrices
/// @{
///
class MatLapack : public Mat_abstraite
{ // surcharge de l'operator de lecture typée
friend istream & operator >> (istream &, MatLapack &);
// surcharge de l'operator d'ecriture typée
friend ostream & operator << (ostream &, const MatLapack &);
public :
// Constructeur par défaut
// il est nécessaire ensuite de définir les données pour l'utiliser
MatLapack ( );
// Constructeur pour une matrice carree quelconque, ou symétrique
// se servant d'un nombre de lignes, de = colonnes, et eventuellement
// d'une valeur d'initialisation
MatLapack (int nb_ligne,int nb_colonne, bool symetrique , double val_init
,Enum_type_resolution_matri type_resol //=RIEN_TYPE_RESOLUTION_MATRI
,Enum_preconditionnement type_precondi //= RIEN_PRECONDITIONNEMENT
);
// Constructeur pour une matrice bande quelconque, ou symétrique
// se servant d'un nombre de lignes (dim) , d'une largeur de bande, et eventuellement
// d'une valeur d'initialisation et d'un type de résolution
MatLapack (Enum_matrice enu_mat, int lb_totale,int dim , bool symetrique, double val_init //=0.0
,Enum_type_resolution_matri type_resol//=GAUSS
,Enum_preconditionnement type_precondi//= RIEN_PRECONDITIONNEMENT
);
// de copie
MatLapack (const MatLapack& a);
// DESTRUCTEUR VIRTUEL :
~MatLapack ();
// METHODES DÉCOULANT DE 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
// ici s'il s'agit de même matrices de même place, c'est ok, sinon erreur
MatLapack & operator = ( const MatLapack &);
// 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 );
Mat_abstraite & operator = ( const Mat_abstraite &);
//-----------------------------------------------------------------------------------
// --- 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 MatLapack& mat_pl);
void operator+= (const Mat_abstraite& mat_pl);
// Surcharge de l'operateur -= : soustraction d'une matrice a la matrice courante
void operator-= (const MatLapack& mat_pl);
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 ;
int operator== (const MatLapack& mat_pl)const ;
// fonction permettant de creer une nouvelle instance d'element
// dérivé. la nouvelle instance = *this, il y a donc utilisation du constructeur de copie
Mat_abstraite * NouvelElement() const;
// Affichage des valeurs de la matrice
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 ;
// 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
void Change_Choix_resolution(Enum_type_resolution_matri type_resol,
Enum_preconditionnement type_precondi);
// Initialisation des valeurs des composantes
void Initialise (double val_init);
// Liberation de la place memoire
void Libere ();
// Retour du nombre de lignes
int Nb_ligne () const {return nb_lign;};
// Retour du nombre de colonnes
int Nb_colonne () const {return nb_col;};
// Symetrie de la matrice
int Symetrie () const;
// Acces aux valeurs de la matrices en écriture
// i : indice de ligne, j : indice de colonne
double& operator () (int i, int j);
// 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
double operator () ( int i, int j) const;
// 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
bool Existe(int i, int j) const;
// Retourne la ieme ligne de la matrice
// lorsque implemente ( a verifier )
Vecteur& Ligne_set(int i);
// 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 colonne fournie
void RemplaceColonneSpe(int j,const Vecteur & v);
//met une valeur identique sur toute la colonne
void MetValColonne(int j,double val);
// Resolution du systeme Ax=b
//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 Resol_syst
(const Tableau & 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 & Resol_systID
(Tableau & 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
// ici on calcul la factorisation LU de la matrice
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 Simple_Resol_syst
(const Tableau & 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 & Simple_Resol_systID
(Tableau & 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 colonne 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;
// ------------ méthode non virtuelle pur -> générales ----------------
//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 * V_Propres(MatLapack& KG,Tableau & VP );
// ================ Méthodes spécifiques à la classe MatLapack ====================
// changement de taille avec initialisation non nulle éventuelle
// si les tailles existantes sont identiques, on ne fait que initialiser
void Change_taille(int nb_li,int nb_col, const double& val_init = 0.);
// Surcharge de l'operateur + : addition de deux matrices
// il y a création d'une matrice de sortie
MatLapack operator+ (const MatLapack& mat_pl) const;
// Surcharge de l'operateur - : soustraction entre deux matrices
// il y a création d'une matrice de sortie
MatLapack operator- (const MatLapack& mat_pl) const;
// Surcharge de l'operateur - : oppose d'une matrice
// il y a création d'une matrice de sortie
MatLapack operator- () const;
// Surcharge de l'operateur * : multiplication de deux matrices
// il y a création d'une matrice de sortie
MatLapack operator* (const MatLapack& mat_pl) const;
// Surcharge de l'operateur * : multiplication d'une matrice par un scalaire
// il y a création d'une matrice de sortie
MatLapack operator* (const double& coeff) const ;
// Retourne la transposee de la matrice mat
// il y a création d'une matrice de sortie
MatLapack Transpose () const;
// determine l'inverse d'une matrice carre
MatLapack Inverse();
// idem mais en retournant l'inverse dans la matrice passée en paramètre
// qui doit être de même taille
MatLapack& Inverse(MatLapack& mat) ;
protected :
// VARIABLES PROTEGEES :
//__CLPK_integer: long int en 32bit et int en 64bit, définit dans clapack.h
ENTIER_POUR_CLAPACK nb_lign,nb_col; // nombre de ligne, nombre de colonne utilisées par herezh
ENTIER_POUR_CLAPACK ldb; // nombre de ligne du stockage effectif
// --- dans le cas d'une matrice bande, nb_col est la largeur de bande totale
// dans le cas d'une matrice bande symétrique, nb_col=1/2largeur de bande+1
// on utilise des variables intermédiaires:
ENTIER_POUR_CLAPACK kl,ku,lda; // correspond à la largeur en dessous et au dessus de la diagonale
//lda: largeur maxi du premier indice, qu'utilise réellement lapack: = 2*kl+ku+1 en bande non symétrique
// dans le cas d'une matrice rectangulaire = nb_lign
Vecteur mat; // la matrice vue comme un vecteur
// le stockage dépend du type de matrice considérée
//-- cas des matrices tridiagonales:
Vecteur B,C; // ligne sup et inf. La diagonale est stockée dans mat
//------- variables utilisées pour la résolution de systèmes avec lapack -----------
// ces variables sont indépendantes de la matrice donc ne rentre pas en ligne de compte dans les opérateurs
// par exemple (échange, sauvegarde ....)
ENTIER_POUR_CLAPACK * ipiv; // la matrice des pivots d'indices qui défini la matrice de permutation P
// utilisé par le prog clapack: dgesv_
MatLapack* ptB; // array intermédiaires contenant en entrée les seconds membres et en sortie les résultats
Vecteur travail; // espace de travail
// pointeur d'acces aux coordonnées avec i variant de 1 à nb_lign et j variant de 1 à nb_col
double& (MatLapack::*Coor) (int i, int j );
double (MatLapack::*Coor_const) (int i, int j ) const;
// ----- declaration des différents accès aux composantes
// 1) pour une matrice générale non sym, rectangulaire: stockage de type GE
double& Coor_GE(int i, int j)
#ifdef MISE_AU_POINT
{ return mat((j-1)*lda + i); };
#else
{ return mat.v[(j-1)*lda + i - 1]; };
#endif
double Coor_GE_const(int i, int j) const
#ifdef MISE_AU_POINT
{ return mat((j-1)*lda + i); };
#else
{ return mat.v[(j-1)*lda + i - 1]; };
#endif
// 2) pour une matrice bande générale: nb_col = la largeur de bande: stockage de type GB
double& Coor_GB(int i, int j)
#ifdef MISE_AU_POINT
{ return mat( (j-1)*lda + nb_col + i-j ); };
#else
{ return mat.v[( (j-1)*lda + nb_col + i-j ) - 1]; };
#endif
double Coor_GB_const(int i, int j) const
#ifdef MISE_AU_POINT
{ return mat( (j-1)*lda + nb_col + i-j ); };
#else
{ return mat.v[( (j-1)*lda + nb_col + i-j ) - 1]; };
#endif
// 3) pour une matrice bande symétrique: nb_col = la 1/2 largeur de bande: stockage de type PB
// c'est à dire symétrique positif, avec un stockage 'U'
// AB(ku+1+i-j,j) = A(i,j) pour max(1,j-ku) <= i <= j
double& Coor_PB(int i, int j)
#ifdef MISE_AU_POINT
{ if (j >= i) return mat( (j-1)*lda + nb_col + i-j);
else return mat( (i-1)*lda + nb_col + j-i);};
#else
{ if (j >= i) return mat.v[( (j-1)*lda + nb_col + i-j ) - 1];
else return mat.v[( (i-1)*lda + nb_col + j-i ) - 1];};
#endif
double Coor_PB_const(int i, int j) const
#ifdef MISE_AU_POINT
{ if (j >= i) return mat( (j-1)*lda + nb_col + i-j);
else return mat( (i-1)*lda + nb_col + j-i);};
#else
{ if (j >= i) return mat.v[( (j-1)*lda + nb_col + i-j ) - 1];
else return mat.v[( (i-1)*lda + nb_col + j-i ) - 1];};
#endif
// 4) pour une matrice carré symétrique: nb_col = la 1/2 largeur de bande: stockage de type PP
// c'est à dire symétrique positif, avec un stockage 'U' et stockage carré compacté
// AB(i+(j-1)*j/2) = A(i,j) pour 1 <= i <= j
double& Coor_PP(int i, int j)
#ifdef MISE_AU_POINT
{ if (j >= i) return mat( (j-1)*j/2 + i);
else return mat( (i-1)*i/2 + j);};
#else
{ if (j >= i) return mat.v[ (j-1)*j/2 + i - 1];
else return mat.v[ (i-1)*i/2 + j - 1];};
#endif
double Coor_PP_const(int i, int j) const
#ifdef MISE_AU_POINT
{ if (j >= i) return mat( (j-1)*j/2 + i);
else return mat( (i-1)*i/2 + j);};
#else
{ if (j >= i) return mat.v[ (j-1)*j/2 + i - 1];
else return mat.v[ (i-1)*i/2 + j - 1];};
#endif
// 5) pour une matrice tridiagonale quelconque : mat contient la diagonale
// B la diag sup, C la diag inf
// A(i-1,i) = B(i-1)
// A(i,i-1) = C(i-1); A(i,i) = mat(i); A(i,i+1) = B(i)
// A(i+1,i) = C(i)
// la suite = def de l'accès aux composantes: en chantier
// fonction interne de vérification des composantes
// arrêt si pb
void Verif_nb_composante(const Vecteur& b,string nom_routine) const ;
void Verif_nb_composante(const Tableau & b,string nom_routine) const ;
// gestion erreur dgesv et dgbsv
void GestionErreurDg_sv(const int& info) const;
// gestion erreur dpbsv
void GestionErreurDpbsv(const int& info) const;
// gestion erreur Inverse
void GestionErreurInverse(const int& info,int cas,const string nom_routine_lapack) const;
};
/// @} // end of group
// pour l'instant on n'arrive pas à faire des includes en linux ??
#ifndef MISE_AU_POINT
#undef MatLapack_H_deja_inclus
#include "MatLapack.cc"
#define MatLapack_H_deja_inclus
#endif
#endif