532 lines
24 KiB
C++
532 lines
24 KiB
C++
|
|
// 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: 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 <f2c.h> // 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 </System/Library/Frameworks/vecLib.framework/Versions/A/Headers/clapack.h>
|
|
// #include "Headers/clapack.h"
|
|
#include <Accelerate.h>
|
|
|
|
// #include </System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Headers/clapack.h>
|
|
typedef __CLPK_integer ENTIER_POUR_CLAPACK;
|
|
#endif
|
|
#endif
|
|
|
|
#ifdef SYSTEM_MAC_OS_X
|
|
#include <Versions/A/Headers/clapack.h>
|
|
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 <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
|
|
// 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 <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 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 <VeurPropre>* V_Propres(MatLapack& KG,Tableau <VeurPropre> & 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 <Vecteur>& 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
|
|
|
|
|