3808 lines
148 KiB
C++
3808 lines
148 KiB
C++
// FICHIER : MatLapack.cc
|
|
// CLASSE : MatLapack
|
|
|
|
// 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/>.
|
|
|
|
|
|
# include <iostream>
|
|
using namespace std; //introduces namespace std
|
|
|
|
|
|
#include "MatLapack.h"
|
|
#include "Sortie.h"
|
|
#include "ConstMath.h"
|
|
#include "MathUtil.h"
|
|
#include "MatDiag.h"
|
|
#include <math.h>
|
|
#include <stdlib.h>
|
|
#include "CharUtil.h"
|
|
#include "ExceptionsMatrices.h"
|
|
#include "ParaGlob.h"
|
|
|
|
#ifdef ENLINUX_2009
|
|
#include "cblas.h"
|
|
#else
|
|
#ifdef SYSTEM_MAC_OS_X_unix
|
|
// #include </System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Headers/cblas.h>
|
|
// #include "Headers/cblas.h"
|
|
#include <Accelerate.h>
|
|
#endif
|
|
#endif
|
|
|
|
#ifdef SYSTEM_MAC_OS_X
|
|
#include <Versions/A/Headers/cblas.h>
|
|
#endif
|
|
|
|
#ifndef MatLapack_H_deja_inclus
|
|
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// il est nécessaire ensuite de définir les données pour l'utiliser
|
|
MatLapack::MatLapack ():
|
|
Mat_abstraite(RECTANGLE_LAPACK,GAUSS,RIEN_PRECONDITIONNEMENT)
|
|
,mat(),ipiv(NULL),ptB(NULL),travail(1)
|
|
,nb_lign(0),nb_col(0),ldb(1),kl(0),ku(0),lda(1) // seules les grandeurs pour le stockage rectangle sont bien initialisées
|
|
{ Coor = &MatLapack::Coor_GE; // association pour la récupération des coordonnées
|
|
Coor_const = &MatLapack::Coor_GE_const;
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Constructeur pour une matrice carree quelconque, ou symétrique
|
|
// se servant d'un nombre de lignes = colonnes, et eventuellement
|
|
// d'une valeur d'initialisation
|
|
MatLapack::MatLapack (int nb_l,int nb_c, bool symetrique, double val_init
|
|
,Enum_type_resolution_matri type_resol
|
|
,Enum_preconditionnement type_precondi):
|
|
Mat_abstraite(CARREE_LAPACK,type_resol,type_precondi)
|
|
,mat(0) // on le met à 0 au début
|
|
,ipiv(NULL),ptB(NULL),travail(1)
|
|
,nb_lign(nb_l),nb_col(nb_c),ldb(MaX(1,nb_l)),lda(MaX(1,nb_l))
|
|
,kl(0),ku(0) // ne servent pas
|
|
{ // on rectifie le type de matrice
|
|
if (symetrique)
|
|
{type_matrice = CARREE_SYMETRIQUE_LAPACK;
|
|
Coor = &MatLapack::Coor_PP; // association pour la récupération des coordonnées
|
|
Coor_const = &MatLapack::Coor_PP_const;
|
|
if (type_resol == RIEN_TYPE_RESOLUTION_MATRI)
|
|
type_resolution = CHOLESKY; // par défaut si rien n'est défini
|
|
mat.Change_taille(nb_l*(nb_c-1)/2+nb_c,val_init);
|
|
}
|
|
else
|
|
// sinon on n'est pas en symétrique
|
|
{ Coor = &MatLapack::Coor_GE; // association pour la récupération des coordonnées
|
|
Coor_const = &MatLapack::Coor_GE_const;
|
|
mat.Change_taille(nb_l*nb_c,val_init);
|
|
};
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Constructeur pour une matrice bande quelconque, ou symétrique
|
|
// se servant d'un nombre de lignes, de = colonnes, et eventuellement
|
|
// d'une valeur d'initialisation et d'un type de résolution
|
|
MatLapack::MatLapack (Enum_matrice enu_mat, int lb_totale,int dim ,bool symetrique, double val_init
|
|
,Enum_type_resolution_matri type_resol
|
|
,Enum_preconditionnement type_precondi
|
|
):
|
|
Mat_abstraite(enu_mat,type_resol,type_precondi)
|
|
,mat() ,nb_lign(dim),nb_col(lb_totale)
|
|
,ipiv(NULL),ptB(NULL),travail(1)
|
|
,ldb(MaX(1,dim)),kl(0),ku(0),lda(0)
|
|
{
|
|
if (enu_mat == BANDE_NON_SYMETRIQUE_LAPACK)
|
|
{// les colonnes de 1 à lb= (lb_totale-1)/2 sont utilisée par lapack
|
|
// pour les routines d'herezh, seules les valeurs à partir de la colonne 1+lb sont utilisées
|
|
// int lb = (lb_totale-1)/2;
|
|
kl=ku=((lb_totale-1)/2);
|
|
lda=2*kl+ku+1;
|
|
Coor = &MatLapack::Coor_GB; // association pour la récupération des coordonnées
|
|
Coor_const = &MatLapack::Coor_GB_const;
|
|
if (type_resolution == CHOLESKY)
|
|
#ifdef MISE_AU_POINT
|
|
{ cout << "\n ** attention, le type de resolution par defaut avec BANDE_NON_SYMETRIQUE_LAPACK "
|
|
<< " est GAUSS " ;
|
|
};
|
|
#endif
|
|
if ((type_resolution == RIEN_TYPE_RESOLUTION_MATRI)||(type_resolution == CHOLESKY))
|
|
type_resolution = GAUSS; // par défaut si rien n'est défini
|
|
}
|
|
else if (enu_mat == BANDE_SYMETRIQUE_LAPACK)
|
|
{// là on ne stocke que la partie symétrique, on suppose de plus que la matrice
|
|
// est définie positive, enfin il n'y a pas de stockage supplémentaire utilisé par lapack
|
|
// pour les routines d'herezh, les valeurs à partir de la colonne 1 sont utilisées
|
|
// int lb = lb_totale;
|
|
kl=ku=lb_totale-1;
|
|
lda=ku+1;
|
|
Coor = &MatLapack::Coor_PB; // association pour la récupération des coordonnées
|
|
Coor_const = &MatLapack::Coor_PB_const;
|
|
if (type_resolution == RIEN_TYPE_RESOLUTION_MATRI)
|
|
type_resolution = CHOLESKY; // par défaut si rien n'est défini
|
|
}
|
|
//*** en chantier
|
|
// else if (enu_mat == TRIDIAGONALE_GENE_LAPACK)
|
|
// {// matrice quelconque tri: là on stocke toute la matrice
|
|
// // est définie positive, enfin il n'y a pas de stockage supplémentaire utilisé par lapack
|
|
// // pour les routines d'herezh, les valeurs à partir de la colonne 1 sont utilisées
|
|
// // int lb = lb_totale;
|
|
// kl=ku=lb_totale-1;
|
|
// lda=ku+1;
|
|
// Coor = &MatLapack::Coor_PB; // association pour la récupération des coordonnées
|
|
// Coor_const = &MatLapack::Coor_PB_const;
|
|
// if (type_resolution == RIEN_TYPE_RESOLUTION_MATRI)
|
|
// type_resolution = CHOLESKY; // par défaut si rien n'est défini
|
|
// }
|
|
// else if (enu_mat == BANDE_SYMETRIQUE_LAPACK)
|
|
// {// là on ne stocke que la partie symétrique, on suppose de plus que la matrice
|
|
// // est définie positive, enfin il n'y a pas de stockage supplémentaire utilisé par lapack
|
|
// // pour les routines d'herezh, les valeurs à partir de la colonne 1 sont utilisées
|
|
// // int lb = lb_totale;
|
|
// kl=ku=lb_totale-1;
|
|
// lda=ku+1;
|
|
// Coor = &MatLapack::Coor_PB; // association pour la récupération des coordonnées
|
|
// Coor_const = &MatLapack::Coor_PB_const;
|
|
// if (type_resolution == RIEN_TYPE_RESOLUTION_MATRI)
|
|
// type_resolution = CHOLESKY; // par défaut si rien n'est défini
|
|
// }
|
|
//*** fin chantier
|
|
;
|
|
// def du conteneur de la matrice
|
|
mat.Change_taille(lda*ldb,val_init);
|
|
#ifdef MISE_AU_POINT
|
|
// pour l'instant on ne considère que des bandes d'étendues égales de part et autres de la diagonale
|
|
if ((enu_mat == BANDE_NON_SYMETRIQUE_LAPACK) && ((lb_totale % 2) == 0))
|
|
{ cout << "erreur de dimensionnement d'une matrice bande lapack ";
|
|
cout << " largeur de bande = " << lb_totale << " devrait etre un nombre impair " ;
|
|
cout << "\n MatLapack::MatLapack (... " << endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
#ifdef MISE_AU_POINT
|
|
// vérif que l'on est bien en bande lapack
|
|
if ((enu_mat != BANDE_SYMETRIQUE_LAPACK) && (enu_mat != BANDE_NON_SYMETRIQUE_LAPACK))
|
|
{ cout << "\n **** erreur: mauvais choix de matrice: " << Nom_matrice(enu_mat)
|
|
<< "\n MatLapack::MatLapack (Enum_matrice enu_mat,...";
|
|
Sortie(1);
|
|
};
|
|
if ((nb_lign < 0) || (nb_col < 0))
|
|
{ cout << "erreur de dimensionnement d une matrice bande";
|
|
cout << " LB = " << nb_col << " dim = " << nb_lign << '\n';
|
|
cout << "\n MatLapack::MatLapack (Enum_matrice enu_mat..." << endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
};
|
|
|
|
|
|
// Constructeur de copie
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
MatLapack::MatLapack (const MatLapack& mat_pl) :
|
|
Mat_abstraite(mat_pl)
|
|
,mat(mat_pl.mat)
|
|
,nb_lign(mat_pl.nb_lign),nb_col(mat_pl.nb_col)
|
|
,ipiv(NULL),ptB(NULL),travail(1)
|
|
,ldb(mat_pl.ldb),lda(mat_pl.lda),ku(mat_pl.ku),kl(mat_pl.kl)
|
|
{
|
|
if (type_matrice == CARREE_SYMETRIQUE_LAPACK)
|
|
{Coor = &MatLapack::Coor_PP; // association pour la récupération des coordonnées
|
|
Coor_const = &MatLapack::Coor_PP_const;
|
|
}
|
|
else if(type_matrice == CARREE_LAPACK)
|
|
// sinon on n'est pas en symétrique
|
|
{ Coor = &MatLapack::Coor_GE; // association pour la récupération des coordonnées
|
|
Coor_const = &MatLapack::Coor_GE_const;
|
|
}
|
|
else if (type_matrice == BANDE_NON_SYMETRIQUE_LAPACK)
|
|
{Coor = &MatLapack::Coor_GB; // association pour la récupération des coordonnées
|
|
Coor_const = &MatLapack::Coor_GB_const;
|
|
}
|
|
else if (type_matrice == BANDE_SYMETRIQUE_LAPACK)
|
|
{Coor = &MatLapack::Coor_PB; // association pour la récupération des coordonnées
|
|
Coor_const = &MatLapack::Coor_PB_const;
|
|
}
|
|
//*** en chantier
|
|
// else if (enu_mat == TRIDIAGONALE_GENE_LAPACK)
|
|
// {// matrice quelconque tri: là on stocke toute la matrice
|
|
// // est définie positive, enfin il n'y a pas de stockage supplémentaire utilisé par lapack
|
|
// // pour les routines d'herezh, les valeurs à partir de la colonne 1 sont utilisées
|
|
// // int lb = lb_totale;
|
|
// kl=ku=lb_totale-1;
|
|
// lda=ku+1;
|
|
// Coor = &MatLapack::Coor_PB; // association pour la récupération des coordonnées
|
|
// Coor_const = &MatLapack::Coor_PB_const;
|
|
// if (type_resolution == RIEN_TYPE_RESOLUTION_MATRI)
|
|
// type_resolution = CHOLESKY; // par défaut si rien n'est défini
|
|
// }
|
|
// else if (enu_mat == BANDE_SYMETRIQUE_LAPACK)
|
|
// {// là on ne stocke que la partie symétrique, on suppose de plus que la matrice
|
|
// // est définie positive, enfin il n'y a pas de stockage supplémentaire utilisé par lapack
|
|
// // pour les routines d'herezh, les valeurs à partir de la colonne 1 sont utilisées
|
|
// // int lb = lb_totale;
|
|
// kl=ku=lb_totale-1;
|
|
// lda=ku+1;
|
|
// Coor = &MatLapack::Coor_PB; // association pour la récupération des coordonnées
|
|
// Coor_const = &MatLapack::Coor_PB_const;
|
|
// if (type_resolution == RIEN_TYPE_RESOLUTION_MATRI)
|
|
// type_resolution = CHOLESKY; // par défaut si rien n'est défini
|
|
// }
|
|
//*** fin chantier
|
|
;
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
MatLapack::~MatLapack ()
|
|
// Destructeur
|
|
{ if (ipiv != NULL) delete [] ipiv;
|
|
if (ptB != NULL) delete ptB;
|
|
};
|
|
|
|
|
|
// fonction permettant de creer une nouvelle instance d'element
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
Mat_abstraite * MatLapack::NouvelElement() const
|
|
{ Mat_abstraite * a;
|
|
a = new MatLapack(*this);
|
|
return a;
|
|
};
|
|
|
|
|
|
// surcharge de l'opérateur d'affectation: cas des matrices abstraites
|
|
// il faut que les matrices soient ici de type et de tailles identiques sinon erreur
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
Mat_abstraite & MatLapack::operator = ( const Mat_abstraite & b)
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
// dans le cas où les matrices ne sont pas de caractéristiques identiques
|
|
// il y a un message d'erreur
|
|
if (type_matrice != b.Type_matrice())
|
|
{cout << "erreur dans l'operation d'affectation";
|
|
cout << " les matrices sont de types différents "
|
|
<< Nom_matrice(type_matrice) << " " << Nom_matrice(b.Type_matrice()) << '\n';
|
|
cout << "MatLapack::operator = ( const Mat_abstraite & b)" << endl;
|
|
Sortie(1);
|
|
}
|
|
#endif
|
|
const MatLapack & mat_pl = *((MatLapack*) & b);
|
|
(*this)=mat_pl;
|
|
// retour
|
|
return (*this);
|
|
};
|
|
|
|
// transfert des informations de *this dans la matrice passée en paramètre
|
|
// la matrice paramètre est au préalable, mise à 0.
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
void MatLapack::Transfert_vers_mat( Mat_abstraite & b )
|
|
{b.Initialise(0.); // init de la matrice
|
|
// puis on transfert
|
|
switch (type_matrice)
|
|
{case RECTANGLE_LAPACK: case CARREE_LAPACK: case CARREE_SYMETRIQUE_LAPACK:
|
|
{int le_max = nb_lign+1; // +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int i=1;i< le_max;i++)
|
|
{int le_maxj = nb_col+1; // +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int j=1;j< le_maxj;j++)
|
|
b(i,j)= (this->*Coor_const)(i,j);
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{int lb = ku;
|
|
int le_max = nb_lign+1; // +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int i=1;i< le_max;i++)
|
|
{int le_maxj = MiN(nb_lign,i+lb)+1;// +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int j= MaX(1,i-lb);j< le_maxj;j++)
|
|
b(i,j)=(*this)(i,j);
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{int lb = ku; // largeur de bande hors diagonale
|
|
int le_max = nb_lign+1; // +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int i=1;i< le_max;i++)
|
|
{int le_maxj = MiN(nb_lign,i+lb)+1;// +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int j= MaX(1,i-lb);j< le_maxj;j++)
|
|
b(i,j)=(*this)(i,j);
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::Affiche ()const "<<endl;
|
|
Sortie(1);
|
|
};
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
MatLapack& MatLapack::operator= ( const MatLapack& mat_pl)
|
|
// Surcharge de l'operateur = : affectation d'une matrice a une autre
|
|
// il faut que les matrices soient ici de tailles identiques sinon erreur
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
// si la place à affecter n'est pas identique -> erreur
|
|
if ((nb_lign != mat_pl.nb_lign) || (nb_col != mat_pl.nb_col))
|
|
{cout << "erreur dans l'operation d'affectation";
|
|
cout << " les matrices ont un nombre de lignes ou colonnes différents: nb ligne= "
|
|
<< nb_lign << " " << mat_pl.nb_lign << " nb colonne= " << nb_col << " " << mat_pl.nb_col << '\n';
|
|
cout << "MatLapack::operator = ( const MatLapack & )" << endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
// les deux matrices sont identiques en type et taille on affecte basiquement
|
|
// donc toutes les variables intermédiaires de dimensionnement sont identiques, reste le contenu de la matrice
|
|
mat = mat_pl.mat;
|
|
// en fait dans le cas d'un stockage bande, on applique également pour la partie travail
|
|
// qui sert pour la résolution uniquement, cependant, comme c'est difficile de séparer les grandeurs utiles (i,j)
|
|
// et celles de travail, on applique à toute la matrice (la partie travail sera de toute manière recalculée avant emploi
|
|
// retour
|
|
return (*this);
|
|
};
|
|
|
|
// Surcharge de l'operateur += : addition d'une matrice a la matrice courante
|
|
// la matrice à additionner doit être de même type ou diagonale, de même nombre
|
|
// de ligne
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
void MatLapack::operator+= (const Mat_abstraite& b)
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
// dans le cas où les matrices ne sont pas de caractéristiques identiques
|
|
// il y a un message d'erreur
|
|
if ((type_matrice != b.Type_matrice())&&(b.Type_matrice() != DIAGONALE))
|
|
{cout << " erreur les matrices sont de types différents "
|
|
<< Nom_matrice(type_matrice) << " " << Nom_matrice(b.Type_matrice()) << '\n';
|
|
cout << "MatLapack::operator+= (const Mat_abstraite& b)" << endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
if (b.Type_matrice() != DIAGONALE)
|
|
{ const MatLapack & mat_pl = *((MatLapack*) & b);
|
|
this->operator+= (mat_pl); // dans le cas de bande : cf. NBM operator= ( const MatLapack& mat_pl)
|
|
}
|
|
else // cas d'une matrice argument diagonale
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
// il faut que le nombre de ligne soit identique
|
|
if (nb_lign != b.Nb_ligne())
|
|
{cout << " les matrices ont un nombre de ligne différent "
|
|
<< nb_lign << " " << b.Nb_ligne() << '\n';
|
|
cout << "MatLapack::operator+= (const Mat_abstraite& b)" << endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
const MatDiag & mat_d = *((MatDiag*) & b);
|
|
int le_max = nb_lign+1;// +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int i=1;i< le_max;i++)
|
|
(*this)(i,i) += mat_d(i,i);
|
|
// mat += mat_d.Vecteur_MatDiag(); ?? une erreur 2 mai 2012
|
|
};
|
|
|
|
};
|
|
|
|
// Surcharge de l'operateur -= : soustraction d'une matrice a la matrice courante
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
void MatLapack::operator-= (const Mat_abstraite& b)
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
// dans le cas où les matrices ne sont pas de caractéristiques identiques
|
|
// il y a un message d'erreur
|
|
if ((type_matrice != b.Type_matrice())&&(b.Type_matrice() != DIAGONALE))
|
|
{cout << " erreur les matrices sont de types différents "
|
|
<< Nom_matrice(type_matrice) << " " << Nom_matrice(b.Type_matrice()) << '\n';
|
|
cout << "MatLapack::operator-= (const Mat_abstraite& b)" << endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
if (b.Type_matrice() != DIAGONALE)
|
|
{ const MatLapack & mat_pl = *((MatLapack*) & b);
|
|
this->operator-= (mat_pl); // dans le cas de bande : cf. NBM operator= ( const MatLapack& mat_pl)
|
|
}
|
|
else // cas d'une matrice argument diagonale
|
|
{ const MatDiag & mat_d = *((MatDiag*) & b);
|
|
#ifdef MISE_AU_POINT
|
|
// il faut que le nombre de ligne soit identique
|
|
if (nb_lign != b.Nb_ligne())
|
|
{cout << " les matrices ont un nombre de ligne différent "
|
|
<< nb_lign << " " << b.Nb_ligne() << '\n';
|
|
cout << "MatLapack::operator-= (const Mat_abstraite& b)" << endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
int le_max = nb_lign; // +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int i=1;i< le_max;i++)
|
|
(*this)(i,i) -= mat_d(i,i);
|
|
// mat -= mat_d.Vecteur_MatDiag(); ?? une erreur 2 mai 2012
|
|
};
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Addition a la matrice courante de la matrice mat_pl
|
|
void MatLapack::operator+= ( const MatLapack& mat_pl)
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
// il faut que les deux matrices de type lapack soient identiques
|
|
if (type_matrice != mat_pl.type_matrice)
|
|
{cout << " erreur les matrices sont de types différents "
|
|
<< Nom_matrice(type_matrice) << " " << Nom_matrice(mat_pl.type_matrice) << '\n';
|
|
cout << "MatLapack::operator+= (const MatLapack& )" << endl;
|
|
Sortie(1);
|
|
};
|
|
if ( ( nb_lign!=mat_pl.nb_lign )|| ( nb_col!=mat_pl.nb_col ) )
|
|
{cout << "Erreur : dimensions des matrices non egales\n";
|
|
cout << "MatLapack::OPERATOR+= (const MatLapack& )\n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
// les deux matrices sont identiques en type et taille on affecte basiquement
|
|
mat += mat_pl.mat; // dans le cas de bande : cf. NBM operator= ( const MatLapack& mat_pl)
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Soustrait a la matrice courante la matrice mat_pl
|
|
void MatLapack::operator-= ( const MatLapack& mat_pl)
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (type_matrice != mat_pl.type_matrice)
|
|
{cout << " erreur les matrices sont de types différents "
|
|
<< Nom_matrice(type_matrice) << " " << Nom_matrice(mat_pl.type_matrice) << '\n';
|
|
cout << "MatLapack::operator-= (const MatLapack& )" << endl;
|
|
Sortie(1);
|
|
};
|
|
if ( ( nb_lign!=mat_pl.nb_lign )|| ( nb_col!=mat_pl.nb_col ) )
|
|
{cout << "Erreur : dimensions des matrices non egales\n";
|
|
cout << "MatLapack::OPERATOR-= (const MatLapack& )\n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
// les deux matrices sont identiques en type et taille on affecte basiquement
|
|
mat -= mat_pl.mat; // dans le cas de bande : cf. NBM operator= ( const MatLapack& mat_pl)
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
//multiplication de la matrice courante par un scalaire
|
|
void MatLapack::operator*= (double coeff)
|
|
{ mat *= coeff; };
|
|
|
|
// Surcharge de l'operateur == : test d'egalite entre deux matrices
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
int MatLapack::operator== (const Mat_abstraite& b) const
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
// dans le cas où les matrices ne sont pas de caractéristiques identiques
|
|
// il y a un message d'erreur
|
|
if (type_matrice != b.Type_matrice())
|
|
{cout << " erreur les matrices sont de types différents "
|
|
<< Nom_matrice(type_matrice) << " " << Nom_matrice(b.Type_matrice()) << '\n';
|
|
cout << "MatLapack::operator+= (const Mat_abstraite& b)" << endl;
|
|
Sortie(1);
|
|
}
|
|
#endif
|
|
const MatLapack & mat_pl = *((MatLapack*) & b);
|
|
return this->operator==(mat_pl);
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Acces aux valeurs de la matrices en écriture
|
|
// i : indice de ligne, j : indice de colonne
|
|
double& MatLapack::operator () (int i, int j)
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK: case CARREE_SYMETRIQUE_LAPACK:
|
|
{if ((i < 1) || (i>nb_lign))
|
|
{cout << "erreur d acces aux composantes d une matrice rectangulaire";
|
|
cout << "\n erreur l'indice en i; " << i << " est impossible "
|
|
<< " car nb_ligne de la matrice = " << nb_lign
|
|
<< "\n double& MatLapack::operator () ( int i, int j)" << endl;
|
|
Sortie(1);
|
|
};
|
|
if ((j < 1) || (j > nb_col))
|
|
{cout << "erreur d acces aux composantes d une matrice rectangulaire";
|
|
cout << "\n erreur l'indice en j; " << j << " est impossible "
|
|
<< " car nb_col de la matrice = " << nb_col
|
|
<< "\n double& MatLapack::operator () ( int i, int j)" << endl;
|
|
Sortie(1);
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:case BANDE_SYMETRIQUE_LAPACK:
|
|
{ if (i<1 || (abs(i-j) > ku))
|
|
{cout << "erreur d acces aux composantes d une matrice bande lapack";
|
|
cout << " indice hors borne (" << i << "," << j <<"), lb = " << ku;
|
|
cout << ", condition inacceptable: i<1 || abs(i-j)="<< abs(i-j)<< " >= lb " << '\n';
|
|
cout << "double& MatLapack::operator () (int i, int j )" << endl;
|
|
};
|
|
if (j<1 || j>nb_lign)
|
|
{cout << "erreur d acces aux composantes d une matrice bande";
|
|
cout << " second indice hors borne j = " << j <<" dim = " << nb_lign << '\n';
|
|
cout << "double& MatLapack::operator () (int i, int j )" << endl;
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::operator () (int i, int j) "<<endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
return (this->*Coor)(i,j);
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// 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 MatLapack::operator () ( int i, int j) const
|
|
{double toto;
|
|
#ifdef MISE_AU_POINT
|
|
if ((i < 1) || (i>nb_lign))
|
|
{cout << "erreur d acces aux composantes d une matrice ";
|
|
cout << "\n erreur l'indice en i; " << i << " est impossible "
|
|
<< " car nb_ligne de la matrice = " << nb_lign
|
|
<< "\n MatLapack::operator () ( int i, int j) const " << endl;
|
|
Sortie(1);
|
|
};
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK: case CARREE_SYMETRIQUE_LAPACK:
|
|
{if ((j < 1) || (j > nb_col))
|
|
{cout << "erreur d acces aux composantes d une matrice ";
|
|
cout << "\n erreur l'indice en j; " << j << " est impossible "
|
|
<< " car nb_col de la matrice = " << nb_col
|
|
<< "\n MatLapack::operator () ( int i, int j) const " << endl;
|
|
Sortie(1);
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK: case BANDE_SYMETRIQUE_LAPACK:
|
|
// {if (i<1 || (abs(i-j) >= (nb_col-1)/2 +1))
|
|
{if (i<1)
|
|
{cout << "erreur d acces aux composantes d une matrice bande lapack";
|
|
cout << " indice hors borne (" << i << "," << j <<"), lb = " << ku;
|
|
cout << ", condition acceptable: i<1 || abs(i-j) >= 2*lb " << '\n';
|
|
cout << "MatLapack::operator () ( int i, int j) const" << endl;
|
|
Sortie(1);
|
|
};
|
|
if (j<1 || j>nb_lign)
|
|
{cout << "erreur d acces aux composantes d une matrice bande";
|
|
cout << " second indice hors borne j = " << j <<" dim = " << nb_lign << '\n';
|
|
cout << "MatLapack::operator () ( int i, int j) const" << endl;
|
|
Sortie(1);
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::operator () (int i, int j)const "<<endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK: case CARREE_SYMETRIQUE_LAPACK:
|
|
{ return (this->*Coor_const)(i,j); break;}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK: case BANDE_SYMETRIQUE_LAPACK:
|
|
{ if (abs(i-j) >= 2*kl) { return 0;}
|
|
else {return (this->*Coor_const)(i,j);};
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::operator () (int i, int j)const "<<endl;
|
|
Sortie(1);
|
|
};
|
|
return toto; // pour taire le compilo
|
|
};
|
|
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// 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 MatLapack::Existe(int i, int j) const
|
|
{switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK: return true;break;
|
|
case CARREE_SYMETRIQUE_LAPACK:
|
|
return ((i-j)<1);
|
|
// erreur: corr 23 sep 2014 case BANDE_NON_SYMETRIQUE_LAPACK: return ((abs(i-j)-(nb_col-1)/2) < 0);break;
|
|
case BANDE_NON_SYMETRIQUE_LAPACK: return ((abs(i-j)-(nb_col-1)/2) < 1);break;
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
return ( ((i-j)<1) && ((i-j)<(ku+1)) );
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::Existe(.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
return false; // pour traire le compilo
|
|
};
|
|
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
void MatLapack::Affiche () const
|
|
// Affiche les donnees de la matrice : dimensions, composantes
|
|
{switch (type_matrice)
|
|
{case RECTANGLE_LAPACK: case CARREE_LAPACK: case CARREE_SYMETRIQUE_LAPACK:
|
|
{cout << "\tNombre de lignes : " << Nb_ligne() << "\n";
|
|
cout << "\tNombre de colonnes : " << Nb_colonne() << "\n";
|
|
int le_max = nb_lign+1; // +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int i=1;i< le_max;i++)
|
|
{cout << "Ligne " << i << " :\t{\t";
|
|
int le_maxj = nb_col+1; // +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int j=1;j< le_maxj;j++)
|
|
cout << (this->*Coor_const)(i,j) << "\t";
|
|
cout << "}\n";
|
|
};
|
|
cout << "\n";
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{cout << "\n affichage d une matrice bande non symetrique de largeur= "
|
|
<< (2*ku+1) << " dimension= " << nb_lign << '\n';
|
|
int lb = ku;
|
|
int le_max = nb_lign+1; // +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int i=1;i< le_max;i++)
|
|
{int le_maxj = MiN(nb_lign,i+lb)+1;// +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int j= MaX(1,i-lb);j< le_maxj;j++)
|
|
{cout << setw(4) << " i= " << i << " j= " << j << " * "
|
|
<< setw(14) << setprecision(7) << (*this)(i,j) << '\n';
|
|
};
|
|
};
|
|
cout << endl ;
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{cout << "\n affichage d une matrice bande symetrique de largeur= "
|
|
<< 2*ku+1 << " dimension= " << nb_lign << '\n';
|
|
int lb = ku; // largeur de bande hors diagonale
|
|
int le_max = nb_lign+1; // +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int i=1;i< le_max;i++)
|
|
{int le_maxj = MiN(nb_lign,i+lb)+1;// +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int j= MaX(1,i-lb);j< le_maxj;j++)
|
|
{cout << setw(4) << " i= " << i << " j= " << j << " * "
|
|
<< setw(14) << setprecision(7) << (*this)(i,j) << '\n';
|
|
};
|
|
};
|
|
cout << endl ;
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::Affiche ()const "<<endl;
|
|
Sortie(1);
|
|
};
|
|
};
|
|
|
|
// Affiche une partie de la matrice (utile pour le debug)
|
|
// MiN_ et max_ sont les bornes de la sous_matrice
|
|
// pas_ indique le pas en i et j pour les indices
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
void MatLapack::Affiche1(int min_i,int max_i,int pas_i,int min_j,int max_j,int pas_j) const
|
|
{switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:case CARREE_SYMETRIQUE_LAPACK:
|
|
{ Mat_abstraite::Affiche1(min_i,max_i,pas_i,min_j,max_j,pas_j);
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ cout << "\n affichage d une matrice bande non symetrique de largeur= " <<
|
|
nb_col << " dimension= " << nb_lign << '\n';
|
|
int lb = (nb_col-1)/2; // largeur de bande hors diagonale
|
|
cout << " ";
|
|
for (int j=min_j;j<=max_j;j+=pas_j) cout << "col " << setw(4) << j << " " ;
|
|
cout << '\n';
|
|
for (int i=min_i;i<=max_i;i+=pas_i)
|
|
{ cout << "lig " << setw(4) << i << " * ";
|
|
for (int j=min_j;j<=max_j;j+=pas_j)
|
|
if(abs(i-j)<=lb)
|
|
cout << setw(11) << setprecision(7) << (*this)(i,j);
|
|
else
|
|
cout << "...........";
|
|
cout << '\n';
|
|
};
|
|
cout << endl;
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{ cout << "\n affichage d une matrice bande symetrique de largeur= " <<
|
|
nb_col << " dimension= " << nb_lign << '\n';
|
|
int lb = (nb_col-1); // largeur de bande hors diagonale
|
|
cout << " ";
|
|
for (int j=min_j;j<=max_j;j+=pas_j) cout << "col " << setw(4) << j << " " ;
|
|
cout << '\n';
|
|
for (int i=min_i;i<=max_i;i+=pas_i)
|
|
{ cout << "lig " << setw(4) << i << " * ";
|
|
for (int j=min_j;j<=max_j;j+=pas_j)
|
|
if(abs(i-j)<=lb)
|
|
cout << setw(11) << setprecision(7) << (*this)(i,j);
|
|
else
|
|
cout << "...........";
|
|
cout << '\n';
|
|
};
|
|
cout << endl;
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::Affiche1(.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
};
|
|
// Affiche une partie de la matrice (utile 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
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
void MatLapack::Affiche2
|
|
(int min_i,int max_i,int pas_i,int min_j,int max_j,int pas_j,int nd) const
|
|
{switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:case CARREE_SYMETRIQUE_LAPACK:
|
|
{ Mat_abstraite::Affiche2(min_i,max_i,pas_i,min_j,max_j,pas_j,nd);
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK: case BANDE_SYMETRIQUE_LAPACK:
|
|
{ cout << "\n affichage d une matrice bande de largeur= "
|
|
<< nb_col << " dimension= " << nb_lign << '\n';
|
|
#ifdef MISE_AU_POINT
|
|
if (nd < 7)
|
|
{ cout << " dimension des colonnes trop faible (<7 ! ) " << endl;
|
|
return;
|
|
};
|
|
#endif
|
|
int lb = ku; // largeur de bande hors diagonale
|
|
SortBlanc(7+9);
|
|
|
|
for (int j=min_j;j<=max_j;j+=pas_j)
|
|
{ cout << "col " << setw(2) << j ;
|
|
SortBlanc(nd-6);
|
|
};
|
|
cout << '\n';
|
|
for (int i=min_i;i<=max_i;i+=pas_i)
|
|
{ cout << "li " << setw(3) << i << "* ";
|
|
for (int j=min_j;j<=max_j;j+=pas_j)
|
|
if(abs(i-j)<=lb)
|
|
cout << setw(nd) << setprecision(7) << (*this)(i,j);
|
|
else
|
|
Sortcar(nd,'.');
|
|
cout << '\n';
|
|
};
|
|
cout << endl;
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::Affiche2 ()const "<<endl;
|
|
Sortie(1);
|
|
};
|
|
};
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// 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 MatLapack::Change_Choix_resolution(Enum_type_resolution_matri type_resol,
|
|
Enum_preconditionnement type_precondi)
|
|
{
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ // pour l'instant seul GAUSS et CHOLESKY sont possibles !! non !!
|
|
// if ((type_resol != GAUSS) && (type_resol != CHOLESKY))
|
|
// { cout << "\n ** attention le type de resolution " << Nom_resolution(type_resol)
|
|
// << " n'est pas possible avec la matrice " << Nom_matrice(type_matrice)
|
|
// << " on change en la valeur par defaut : GAUSS ";
|
|
// type_resolution = GAUSS; type_preconditionnement = RIEN_PRECONDITIONNEMENT;
|
|
// }
|
|
// else
|
|
{ type_resolution = type_resol; type_preconditionnement = type_precondi; };
|
|
break;
|
|
}
|
|
case CARREE_SYMETRIQUE_LAPACK:
|
|
{ if (type_resol == GAUSS) // gauss n'est pas possible
|
|
{ cout << "\n ** attention le type de resolution " << Nom_resolution(type_resol)
|
|
<< " n'est pas possible avec la matrice " << Nom_matrice(type_matrice)
|
|
<< " on change en la valeur par defaut : CHOLESKY ";
|
|
type_resolution = CHOLESKY; type_preconditionnement = RIEN_PRECONDITIONNEMENT;
|
|
}
|
|
else
|
|
{ type_resolution = type_resol; type_preconditionnement = type_precondi; };
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ // CHOLESKY n'est pas possible
|
|
if (type_resol == CHOLESKY)
|
|
{ cout << "\n ** attention le type de resolution " << Nom_resolution(type_resol)
|
|
<< " n'est pas possible avec la matrice " << Nom_matrice(type_matrice)
|
|
<< " on change en la valeur par defaut : GAUSS ";
|
|
type_resolution = GAUSS; type_preconditionnement = RIEN_PRECONDITIONNEMENT;
|
|
}
|
|
else
|
|
{ type_resolution = type_resol; type_preconditionnement = type_precondi; };
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{ // GAUSS n'est pas possible
|
|
if (type_resol == GAUSS)
|
|
{ cout << "\n ** attention le type de resolution " << Nom_resolution(type_resol)
|
|
<< " n'est pas possible avec la matrice " << Nom_matrice(type_matrice)
|
|
<< " on change en la valeur par defaut : CHOLESKY ";
|
|
type_resolution = CHOLESKY; type_preconditionnement = RIEN_PRECONDITIONNEMENT;
|
|
}
|
|
else
|
|
{ type_resolution = type_resol; type_preconditionnement = type_precondi; };
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n **** erreur: mauvais choix de matrice: " << Nom_matrice(type_matrice)
|
|
<< "\n MatLapack::Change_Choix_resolution(...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
void MatLapack::Initialise (double val_init)
|
|
// Initialise les composantes de la matrice a val_init
|
|
{ // initialisation
|
|
mat.Inita(val_init); // dans le cas de bande : voir NB operator= ( const MatLapack& mat_pl)
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Liberation de la place memoire
|
|
void MatLapack::Libere ()
|
|
{ mat.Libere();
|
|
nb_lign=nb_col=0;
|
|
if (ipiv != NULL) {delete [] ipiv;ipiv=NULL;};
|
|
if (ptB != NULL) {delete ptB ; ptB=NULL;};//->Libere();
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// changement de taille avec initialisation non nulle éventuelle
|
|
// si les tailles existantes sont identiques, on ne fait que initialiser
|
|
void MatLapack::Change_taille(int nb_li,int nb_co, const double& val_init)
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if ((nb_li < 0) || (nb_co < 0))
|
|
{cout << "erreur de dimensionnement d une matrice ";
|
|
cout << " nb_li = " << nb_li << " nb_co = " << nb_co << '\n';
|
|
cout << "void MatLapack::Change_taille(... " << endl;
|
|
}
|
|
#endif
|
|
if ((nb_li != nb_lign) || (nb_co != nb_col))
|
|
{ ldb = MaX(1,nb_li);
|
|
this->Libere();
|
|
nb_lign = nb_li; nb_col = nb_co;
|
|
bool matrice_a_dimensionner = false; // init
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{// on change de taille
|
|
lda= MaX(1,nb_lign);
|
|
kl=ku=0;
|
|
matrice_a_dimensionner=true;
|
|
break;
|
|
}
|
|
case CARREE_SYMETRIQUE_LAPACK:
|
|
{// on change de taille, on utilise un stockage compacté, qui utilise la moitié de la taille théorique
|
|
lda= MaX(1,nb_lign);
|
|
kl=ku=0;
|
|
mat.Change_taille(nb_lign*(nb_col+1)/2,val_init); // division entière !
|
|
matrice_a_dimensionner=false; // pour signaler que le dimensionnement est déjà fait
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{// calcul des grandeurs relatives à la largeur de bande
|
|
#ifdef MISE_AU_POINT
|
|
// pour l'instant on ne considère que des bandes d'étendues symétriques de part et autres de la diagonale
|
|
if ((nb_co % 2) == 0)
|
|
{ cout << "erreur de dimensionnement d'une matrice bande non symetrique lapack ";
|
|
cout << " largeur de bande = " << nb_co << " devrait etre un nombre impair " ;
|
|
cout << "\n MatLapack::Change_taille(... " << endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
kl = (nb_col-1)/2; ku=kl;
|
|
lda = 2*kl+ku+1;
|
|
matrice_a_dimensionner=true;
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{// là on ne stocke que la partie symétrique, on suppose de plus que la matrice
|
|
// est définie positive, enfin il n'y a pas de stockage supplémentaire utilisé par lapack
|
|
// pour les routines d'herezh, les valeurs à partir de la colonne 1 sont utilisées
|
|
// int lb = lb_totale;
|
|
kl = ku = nb_col-1;
|
|
lda = nb_col;
|
|
matrice_a_dimensionner=true;
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::Change_taille(.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
// dimensionnement conditionnel
|
|
if (matrice_a_dimensionner)
|
|
mat.Change_taille(lda*ldb,val_init);
|
|
}
|
|
else
|
|
// sinon c'est seulement un dimensionnement qui est demandé
|
|
{ mat.Inita(val_init); }; // dans le cas de bande : voir NB operator= ( const MatLapack& mat_pl)
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Test sur la symetrie de la matrice
|
|
// Retourne 1 si la matrice est symetrique
|
|
// 0 sinon
|
|
int MatLapack::Symetrie () const
|
|
{ switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ if ( nb_lign != nb_col )
|
|
return 0;
|
|
for (int i=1;i<=nb_lign;i++)
|
|
for (int j=1;j<i;j++)
|
|
if ( diffpetit((this->*Coor_const)(i,j),(this->*Coor_const)(j,i))) // test d'egalite de la ieme jieme composante
|
|
return 0;
|
|
return 1;
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ int lb = (nb_col-1)/2;
|
|
for (int i=1;i<=nb_lign;i++)
|
|
for (int j=1;j<i;j++)
|
|
if (abs(i-j) <= lb)
|
|
if ( diffpetit((this->*Coor_const)(i,j),(this->*Coor_const)(j,i))) // test d'egalite de la ieme jieme composante
|
|
return 0;
|
|
return 1;
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:case CARREE_SYMETRIQUE_LAPACK:
|
|
{ return 1; // par construction la matrice est symétrique
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::Symetrie(.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
return 0; // pour taire le compilo
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Retourne la ieme ligne de la matrice
|
|
// lorsque implemente ( a verifier )
|
|
Vecteur& MatLapack::Ligne_set(int )
|
|
{ cout << "\n erreur methode non implante "
|
|
<< "\n MatLapack::Ligne_set(int i)";
|
|
Sortie(1);
|
|
// on passe jamais dans la suite !
|
|
return mat; // pour taire le compilo
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Retourne la ieme ligne de la matrice
|
|
Vecteur MatLapack::Ligne(int i) const
|
|
{ Vecteur ret(nb_lign); // init d'un vecteur total (indépendament du stockage)
|
|
// nb_lign = la dimension de la matrice
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:case CARREE_SYMETRIQUE_LAPACK:
|
|
{ int le_max = nb_col+1; // +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int j=1;j< le_max;j++)
|
|
ret(j)=(this->*Coor_const)(i,j);
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{// nb_lign = la dimension de la matrice théorique carrée
|
|
int lb=(nb_col-1)/2; int dim = nb_lign;
|
|
int maxJ = MiN(dim,i+lb)+1;// +1 pour inférieur pure sur la ligne qui suit
|
|
for (int j= MaX(1,i-lb);j< maxJ;j++)
|
|
ret(j) = (this->*Coor_const)(i,j);
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{// nb_lign = la dimension de la matrice théorique carrée
|
|
int lb=ku;
|
|
int maxJ=MiN(nb_lign,i+lb)+1;// +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int j= MaX(1,i-lb);j< maxJ;j++)
|
|
ret(j) = (this->*Coor_const)(i,j);
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::Ligne(.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
return ret;
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Retourne la ieme ligne de la matrice
|
|
// sous le format de stokage propre a la matrice
|
|
// donc a n'utiliser que comme sauvegarde en parallele
|
|
// avec la fonction RemplaceLigne
|
|
Vecteur MatLapack::LigneSpe(int i) const
|
|
{ switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{Vecteur ret(nb_col);
|
|
for (int j=1;j<=nb_col;j++)
|
|
ret(j)=(this->*Coor_const)(i,j);
|
|
return ret;
|
|
break;
|
|
}
|
|
case CARREE_SYMETRIQUE_LAPACK:
|
|
{// comme la ligne adresse deux mêmes parties symétriques, il vaut mieux
|
|
// sauvegarder toute la ligne
|
|
Vecteur ret(nb_col);
|
|
for (int j=1 ;j<=nb_col;j++)
|
|
ret(j)=(this->*Coor_const)(i,j);
|
|
return ret;
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{Vecteur ligne(nb_col); // mise a zero de ligne
|
|
int lb=(nb_col-1)/2; int dim = nb_lign;
|
|
int maxJ = MiN(dim,i+lb)+1;// +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int j= MaX(1,i-lb);j< maxJ;j++)
|
|
ligne(i-j+lb+1) = (this->*Coor_const)(i,j);
|
|
return ligne;
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{// comme la ligne adresse deux mêmes parties symétriques, il vaut mieux
|
|
// sauvegarder toute la ligne
|
|
int lb=nb_col-1; int dim = nb_lign;
|
|
Vecteur ligne(2*lb+1); // mise a zero de ligne
|
|
int le_max = MiN(dim,i+lb)+1;// +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
int le_min = MaX(1,i-lb);
|
|
for (int j= le_min;j < le_max;j++)
|
|
// for (int j= MaX(1,i);j<= MiN(dim,i+lb);j++)
|
|
ligne(i-j+lb+1) = (this->*Coor_const)(i,j);
|
|
return ligne;
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::LigneSpe(.. "<<endl;
|
|
Sortie(1);
|
|
Vecteur ret; // pour taire le compilo
|
|
return ret; // "
|
|
};
|
|
Vecteur ret; // pour taire le compilo
|
|
return ret; // "
|
|
};
|
|
|
|
// remplace la ligne de la matrice par la ligne fournie
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
void MatLapack::RemplaceLigneSpe(int i,const Vecteur & v)
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{if (nb_col != v.Taille())
|
|
{cout << " \nla taille de la ligne i= " << i <<" a remplacer n est pas bonne " ;
|
|
cout << " taille = " << v.Taille() << " tailleLigneMatrice = " << nb_col;
|
|
cout << " void MatLapack::RemplaceLigne(int i,Vecteur & v)" << endl;
|
|
Sortie (1);
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{if (v.Taille() != nb_col)
|
|
{cout << "\nerreur d affectation pour le vecteur";
|
|
cout << " dim vecteur = " << v.Taille() << " au lieu de " << nb_col <<'\n';
|
|
cout << "void MatLapack::RemplaceLigne(int i,Vecteur & v)" << endl;
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{if (v.Taille() != nb_col)
|
|
{cout << "\nerreur d affectation pour le vecteur";
|
|
cout << " dim vecteur = " << v.Taille() << " au lieu de " << nb_col <<'\n';
|
|
cout << "void MatLapack::RemplaceLigne(int i,Vecteur & v)" << endl;
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::RemplaceLigneSpe(.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ for(int j=1; j <= nb_col; j++)
|
|
(this->*Coor)(i,j)=v(j);
|
|
break;
|
|
}
|
|
case CARREE_SYMETRIQUE_LAPACK:
|
|
// comme la ligne adresse deux mêmes parties symétriques, il vaut mieux
|
|
// sauvegarder toute la ligne
|
|
{for (int j=1 ;j<=nb_col;j++)
|
|
(this->*Coor)(i,j)=v(j);
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{int lb=(nb_col-1)/2; int dim = nb_lign;
|
|
int maxJ = MiN(dim,i+lb)+1;// +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int j= MaX(1,i-lb);j< maxJ;j++)
|
|
(this->*Coor)(i,j)=v(i-j+lb+1);
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{// comme la ligne adresse deux mêmes parties symétriques, il vaut mieux
|
|
// sauvegarder toute la ligne
|
|
int lb=nb_col-1; int dim = nb_lign;
|
|
Vecteur ligne(2*lb+1); // mise a zero de ligne
|
|
int le_max = MiN(dim,i+lb)+1;// +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
int le_min = MaX(1,i-lb);
|
|
for (int j= le_min;j< le_max;j++)
|
|
(this->*Coor)(i,j)=v(i-j+lb+1);
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::RemplaceLigneSpe(.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Resolution du systeme Ax=b
|
|
Vecteur MatLapack::Resol_syst(const Vecteur& b,const double &tole, const int maxi, const int rest)
|
|
{// init général
|
|
ENTIER_POUR_CLAPACK nb_SM = 1; // un seul second membre
|
|
if (ipiv == NULL) {ipiv = new ENTIER_POUR_CLAPACK [nb_lign];};
|
|
ENTIER_POUR_CLAPACK info = 0; // info de sortie en fait c'est le retour de la fonction
|
|
Vecteur resul(b); // Vecteur solution initialisé = second membre
|
|
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case GAUSS:
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
Verif_nb_composante(b,"Resol_syst");
|
|
#endif
|
|
// appel de la routine clapack
|
|
// /* Subroutine */ int dgesv_(__CLPK_integer *n, __CLPK_integer *nrhs, __CLPK_doublereal *a, __CLPK_integer
|
|
// *lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer *ldb, __CLPK_integer *info);
|
|
dgesv_(&nb_lign,&nb_SM,mat.v,&lda,ipiv,resul.v,&ldb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{ Vecteur resul(b);
|
|
Mat_abstraite::Resolution_syst(b,resul,tole,maxi,rest);
|
|
return resul;
|
|
}
|
|
};
|
|
break;
|
|
}
|
|
case CARREE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case CHOLESKY:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dppsv_(char *uplo, __CLPK_integer *n, __CLPK_integer *nrhs,
|
|
// __CLPK_doublereal *ap, __CLPK_doublereal *b, __CLPK_integer *ldb, __CLPK_integer *info)
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
dppsv_(&uplo,&nb_lign,&nb_SM,mat.v,resul.v,&ldb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{ Vecteur resul(b);
|
|
Mat_abstraite::Resolution_syst(b,resul,tole,maxi,rest);
|
|
return resul;
|
|
}
|
|
}
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case GAUSS:
|
|
{// ici les parties sup et inf de la largeur de bande sont supposée identique
|
|
// appel de la routine clapack
|
|
// /* Subroutine */ int dgbsv_(__CLPK_integer *n, __CLPK_integer *kl, __CLPK_integer *ku, __CLPK_integer *
|
|
// nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_integer *ipiv, __CLPK_doublereal *b,
|
|
// __CLPK_integer *ldb, __CLPK_integer *info);
|
|
dgbsv_(&nb_lign,&kl,&ku,&nb_SM,mat.v,&lda,ipiv,resul.v,&ldb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{ Mat_abstraite::Resolution_syst(b,resul,tole,maxi,rest);
|
|
}
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case CHOLESKY:
|
|
{// ici seule la partie sup de la largeur de bande est définie
|
|
// appel de la routine clapack
|
|
///* Subroutine */ int dpbsv_(char *uplo, __CLPK_integer *n, __CLPK_integer *kd, __CLPK_integer *
|
|
// nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_doublereal *b, __CLPK_integer *ldb,
|
|
// __CLPK_integer *info) __OSX_AVAILABLE_STARTING(__MAC_10_2,__IPHONE_4_0);
|
|
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
dpbsv_(&uplo,&nb_lign,&ku,&nb_SM,mat.v,&lda,resul.v,&ldb,&info);
|
|
GestionErreurDpbsv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{ Mat_abstraite::Resolution_syst(b,resul,tole,maxi,rest);
|
|
}
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n **** erreur: mauvais choix de matrice: " << Nom_matrice(type_matrice)
|
|
<< "\n MatLapack::Resol_syst(...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
// retour
|
|
return resul;
|
|
};
|
|
|
|
//2) avec en sortie le vecteur d'entree
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
Vecteur& MatLapack::Resol_systID
|
|
(Vecteur& b,const double &tole, const int maxi, const int rest)
|
|
{
|
|
ENTIER_POUR_CLAPACK nb_SM = 1; // un seul second membre
|
|
if (ipiv == NULL) {ipiv = new ENTIER_POUR_CLAPACK [nb_lign];};
|
|
ENTIER_POUR_CLAPACK info = 0; // info de sortie en fait c'est le retour de la fonction
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case GAUSS:
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
Verif_nb_composante(b,"Resol_systID");
|
|
#endif
|
|
// appel de la routine clapack
|
|
dgesv_(&nb_lign,&nb_SM,mat.v,&lda,ipiv,b.v,&ldb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
// le résultat est dans b.v
|
|
break;
|
|
}
|
|
default:
|
|
{Vecteur bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tole,maxi,rest);
|
|
b = bb;
|
|
};
|
|
};
|
|
break;
|
|
}
|
|
case CARREE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case CHOLESKY:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dppsv_(char *uplo, __CLPK_integer *n, __CLPK_integer *nrhs, __CLPK_doublereal
|
|
//*ap, __CLPK_doublereal *b, __CLPK_integer *ldb, __CLPK_integer *info) __OSX_AVAILABLE_STARTING(__MAC_10_2,__IPHONE_4_0);
|
|
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
dppsv_(&uplo,&nb_lign,&nb_SM,mat.v,b.v,&ldb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{Vecteur bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tole,maxi,rest);
|
|
b = bb;
|
|
};
|
|
}
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case GAUSS:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dgbsv_(__CLPK_integer *n, __CLPK_integer *kl, __CLPK_integer *ku, __CLPK_integer *
|
|
// nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_integer *ipiv, __CLPK_doublereal *b,
|
|
// __CLPK_integer *ldb, __CLPK_integer *info);
|
|
dgbsv_(&nb_lign,&kl,&ku,&nb_SM,mat.v,&lda,ipiv,b.v,&ldb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
// le résultat est dans b.v
|
|
break;
|
|
}
|
|
default:
|
|
{Vecteur bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tole,maxi,rest);
|
|
b = bb;
|
|
};
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case CHOLESKY:
|
|
{// ici seule la partie sup de la largeur de bande est définie
|
|
// appel de la routine clapack
|
|
///* Subroutine */ int dpbsv_(char *uplo, __CLPK_integer *n, __CLPK_integer *kd, __CLPK_integer *
|
|
// nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_doublereal *b, __CLPK_integer *ldb,
|
|
// __CLPK_integer *info) __OSX_AVAILABLE_STARTING(__MAC_10_2,__IPHONE_4_0);
|
|
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
dpbsv_(&uplo,&nb_lign,&ku,&nb_SM,mat.v,&lda,b.v,&ldb,&info);
|
|
GestionErreurDpbsv(info); // affichage d'erreurs éventuelles
|
|
// le résultat est dans b.v
|
|
break;
|
|
}
|
|
default:
|
|
{Vecteur bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tole,maxi,rest);
|
|
b = bb;
|
|
};
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n **** erreur: mauvais choix de matrice: " << Nom_matrice(type_matrice)
|
|
<< "\n MatLapack::Resol_syst(...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
// retour
|
|
return b;
|
|
};
|
|
|
|
//3) avec en entrée un tableau de vecteur second membre et
|
|
// en sortie un nouveau tableau de vecteurs
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
Tableau <Vecteur> MatLapack::Resol_syst
|
|
(const Tableau <Vecteur>& b,const double &tole, const int maxi, const int rest)
|
|
{
|
|
ENTIER_POUR_CLAPACK nb_SM = b.Taille(); // a priori plusieurs second membre
|
|
if (ipiv == NULL) {ipiv = new ENTIER_POUR_CLAPACK [nb_lign];};
|
|
ENTIER_POUR_CLAPACK info=0; // info de sortie
|
|
// on s'occupe de la matrice rectangulaire qui contient tous les seconds membres: init général
|
|
// ici le type de résolution importe peu, c'est du stockage
|
|
if (ptB == NULL)
|
|
{ptB = new MatLapack(nb_lign,nb_SM,false,0.,GAUSS,RIEN_PRECONDITIONNEMENT);}
|
|
else // on regarde si la taille est correcte
|
|
{ if ((nb_lign != ptB->Nb_ligne()) || (nb_SM != ptB->Nb_colonne()))
|
|
// c'est différent donc on redimentionne
|
|
{ ptB->Change_taille(nb_lign,nb_SM); };
|
|
};
|
|
// init particulier
|
|
for (int k=1;k<=nb_SM;k++)
|
|
{const Vecteur & bk=b(k);
|
|
for (int i=1;i<=nb_lign;i++)
|
|
(ptB->*Coor)(i,k)=bk(i);
|
|
};
|
|
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case GAUSS:
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
Verif_nb_composante(b,"Resol_syst");
|
|
#endif
|
|
// appel de la routine clapack
|
|
dgesv_(&nb_lign,&nb_SM,mat.v,&lda,ipiv,(*ptB).mat.v,&ldb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{Tableau <Vecteur> bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tole,maxi,rest);
|
|
return bb;
|
|
}
|
|
};
|
|
break;
|
|
}
|
|
case CARREE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case CHOLESKY:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dppsv_(char *uplo, __CLPK_integer *n, __CLPK_integer *nrhs,
|
|
// __CLPK_doublereal *ap, __CLPK_doublereal *b, __CLPK_integer *ldb, __CLPK_integer *info)
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
dppsv_(&uplo,&nb_lign,&nb_SM,mat.v,(*ptB).mat.v,&ldb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{Tableau <Vecteur> bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tole,maxi,rest);
|
|
return bb;
|
|
}
|
|
}
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case GAUSS:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dgbsv_(__CLPK_integer *n, __CLPK_integer *kl, __CLPK_integer *ku, __CLPK_integer *
|
|
// nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_integer *ipiv, __CLPK_doublereal *b,
|
|
// __CLPK_integer *ldb, __CLPK_integer *info);
|
|
dgbsv_(&nb_lign,&kl,&ku,&nb_SM,mat.v,&lda,ipiv,(*ptB).mat.v,&ldb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
// le résultat est dans b.v
|
|
break;
|
|
}
|
|
default:
|
|
{Tableau <Vecteur> bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tole,maxi,rest);
|
|
return bb;
|
|
}
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case CHOLESKY:
|
|
{// ici seule la partie sup de la largeur de bande est définie
|
|
// appel de la routine clapack
|
|
///* Subroutine */ int dpbsv_(char *uplo, __CLPK_integer *n, __CLPK_integer *kd, __CLPK_integer *
|
|
// nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_doublereal *b, __CLPK_integer *ldb,
|
|
// __CLPK_integer *info) __OSX_AVAILABLE_STARTING(__MAC_10_2,__IPHONE_4_0);
|
|
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
dpbsv_(&uplo,&nb_lign,&ku,&nb_SM,mat.v,&lda,(*ptB).mat.v,&ldb,&info);
|
|
GestionErreurDpbsv(info); // affichage d'erreurs éventuelles
|
|
// le résultat est dans b.v
|
|
break;
|
|
}
|
|
default:
|
|
{Tableau <Vecteur> bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tole,maxi,rest);
|
|
return bb;
|
|
}
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n **** erreur: mauvais choix de matrice: " << Nom_matrice(type_matrice)
|
|
<< "\n MatLapack::Resol_syst(...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
|
|
// récup du résultat
|
|
Tableau <Vecteur> bksortie(nb_SM);
|
|
for (int k=1;k<=nb_SM;k++)
|
|
{ bksortie(k).Change_taille(nb_lign);
|
|
Vecteur & bksortie_k=bksortie(k);
|
|
for (int i=1;i<=nb_lign;i++)
|
|
bksortie_k(i) = (ptB->*Coor_const)(i,k);
|
|
};
|
|
// retour des résultats
|
|
return bksortie;
|
|
};
|
|
|
|
|
|
//4) avec en entrée un tableau de vecteur second membre et
|
|
// en sortie le tableau de vecteurs d'entree
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
Tableau <Vecteur>& MatLapack::Resol_systID
|
|
(Tableau <Vecteur>& b,const double &tole, const int maxi, const int rest)
|
|
{
|
|
ENTIER_POUR_CLAPACK nb_SM = b.Taille(); // a priori plusieurs second membre
|
|
// init général
|
|
if (ipiv == NULL) {ipiv = new ENTIER_POUR_CLAPACK [nb_lign];};
|
|
// pour les seconds membres, le type de résolution importe peu, on met une grandeur par défaut
|
|
if (ptB == NULL)
|
|
{ ptB = new MatLapack(nb_lign,nb_SM,false,0.,GAUSS,RIEN_PRECONDITIONNEMENT);}
|
|
else // on regarde si la taille est correcte
|
|
{ if ((nb_lign != ptB->Nb_ligne()) || (nb_SM != ptB->Nb_colonne()))
|
|
// c'est différent donc on redimentionne
|
|
{ ptB->Change_taille(nb_lign,nb_SM); };
|
|
};
|
|
ENTIER_POUR_CLAPACK info=0; // info de sortie
|
|
// init particulier
|
|
for (int k=1;k<=nb_SM;k++)
|
|
{const Vecteur & bk=b(k);
|
|
for (int i=1;i<=nb_lign;i++)
|
|
(ptB->*Coor)(i,k)=bk(i);
|
|
};
|
|
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case GAUSS:
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
Verif_nb_composante(b,"Resol_systID");
|
|
#endif
|
|
// appel de la routine clapack
|
|
dgesv_(&nb_lign,&nb_SM,mat.v,&lda,ipiv,(*ptB).mat.v,&ldb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{Tableau <Vecteur> bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tole,maxi,rest);
|
|
b = bb;
|
|
return b;
|
|
};
|
|
};
|
|
break;
|
|
}
|
|
case CARREE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case CHOLESKY:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dppsv_(char *uplo, __CLPK_integer *n, __CLPK_integer *nrhs,
|
|
// __CLPK_doublereal *ap, __CLPK_doublereal *b, __CLPK_integer *ldb, __CLPK_integer *info)
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
dppsv_(&uplo,&nb_lign,&nb_SM,mat.v,(*ptB).mat.v,&ldb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{Tableau <Vecteur> bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tole,maxi,rest);
|
|
b = bb;
|
|
return b;
|
|
}
|
|
}
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case GAUSS:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dgbsv_(__CLPK_integer *n, __CLPK_integer *kl, __CLPK_integer *ku, __CLPK_integer *
|
|
// nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_integer *ipiv, __CLPK_doublereal *b,
|
|
// __CLPK_integer *ldb, __CLPK_integer *info);
|
|
dgbsv_(&nb_lign,&kl,&ku,&nb_SM,mat.v,&lda,ipiv,(*ptB).mat.v,&ldb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
// le résultat est dans b.v
|
|
break;
|
|
}
|
|
default:
|
|
{Tableau <Vecteur> bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tole,maxi,rest);
|
|
b = bb;
|
|
return b;
|
|
};
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case CHOLESKY:
|
|
{// ici seule la partie sup de la largeur de bande est définie
|
|
// appel de la routine clapack
|
|
///* Subroutine */ int dpbsv_(char *uplo, __CLPK_integer *n, __CLPK_integer *kd, __CLPK_integer *
|
|
// nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_doublereal *b, __CLPK_integer *ldb,
|
|
// __CLPK_integer *info) __OSX_AVAILABLE_STARTING(__MAC_10_2,__IPHONE_4_0);
|
|
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
dpbsv_(&uplo,&nb_lign,&ku,&nb_SM,mat.v,&lda,(*ptB).mat.v,&ldb,&info);
|
|
GestionErreurDpbsv(info); // affichage d'erreurs éventuelles
|
|
// le résultat est dans b.v
|
|
break;
|
|
}
|
|
default:
|
|
{Tableau <Vecteur> bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tole,maxi,rest);
|
|
b = bb;
|
|
return b;
|
|
}
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n **** erreur: mauvais choix de matrice: " << Nom_matrice(type_matrice)
|
|
<< "\n MatLapack::Resol_syst(...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
|
|
// récup du résultat
|
|
for (int k=1;k<=nb_SM;k++)
|
|
{Vecteur & bk=b(k);
|
|
for (int i=1;i<=nb_lign;i++)
|
|
bk(i) = (ptB->*Coor_const)(i,k);
|
|
};
|
|
return b;
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Realise le produit du vecteur vec par la matrice
|
|
Vecteur MatLapack::Prod_vec_mat (const Vecteur& vec) const
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (vec.Taille() != nb_lign)
|
|
{cout << "erreur de taille, la dimension du vecteur: " << vec.Taille()
|
|
<< " devrait etre identique au nombre de ligne de la matrice= " << nb_lign << " "
|
|
<< "\n MatLapack::Prod_vec_mat (const Vecteur& vec)" << endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:case CARREE_SYMETRIQUE_LAPACK:
|
|
{ Vecteur result(nb_col);
|
|
for (int j=1;j<= nb_col;j++)
|
|
{double* pt=&result(j);
|
|
for (int k=1; k<= nb_lign; k++)
|
|
(*pt) += vec(k) * (this->*Coor_const)(k,j);
|
|
};
|
|
return result;
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:case BANDE_SYMETRIQUE_LAPACK:
|
|
{ // on considère que la matrice doit représenter une matrice carrée
|
|
int lb=ku; int dim = nb_lign;
|
|
Vecteur res(dim);
|
|
for (int j=1;j<=dim;j++)
|
|
for (int i= MaX(1,j-lb);i<= MiN(dim,j+lb);i++)
|
|
res(j) += vec(i) * (this->*Coor_const)(i,j);
|
|
return res;
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n **** erreur: mauvais choix de matrice: " << Nom_matrice(type_matrice)
|
|
<< "\n MatLapack::Prod_vec_mat(...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
Vecteur ret; // pour taire le compilo
|
|
return ret; // "
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
Vecteur& MatLapack::Prod_vec_mat (const Vecteur& vec,Vecteur& res) const
|
|
// Realise le produit du vecteur vec par la matrice
|
|
// ici on se sert du second vecteur pour le résultat
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (vec.Taille() != nb_lign)
|
|
{cout << "erreur de taille, la dimension du vecteur: " << vec.Taille()
|
|
<< " devrait etre identique au nombre de ligne de la matrice= " << nb_lign << " "
|
|
<< "\n MatLapack::Prod_vec_mat(const Vecteur& vec,Vecteur& result)" << endl;
|
|
Sortie(1);
|
|
}
|
|
#endif
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:case CARREE_SYMETRIQUE_LAPACK:
|
|
{ res.Change_taille(nb_col);
|
|
for (int j=1;j<= nb_col;j++)
|
|
{ double* pt=&res(j);
|
|
for (int k=1; k<= nb_lign; k++)
|
|
(*pt) += vec(k) * (this->*Coor_const)(k,j);
|
|
};
|
|
return res;
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:case BANDE_SYMETRIQUE_LAPACK:
|
|
{ // on considère que la matrice doit représenter une matrice carrée
|
|
int lb=ku; int dim = nb_lign;
|
|
res.Change_taille(dim);
|
|
for (int j=1;j<=dim;j++)
|
|
{int maxj = MiN(dim,j+lb)+1; // +1 pour < strict dans la boucle
|
|
for (int i= MaX(1,j-lb);i< maxj;i++)
|
|
res(j) += vec(i) * (this->*Coor_const)(i,j);
|
|
};
|
|
return res;
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n **** erreur: mauvais choix de matrice: " << Nom_matrice(type_matrice)
|
|
<< "\n MatLapack::Prod_vec_mat(...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
Vecteur ret; // pour taire le compilo
|
|
return ret; // "
|
|
};
|
|
|
|
//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
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
Vecteur& MatLapack::Resol_systID_2 (const Vecteur& b,Vecteur& vortie
|
|
, const double &tole,const int maxi,const int rest)
|
|
{// init général
|
|
ENTIER_POUR_CLAPACK nb_SM = 1; // un seul second membre
|
|
if (ipiv == NULL) {ipiv = new ENTIER_POUR_CLAPACK [nb_lign];};
|
|
ENTIER_POUR_CLAPACK info = 0; // info de sortie en fait c'est le retour de la fonction
|
|
vortie=b; // Vecteur solution initialisé = second membre
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{switch (type_resolution)
|
|
{case GAUSS:
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
Verif_nb_composante(b,"Resol_systID_2");
|
|
#endif
|
|
// appel de la routine clapack
|
|
dgesv_(&nb_lign,&nb_SM,mat.v,&lda,ipiv,vortie.v,&ldb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{ Mat_abstraite::Resolution_syst(b,vortie,tole,maxi,rest);};
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{switch (type_resolution)
|
|
{case GAUSS:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dgbsv_(__CLPK_integer *n, __CLPK_integer *kl, __CLPK_integer *ku, __CLPK_integer *
|
|
// nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_integer *ipiv, __CLPK_doublereal *b,
|
|
// __CLPK_integer *ldb, __CLPK_integer *info);
|
|
dgbsv_(&nb_lign,&kl,&ku,&nb_SM,mat.v,&lda,ipiv,vortie.v,&ldb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{ Mat_abstraite::Resolution_syst(b,vortie,tole,maxi,rest);};
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case CHOLESKY:
|
|
{// ici seule la partie sup de la largeur de bande est définie
|
|
// appel de la routine clapack
|
|
///* Subroutine */ int dpbsv_(char *uplo, __CLPK_integer *n, __CLPK_integer *kd, __CLPK_integer *
|
|
// nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_doublereal *b, __CLPK_integer *ldb,
|
|
// __CLPK_integer *info) __OSX_AVAILABLE_STARTING(__MAC_10_2,__IPHONE_4_0);
|
|
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
dpbsv_(&uplo,&nb_lign,&ku,&nb_SM,mat.v,&lda,vortie.v,&ldb,&info);
|
|
GestionErreurDpbsv(info); // affichage d'erreurs éventuelles
|
|
// le résultat est dans vortie.v
|
|
break;
|
|
}
|
|
default:
|
|
{ Mat_abstraite::Resolution_syst(b,vortie,tole,maxi,rest);};
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n **** erreur: mauvais choix de matrice: " << Nom_matrice(type_matrice)
|
|
<< "\n MatLapack::Resol_syst(...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
// retour
|
|
return vortie;
|
|
};
|
|
|
|
// ===== 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
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
void MatLapack::Preparation_resol()
|
|
{
|
|
// init général
|
|
if (ipiv == NULL) {ipiv = new ENTIER_POUR_CLAPACK [nb_lign];}; // normalement ipiv doit avoir ici
|
|
ENTIER_POUR_CLAPACK info = 0; // info de sortie en fait c'est le retour de la fonction
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{ case GAUSS:
|
|
{ // appel de la routine lapack
|
|
dgetrf_(&nb_lign,&nb_col,mat.v,&lda,ipiv,&info);
|
|
break;
|
|
}
|
|
// dans les autres cas on ne fait rien
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::Preparation_resol (.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
break;
|
|
}
|
|
case CARREE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case CHOLESKY:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dpptrf_(char *uplo, __CLPK_integer *n, __CLPK_doublereal *ap, __CLPK_integer *info)
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
dpptrf_(&uplo,&nb_lign,mat.v,&info) ;
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
// dans les autres cas on ne fait rien
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::Preparation_resol (.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{ case GAUSS:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dgbtrf_(__CLPK_integer *m, __CLPK_integer *n, __CLPK_integer *kl, __CLPK_integer *ku,
|
|
// __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_integer *ipiv, __CLPK_integer *info);
|
|
dgbtrf_(&nb_lign,&nb_lign,&kl,&ku,mat.v,&lda,ipiv,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::Preparation_resol (.. "<<endl;
|
|
Sortie(1);
|
|
// dans les autres cas on ne fait rien
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{ case CHOLESKY:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dpbtrf_(char *uplo, __CLPK_integer *n, __CLPK_integer *kd, __CLPK_doublereal *
|
|
// ab, __CLPK_integer *ldab, __CLPK_integer *info) __OSX_AVAILABLE_STARTING(__MAC_10_2,__IPHONE_4_0);
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
dpbtrf_(&uplo,&nb_lign,&ku,mat.v,&lda,&info);
|
|
GestionErreurDpbsv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::Preparation_resol (.. "<<endl;
|
|
Sortie(1);
|
|
// dans les autres cas on ne fait rien
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n **** erreur: mauvais choix de matrice: " << Nom_matrice(type_matrice)
|
|
<< "\n MatLapack::Preparation_resol(...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
// retour
|
|
};
|
|
|
|
// 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
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
Vecteur MatLapack::Simple_Resol_syst (const Vecteur& b,const double &tol
|
|
,const int maxi,const int rest) const
|
|
{// init général
|
|
ENTIER_POUR_CLAPACK nb_SM = 1; // un seul second membre
|
|
if ((type_matrice != BANDE_SYMETRIQUE_LAPACK) && (ipiv == NULL))
|
|
{cout << "\n erreur: la methode MatLapack::Preparation_resol() n'a pas ete appelee "
|
|
<< " le tableau ipiv n'est pas correcte, on ne peut pas continuer ! "
|
|
<< "\n MatLapack::Simple_Resol_syst (.... ";
|
|
Sortie(1);
|
|
};
|
|
ENTIER_POUR_CLAPACK info = 0; // info de sortie en fait c'est le retour de la fonction
|
|
Vecteur resul(b); // Vecteur solution initialisé = second membre
|
|
char trans = 'N'; // "No transpose"
|
|
// comme la méthode est const, on doit utiliser des cast, se qui n'est pas sûr mais difficile de faire autrement !
|
|
double * matric_v = (double *) mat.v; // on impose un cast
|
|
ENTIER_POUR_CLAPACK * ipiva = (ENTIER_POUR_CLAPACK*) ipiv; // ""
|
|
// toujour à cause du const, on est obligé d'utiliser des variables intermédiaires
|
|
ENTIER_POUR_CLAPACK nb_ligne=nb_lign;
|
|
ENTIER_POUR_CLAPACK ldaa = lda;
|
|
ENTIER_POUR_CLAPACK ldbb = ldb;
|
|
ENTIER_POUR_CLAPACK kll = kl ;
|
|
ENTIER_POUR_CLAPACK kuu = ku ;
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{ case GAUSS:
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
Verif_nb_composante(b,"Simple_Resol_syst");
|
|
#endif
|
|
// appel de la routine clapack
|
|
// /* Subroutine */ int dgetrs_(char *trans, __CLPK_integer *n, __CLPK_integer *nrhs,
|
|
// __CLPK_doublereal *a, __CLPK_integer *lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer *
|
|
// ldb, __CLPK_integer *info);
|
|
dgetrs_(&trans,&nb_ligne,&nb_SM,matric_v,&ldaa,ipiva,resul.v,&ldbb,&info);
|
|
break;
|
|
}
|
|
default:
|
|
{ Mat_abstraite::Resolution_syst(b,resul,tol,maxi,rest);}
|
|
};
|
|
break;
|
|
}
|
|
case CARREE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case CHOLESKY:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dpptrs_(char *uplo, __CLPK_integer *n, __CLPK_integer *nrhs,
|
|
// __CLPK_doublereal *ap, __CLPK_doublereal *b, __CLPK_integer *ldb, __CLPK_integer *info)
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
dpptrs_(&uplo,&nb_ligne,&nb_SM,matric_v,resul.v,&ldbb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{ Mat_abstraite::Resolution_syst(b,resul,tol,maxi,rest);}
|
|
}
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{ case GAUSS:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dgbtrs_(char *trans, __CLPK_integer *n, __CLPK_integer *kl, __CLPK_integer *
|
|
// ku, __CLPK_integer *nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_integer *ipiv,
|
|
// __CLPK_doublereal *b, __CLPK_integer *ldb, __CLPK_integer *info);
|
|
dgbtrs_(&trans,&nb_ligne,&kll,&kuu,&nb_SM,matric_v,&ldaa,ipiva,resul.v,&ldbb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{ Mat_abstraite::Resolution_syst(b,resul,tol,maxi,rest);}
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{ case CHOLESKY:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dpbtrs_(char *uplo, __CLPK_integer *n, __CLPK_integer *kd, __CLPK_integer *
|
|
// nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_doublereal *b, __CLPK_integer *ldb,
|
|
// __CLPK_integer *info) __OSX_AVAILABLE_STARTING(__MAC_10_2,__IPHONE_4_0);
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
dpbtrs_(&uplo,&nb_ligne,&kuu,&nb_SM,matric_v,&ldaa,resul.v,&ldbb,&info);
|
|
GestionErreurDpbsv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{ Mat_abstraite::Resolution_syst(b,resul,tol,maxi,rest);}
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n **** erreur: mauvais choix de matrice: " << Nom_matrice(type_matrice)
|
|
<< "\n MatLapack::Simple_Resol_syst (...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
// retour
|
|
return resul;
|
|
};
|
|
|
|
// b) avec en sortie le vecteur d'entree
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
Vecteur& MatLapack::Simple_Resol_systID (Vecteur& b,const double &tol
|
|
,const int maxi,const int restart ) const
|
|
{// init général
|
|
ENTIER_POUR_CLAPACK nb_SM = 1; // un seul second membre
|
|
if ((type_matrice != BANDE_SYMETRIQUE_LAPACK) && (ipiv == NULL))
|
|
{cout << "\n erreur: la methode MatLapack::Preparation_resol() n'a pas ete appelee "
|
|
<< " le tableau ipiv n'est pas correcte, on ne peut pas continuer ! "
|
|
<< "\n MatLapack::Simple_Resol_systID (.... ";
|
|
Sortie(1);
|
|
};
|
|
ENTIER_POUR_CLAPACK ldb = MaX(1,nb_lign);
|
|
ENTIER_POUR_CLAPACK info = 0; // info de sortie en fait c'est le retour de la fonction
|
|
char trans = 'N'; // "No transpose"
|
|
// comme la méthode est const, on doit utiliser des cast, se qui n'est pas sûr mais difficile de faire autrement !
|
|
double * matric_v = (double *) mat.v; // on impose un cast
|
|
ENTIER_POUR_CLAPACK * ipiva = (ENTIER_POUR_CLAPACK*) ipiv; // ""
|
|
// toujour à cause du const, on est obligé d'utiliser des variables intermédiaires
|
|
ENTIER_POUR_CLAPACK nb_ligne=nb_lign;
|
|
ENTIER_POUR_CLAPACK ldaa = lda;
|
|
ENTIER_POUR_CLAPACK ldbb = ldb;
|
|
ENTIER_POUR_CLAPACK kll = kl ;
|
|
ENTIER_POUR_CLAPACK kuu = ku ;
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{ case GAUSS:
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
Verif_nb_composante(b,"Simple_Resol_systID");
|
|
#endif
|
|
// appel de la routine clapack
|
|
// /* Subroutine */ int dgetrs_(char *trans, __CLPK_integer *n, __CLPK_integer *nrhs,
|
|
// __CLPK_doublereal *a, __CLPK_integer *lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer *
|
|
// ldb, __CLPK_integer *info);
|
|
dgetrs_(&trans,&nb_ligne,&nb_SM,matric_v,&ldaa,ipiva,b.v,&ldbb,&info);
|
|
break;
|
|
}
|
|
default:
|
|
{ Vecteur bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tol,maxi,restart);
|
|
b = bb;
|
|
};
|
|
};
|
|
break;
|
|
}
|
|
case CARREE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case CHOLESKY:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dpptrs_(char *uplo, __CLPK_integer *n, __CLPK_integer *nrhs,
|
|
// __CLPK_doublereal *ap, __CLPK_doublereal *b, __CLPK_integer *ldb, __CLPK_integer *info)
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
dpptrs_(&uplo,&nb_ligne,&nb_SM,matric_v,b.v,&ldbb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{ Vecteur bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tol,maxi,restart);
|
|
b = bb;
|
|
};
|
|
}
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{ case GAUSS:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dgbtrs_(char *trans, __CLPK_integer *n, __CLPK_integer *kl, __CLPK_integer *
|
|
// ku, __CLPK_integer *nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_integer *ipiv,
|
|
// __CLPK_doublereal *b, __CLPK_integer *ldb, __CLPK_integer *info);
|
|
dgbtrs_(&trans,&nb_ligne,&kll,&kuu,&nb_SM,matric_v,&ldaa,ipiva,b.v,&ldbb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{ Vecteur bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tol,maxi,restart);
|
|
b = bb;
|
|
};
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{ case CHOLESKY:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dpbtrs_(char *uplo, __CLPK_integer *n, __CLPK_integer *kd, __CLPK_integer *
|
|
// nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_doublereal *b, __CLPK_integer *ldb,
|
|
// __CLPK_integer *info) __OSX_AVAILABLE_STARTING(__MAC_10_2,__IPHONE_4_0);
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
dpbtrs_(&uplo,&nb_ligne,&kuu,&nb_SM,matric_v,&ldaa,b.v,&ldbb,&info);
|
|
GestionErreurDpbsv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{ Vecteur bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tol,maxi,restart);
|
|
b = bb;
|
|
};
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n **** erreur: mauvais choix de matrice: " << Nom_matrice(type_matrice)
|
|
<< "\n MatLapack::Simple_Resol_syst (...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
return b;
|
|
};
|
|
|
|
// c) avec en entrée un tableau de vecteur second membre et
|
|
// en sortie un nouveau tableau de vecteurs
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
Tableau <Vecteur> MatLapack::Simple_Resol_syst
|
|
(const Tableau <Vecteur>& b,const double &tol
|
|
,const int maxi,const int restart ) const
|
|
{
|
|
ENTIER_POUR_CLAPACK nb_SM = 1; // un seul second membre à la fois
|
|
if ((type_matrice != BANDE_SYMETRIQUE_LAPACK) && (ipiv == NULL))
|
|
{ cout << "\n erreur: la methode MatLapack::Preparation_resol() n'a pas ete appeled "
|
|
<< " le tableau ipiv n'est pas correcte, on ne peut pas continuer ! "
|
|
<< "\n MatLapack::Simple_Resol_syst (Tableau <Vecteur>&.... ";
|
|
Sortie(1);
|
|
};
|
|
// #ifdef MISE_AU_POINT
|
|
// Verif_nb_composante(b,"Simple_Resol_syst");
|
|
// #endif
|
|
ENTIER_POUR_CLAPACK ldb = MaX(1,nb_lign);
|
|
ENTIER_POUR_CLAPACK info = 0; // info de sortie en fait c'est le retour de la fonction
|
|
char trans = 'N'; // "No transpose"
|
|
// comme la méthode est const, on doit utiliser des cast, ce qui n'est pas sûr mais difficile de faire autrement !
|
|
double * matric_v = (double *) mat.v; // on impose un cast
|
|
ENTIER_POUR_CLAPACK * ipiva = (ENTIER_POUR_CLAPACK*) ipiv; // ""
|
|
// toujour à cause du const, on est obligé d'utiliser des variables intermédiaires
|
|
ENTIER_POUR_CLAPACK nb_ligne=nb_lign;
|
|
ENTIER_POUR_CLAPACK ldaa = lda;
|
|
ENTIER_POUR_CLAPACK ldbb = ldb;
|
|
ENTIER_POUR_CLAPACK kll = kl ;
|
|
ENTIER_POUR_CLAPACK kuu = ku ;
|
|
// appel de la routine clapack
|
|
int nb_vrai_sm=b.Taille();
|
|
Tableau <Vecteur> bretour(b);
|
|
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{ case GAUSS:
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
Verif_nb_composante(b,"Simple_Resol_systID");
|
|
#endif
|
|
// appel de la routine clapack
|
|
// /* Subroutine */ int dgetrs_(char *trans, __CLPK_integer *n, __CLPK_integer *nrhs,
|
|
// __CLPK_doublereal *a, __CLPK_integer *lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer *
|
|
// ldb, __CLPK_integer *info);
|
|
for (int im=1;im<=nb_vrai_sm;im++)
|
|
dgetrs_(&trans,&nb_ligne,&nb_SM,matric_v,&ldaa,ipiva,bretour(im).v,&ldbb,&info);
|
|
break;
|
|
}
|
|
default:
|
|
{ Mat_abstraite::Resolution_syst(b,bretour,tol,maxi,restart);};
|
|
};
|
|
break;
|
|
}
|
|
case CARREE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case CHOLESKY:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dpptrs_(char *uplo, __CLPK_integer *n, __CLPK_integer *nrhs,
|
|
// __CLPK_doublereal *ap, __CLPK_doublereal *b, __CLPK_integer *ldb, __CLPK_integer *info)
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
for (int im=1;im<=nb_vrai_sm;im++)
|
|
{dpptrs_(&uplo,&nb_ligne,&nb_SM,matric_v,bretour(im).v,&ldbb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
{ Mat_abstraite::Resolution_syst(b,bretour,tol,maxi,restart);};
|
|
}
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{ case GAUSS:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dgbtrs_(char *trans, __CLPK_integer *n, __CLPK_integer *kl, __CLPK_integer *
|
|
// ku, __CLPK_integer *nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_integer *ipiv,
|
|
// __CLPK_doublereal *b, __CLPK_integer *ldb, __CLPK_integer *info);
|
|
for (int im=1;im<=nb_vrai_sm;im++)
|
|
dgbtrs_(&trans,&nb_ligne,&kll,&kuu,&nb_SM,matric_v,&ldaa,ipiva,bretour(im).v,&ldbb,&info);
|
|
break;
|
|
}
|
|
default:
|
|
{ Mat_abstraite::Resolution_syst(b,bretour,tol,maxi,restart);};
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{ case CHOLESKY:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dpbtrs_(char *uplo, __CLPK_integer *n, __CLPK_integer *kd, __CLPK_integer *
|
|
// nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_doublereal *b, __CLPK_integer *ldb,
|
|
// __CLPK_integer *info) __OSX_AVAILABLE_STARTING(__MAC_10_2,__IPHONE_4_0);
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
for (int im=1;im<=nb_vrai_sm;im++)
|
|
dpbtrs_(&uplo,&nb_ligne,&kuu,&nb_SM,matric_v,&ldaa,bretour(im).v,&ldbb,&info);
|
|
GestionErreurDpbsv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{ Mat_abstraite::Resolution_syst(b,bretour,tol,maxi,restart);};
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n **** erreur: mauvais choix de matrice: " << Nom_matrice(type_matrice)
|
|
<< "\n MatLapack::Simple_Resol_syst (...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
return bretour;
|
|
};
|
|
// d) avec en entrée un tableau de vecteur second membre et
|
|
// en sortie le tableau de vecteurs d'entree
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
Tableau <Vecteur>& MatLapack::Simple_Resol_systID
|
|
(Tableau <Vecteur>& b,const double &tol
|
|
,const int maxi,const int restart ) const
|
|
{
|
|
ENTIER_POUR_CLAPACK nb_SM = 1; // un seul second membre à la fois
|
|
if ((type_matrice != BANDE_SYMETRIQUE_LAPACK) && (ipiv == NULL))
|
|
{ cout << "\n erreur: la methode MatLapack::Preparation_resol() n'a pas ete appelee "
|
|
<< " le tableau ipiv n'est pas correcte, on ne peut pas continuer ! "
|
|
<< "\n MatLapack::Simple_Resol_systID (Tableau <Vecteur>&.... ";
|
|
Sortie(1);
|
|
};
|
|
// #ifdef MISE_AU_POINT
|
|
// Verif_nb_composante(b,"Simple_Resol_systID");
|
|
// #endif
|
|
ENTIER_POUR_CLAPACK ldb = MaX(1,nb_lign);
|
|
ENTIER_POUR_CLAPACK info = 0; // info de sortie en fait c'est le retour de la fonction
|
|
char trans = 'N'; // "No transpose"
|
|
// comme la méthode est const, on doit utiliser des cast, ce qui n'est pas sûr mais difficile de faire autrement !
|
|
double * matric_v = (double *) mat.v; // on impose un cast
|
|
ENTIER_POUR_CLAPACK * ipiva = (ENTIER_POUR_CLAPACK*) ipiv; // ""
|
|
// toujour à cause du const, on est obligé d'utiliser des variables intermédiaires
|
|
ENTIER_POUR_CLAPACK nb_ligne=nb_lign;
|
|
ENTIER_POUR_CLAPACK ldaa = lda;
|
|
ENTIER_POUR_CLAPACK ldbb = ldb;
|
|
ENTIER_POUR_CLAPACK kll = kl ;
|
|
ENTIER_POUR_CLAPACK kuu = ku ;
|
|
// appel de la routine clapack
|
|
int nb_vrai_sm=b.Taille();
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{ case GAUSS:
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
Verif_nb_composante(b,"Simple_Resol_systID");
|
|
#endif
|
|
// appel de la routine clapack
|
|
// /* Subroutine */ int dgetrs_(char *trans, __CLPK_integer *n, __CLPK_integer *nrhs,
|
|
// __CLPK_doublereal *a, __CLPK_integer *lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer *
|
|
// ldb, __CLPK_integer *info);
|
|
for (int im=1;im<=nb_vrai_sm;im++)
|
|
dgetrs_(&trans,&nb_ligne,&nb_SM,matric_v,&ldaa,ipiva,b(im).v,&ldbb,&info);
|
|
break;
|
|
}
|
|
default:
|
|
{Tableau <Vecteur> bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tol,maxi,restart);
|
|
b = bb;
|
|
};
|
|
};
|
|
break;
|
|
}
|
|
case CARREE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case CHOLESKY:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dpptrs_(char *uplo, __CLPK_integer *n, __CLPK_integer *nrhs,
|
|
// __CLPK_doublereal *ap, __CLPK_doublereal *b, __CLPK_integer *ldb, __CLPK_integer *info)
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
for (int im=1;im<=nb_vrai_sm;im++)
|
|
{dpptrs_(&uplo,&nb_ligne,&nb_SM,matric_v,b(im).v,&ldbb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
{Tableau <Vecteur> bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tol,maxi,restart);
|
|
b = bb;
|
|
};
|
|
}
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{ case GAUSS:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dgbtrs_(char *trans, __CLPK_integer *n, __CLPK_integer *kl, __CLPK_integer *
|
|
// ku, __CLPK_integer *nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_integer *ipiv,
|
|
// __CLPK_doublereal *b, __CLPK_integer *ldb, __CLPK_integer *info);
|
|
for (int im=1;im<=nb_vrai_sm;im++)
|
|
dgbtrs_(&trans,&nb_ligne,&kll,&kuu,&nb_SM,matric_v,&ldaa,ipiva,b(im).v,&ldbb,&info);
|
|
break;
|
|
}
|
|
default:
|
|
{Tableau <Vecteur> bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tol,maxi,restart);
|
|
b = bb;
|
|
};
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{ case CHOLESKY:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dpbtrs_(char *uplo, __CLPK_integer *n, __CLPK_integer *kd, __CLPK_integer *
|
|
// nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_doublereal *b, __CLPK_integer *ldb,
|
|
// __CLPK_integer *info) __OSX_AVAILABLE_STARTING(__MAC_10_2,__IPHONE_4_0);
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
for (int im=1;im<=nb_vrai_sm;im++)
|
|
dpbtrs_(&uplo,&nb_ligne,&kuu,&nb_SM,matric_v,&ldaa,b(im).v,&ldbb,&info);
|
|
GestionErreurDpbsv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{Tableau <Vecteur> bb(b);
|
|
Mat_abstraite::Resolution_syst(b,bb,tol,maxi,restart);
|
|
b = bb;
|
|
};
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n **** erreur: mauvais choix de matrice: " << Nom_matrice(type_matrice)
|
|
<< "\n MatLapack::Simple_Resol_systID (...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
return b;
|
|
};
|
|
// 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
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
Vecteur& MatLapack::Simple_Resol_systID_2 (const Vecteur& b,Vecteur& vortie
|
|
, const double &tol,const int maxi,const int rest) const
|
|
{// init général
|
|
ENTIER_POUR_CLAPACK nb_SM = 1; // un seul second membre
|
|
if ((type_matrice != BANDE_SYMETRIQUE_LAPACK) && (ipiv == NULL))
|
|
{cout << "\n erreur: la methode MatLapack::Preparation_resol() n'a pas ete appelee "
|
|
<< " le tableau ipiv n'est pas correcte, on ne peut pas continuer ! "
|
|
<< "\n MatLapack::Simple_Resol_systID_2 (.... ";
|
|
Sortie(1);
|
|
};
|
|
ENTIER_POUR_CLAPACK ldb = MaX(1,nb_lign);
|
|
ENTIER_POUR_CLAPACK info = 0; // info de sortie en fait c'est le retour de la fonction
|
|
vortie=b; // init
|
|
char trans = 'N'; // "No transpose"
|
|
// comme la méthode est const, on doit utiliser des cast, se qui n'est pas sûr mais difficile de faire autrement !
|
|
double * matric_v = (double *) mat.v; // on impose un cast
|
|
ENTIER_POUR_CLAPACK * ipiva = (ENTIER_POUR_CLAPACK*) ipiv; // ""
|
|
// toujour à cause du const, on est obligé d'utiliser des variables intermédiaires
|
|
ENTIER_POUR_CLAPACK nb_ligne=nb_lign;
|
|
ENTIER_POUR_CLAPACK ldaa = lda;
|
|
ENTIER_POUR_CLAPACK ldbb = ldb;
|
|
ENTIER_POUR_CLAPACK kll = kl ;
|
|
ENTIER_POUR_CLAPACK kuu = ku ;
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{ case GAUSS:
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
Verif_nb_composante(b,"Simple_Resol_systID_2");
|
|
#endif
|
|
// appel de la routine clapack
|
|
// /* Subroutine */ int dgetrs_(char *trans, __CLPK_integer *n, __CLPK_integer *nrhs,
|
|
// __CLPK_doublereal *a, __CLPK_integer *lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer *
|
|
// ldb, __CLPK_integer *info);
|
|
dgetrs_(&trans,&nb_ligne,&nb_SM,matric_v,&ldaa,ipiva,vortie.v,&ldbb,&info);
|
|
break;
|
|
}
|
|
default:
|
|
{ Mat_abstraite::Resolution_syst(b,vortie,tol,maxi,rest);}
|
|
};
|
|
break;
|
|
}
|
|
case CARREE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{case CHOLESKY:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dpptrs_(char *uplo, __CLPK_integer *n, __CLPK_integer *nrhs,
|
|
// __CLPK_doublereal *ap, __CLPK_doublereal *b, __CLPK_integer *ldb, __CLPK_integer *info)
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
dpptrs_(&uplo,&nb_ligne,&nb_SM,matric_v,vortie.v,&ldbb,&info);
|
|
GestionErreurDg_sv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{ Mat_abstraite::Resolution_syst(b,vortie,tol,maxi,rest);}
|
|
}
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{ case GAUSS:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dgbtrs_(char *trans, __CLPK_integer *n, __CLPK_integer *kl, __CLPK_integer *
|
|
// ku, __CLPK_integer *nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_integer *ipiv,
|
|
// __CLPK_doublereal *b, __CLPK_integer *ldb, __CLPK_integer *info);
|
|
dgbtrs_(&trans,&nb_ligne,&kll,&kuu,&nb_SM,matric_v,&ldaa,ipiva,vortie.v,&ldbb,&info);
|
|
break;
|
|
}
|
|
default:
|
|
{ Mat_abstraite::Resolution_syst(b,vortie,tol,maxi,rest);}
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{ switch (type_resolution)
|
|
{ case CHOLESKY:
|
|
{// appel de la routine clapack
|
|
// /* Subroutine */ int dpbtrs_(char *uplo, __CLPK_integer *n, __CLPK_integer *kd, __CLPK_integer *
|
|
// nrhs, __CLPK_doublereal *ab, __CLPK_integer *ldab, __CLPK_doublereal *b, __CLPK_integer *ldb,
|
|
// __CLPK_integer *info) __OSX_AVAILABLE_STARTING(__MAC_10_2,__IPHONE_4_0);
|
|
char uplo = 'U'; // up diagonal: stockage au dessus de la diagonale
|
|
dpbtrs_(&uplo,&nb_ligne,&kuu,&nb_SM,matric_v,&ldaa,vortie.v,&ldbb,&info);
|
|
GestionErreurDpbsv(info); // affichage d'erreurs éventuelles
|
|
break;
|
|
}
|
|
default:
|
|
{ Mat_abstraite::Resolution_syst(b,vortie,tol,maxi,rest);}
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n **** erreur: mauvais choix de matrice: " << Nom_matrice(type_matrice)
|
|
<< "\n MatLapack::Simple_Resol_systID_2 (...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
// retour
|
|
return vortie;
|
|
};
|
|
// ===== FIN RÉSOLUTION EN DEUX TEMPS ================ :
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Realise le produit de la matrice par le vecteur vec
|
|
Vecteur MatLapack::Prod_mat_vec (const Vecteur& vec) const
|
|
{ Vecteur result(nb_lign);
|
|
#ifdef MISE_AU_POINT
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK: case CARREE_SYMETRIQUE_LAPACK:
|
|
{ if (nb_col != vec.Taille())
|
|
{ cout << "erreur de taille, le vecteur doit avoir une dimension identique au nombre"
|
|
<< " de colonne de la matrice";
|
|
cout << " \n Prod_mat_vec, tailleligneMatrice = " <<nb_col
|
|
<< " taillevecteur = " << vec.Taille();
|
|
cout << " Vecteur MatLapack::Prod_mat_vec (... " << endl;
|
|
Sortie (1);
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK: case BANDE_SYMETRIQUE_LAPACK:
|
|
{ if (vec.Taille() != nb_lign)
|
|
{ cout << "erreur de taille, le vecteur doit avoir une dimension identique au nombre"
|
|
<< " de colonne de la matrice";
|
|
cout << " \n Prod_mat_vec, tailleligneMatrice = " <<nb_lign
|
|
<< " taillevecteur = " << vec.Taille();
|
|
cout << " Vecteur MatLapack::Prod_mat_vec (... " << endl;
|
|
Sortie (1);
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n Vecteur MatLapack::Prod_mat_vec (... " << endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
// y := alpha*A*x + beta*y,
|
|
// avec ici beta = 0.,
|
|
double alpha=1.;double beta = 0.;
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ //void cblas_dgemv(const enum CBLAS_ORDER Order,
|
|
// const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
|
|
// const double alpha, const double *A, const int lda,
|
|
// const double *X, const int incX, const double beta,
|
|
// double *Y, const int incY) __OSX_AVAILABLE_STARTING(__MAC_10_2,__IPHONE_4_0);
|
|
cblas_dgemv(CblasRowMajor,CblasNoTrans, nb_lign, lda,alpha,mat.v,lda,vec.v,1,beta,result.v,1);
|
|
//for (int i=1;i<=nb_lign;i++)
|
|
// {double* pt= &result(i);
|
|
// for (int k=1;k<=nb_col;k++)
|
|
// (*pt) += (this->*Coor_const)(i,k) * vec(k);
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ //void cblas_dgbmv(const enum CBLAS_ORDER Order,
|
|
// const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
|
|
// const int KL, const int KU, const double alpha,
|
|
// const double *A, const int lda, const double *X,
|
|
// const int incX, const double beta, double *Y, const int incY)
|
|
|
|
cblas_dgbmv(CblasRowMajor,CblasNoTrans, nb_lign, lda,kl,ku,alpha,mat.v,lda,vec.v,1,beta,result.v,1);
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{ //void cblas_dsbmv(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
|
|
// const int N, const int K, const double alpha, const double *A,
|
|
// const int lda, const double *X, const int incX,
|
|
// const double beta, double *Y, const int incY) __OSX_AVAILABLE_STARTING(__MAC_10_2,
|
|
|
|
cblas_dsbmv(CblasRowMajor, CblasUpper,nb_lign,ku,alpha,mat.v,lda,vec.v,1,beta,result.v,1);
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n **** erreur: mauvais choix de matrice: " << Nom_matrice(type_matrice)
|
|
<< "\n MatLapack::Simple_Resol_systID_2 (...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
return result;
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Realise le produit de la matrice par le vecteur vec
|
|
// ici on se sert du second vecteur pour le résultat
|
|
Vecteur& MatLapack::Prod_mat_vec (const Vecteur& vec,Vecteur& result) const
|
|
{ result.Change_taille(nb_lign);
|
|
#ifdef MISE_AU_POINT
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK: case CARREE_SYMETRIQUE_LAPACK:
|
|
{ if (nb_col != vec.Taille())
|
|
{ cout << "erreur de taille, le vecteur doit avoir une dimension identique au nombre"
|
|
<< " de colonne de la matrice";
|
|
cout << " \n Prod_mat_vec, tailleligneMatrice = " <<nb_col
|
|
<< " taillevecteur = " << vec.Taille();
|
|
cout << " Vecteur& MatLapack::Prod_mat_vec (... " << endl;
|
|
Sortie (1);
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK: case BANDE_SYMETRIQUE_LAPACK:
|
|
{ if (vec.Taille() != nb_lign)
|
|
{ cout << "erreur de taille, le vecteur doit avoir une dimension identique au nombre"
|
|
<< " de colonne de la matrice";
|
|
cout << " \n Prod_mat_vec, tailleligneMatrice = " <<nb_lign
|
|
<< " taillevecteur = " << vec.Taille();
|
|
cout << " Vecteur& MatLapack::Prod_mat_vec (... " << endl;
|
|
Sortie (1);
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n Vecteur& MatLapack::Prod_mat_vec (... " << endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
#ifdef MISE_AU_POINT
|
|
// il faut que les lignes et vec aient la même dimension
|
|
if (vec.Taille() != Ligne(1).Taille())
|
|
{cout << "erreur de taille, la première ligne de la matrice n'a pas la même dimension";
|
|
cout << " que le vecteur à multiplier "
|
|
<< "\n MatLapack::Prod_mat_vec (const Vecteur& vec,Vecteur& result)" << endl;
|
|
Sortie(1);
|
|
}
|
|
#endif
|
|
// for (int i=1;i<=nb_lign;i++)
|
|
// {double* pt= &result(i);
|
|
// for (int k=1;k<=nb_col;k++)
|
|
// (*pt) += (this->*Coor_const)(i,k) * vec(k);
|
|
// };
|
|
|
|
|
|
// y := alpha*A*x + beta*y,
|
|
// avec ici beta = 0.,
|
|
double alpha=1.;double beta = 0.;
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ //void cblas_dgemv(const enum CBLAS_ORDER Order,
|
|
// const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
|
|
// const double alpha, const double *A, const int lda,
|
|
// const double *X, const int incX, const double beta,
|
|
// double *Y, const int incY) __OSX_AVAILABLE_STARTING(__MAC_10_2,__IPHONE_4_0);
|
|
cblas_dgemv(CblasRowMajor,CblasNoTrans, nb_lign, lda,alpha,mat.v,lda,vec.v,1,beta,result.v,1);
|
|
//for (int i=1;i<=nb_lign;i++)
|
|
// {double* pt= &result(i);
|
|
// for (int k=1;k<=nb_col;k++)
|
|
// (*pt) += (this->*Coor_const)(i,k) * vec(k);
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ //void cblas_dgbmv(const enum CBLAS_ORDER Order,
|
|
// const enum CBLAS_TRANSPOSE TransA, const int M, const int N,
|
|
// const int KL, const int KU, const double alpha,
|
|
// const double *A, const int lda, const double *X,
|
|
// const int incX, const double beta, double *Y, const int incY)
|
|
|
|
cblas_dgbmv(CblasRowMajor,CblasNoTrans, nb_lign, lda,kl,ku,alpha,mat.v,lda,vec.v,1,beta,result.v,1);
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
{ //void cblas_dsbmv(const enum CBLAS_ORDER Order, const enum CBLAS_UPLO Uplo,
|
|
// const int N, const int K, const double alpha, const double *A,
|
|
// const int lda, const double *X, const int incX,
|
|
// const double beta, double *Y, const int incY) __OSX_AVAILABLE_STARTING(__MAC_10_2,
|
|
|
|
cblas_dsbmv(CblasRowMajor, CblasUpper,nb_lign,ku,alpha,mat.v,lda,vec.v,1,beta,result.v,1);
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n **** erreur: mauvais choix de matrice: " << Nom_matrice(type_matrice)
|
|
<< "\n MatLapack::Simple_Resol_systID_2 (...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
|
|
return result;
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Retourne la transposee de la matrice mat
|
|
// il y a création d'une matrice de sortie
|
|
MatLapack MatLapack::Transpose () const
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if ( ( nb_lign ==0 ) || ( nb_col ==0 ) )
|
|
{
|
|
cout << "Erreur : une des dimensions de la matrice:NxM= " << nb_lign <<"x"<< nb_col
|
|
<< " est nulle, la transposee n'a pas de sens ici ";
|
|
cout << "\n MatLapack::TRANSPOSE ( )\n";
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
// comme résolution et type de matrice on met la même chose que this
|
|
// for (int i=1;i <= nb_col;i++)
|
|
// for (int j=1;j <= nb_lign;j++)
|
|
// (result.*Coor)(i,j)=(this->*Coor_const)(j,i);
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ // dans le cas d'une matrice carrée ou rectangulaire on inverse les dimensions
|
|
// car on considère que le stockage décri la forme réelle de la matrice
|
|
MatLapack result(nb_col,nb_lign
|
|
,Symetrique_Enum_matrice(this->type_matrice),0.
|
|
,this->type_resolution,this->type_preconditionnement);
|
|
for (int i=1;i <= nb_col;i++)
|
|
for (int j=2;j <= nb_lign;j++) // on commence à 2, car la diag ne change pas
|
|
(result.*Coor)(i,j)=(this->*Coor_const)(j,i);
|
|
return result;
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ // on considère que la matrice réelle est carrée, et que son stockage bande n'est là
|
|
// que pour l'optimisation. Donc la matrice transposée est exactement de même type
|
|
// même largeur de bande etc, et dans notre cas on considère que la bande est de même 1/2 largeur
|
|
// dessus et dessous, donc kl = ku
|
|
MatLapack result(*this);
|
|
for (int i=1;i <= nb_col;i++)
|
|
for (int j=2;j <= nb_lign;j++) // on commence à 2, car la diag ne change pas
|
|
(result.*Coor)(i,j)=(this->*Coor_const)(j,i);
|
|
return result;
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:case CARREE_SYMETRIQUE_LAPACK:
|
|
{ // on considère que la matrice réelle est carrée, et que son stockage bande est symétrique supérieur
|
|
// donc la transposée est identique
|
|
MatLapack result(*this);
|
|
return result;
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n **** erreur: mauvais choix de matrice: " << Nom_matrice(type_matrice)
|
|
<< "\n MatLapack::Simple_Resol_systID_2 (...";
|
|
Sortie(1);
|
|
};
|
|
};
|
|
MatLapack toto;
|
|
return toto; // pour taire le compilo
|
|
};
|
|
|
|
// determine l'inverse d'une matrice carre
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
MatLapack MatLapack::Inverse()
|
|
{switch (type_resolution)
|
|
{case GAUSS:
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
// l'appel n'est possible que si la matrice est carré
|
|
if (nb_lign != nb_col)
|
|
{cout << "\n erreur la matrice n'est pas carree, resolution impossible "
|
|
<< " nb_ligne= " << nb_lign << " nb_colonne= " << nb_col
|
|
<< "\n MatLapack::Inverse(..." << endl;
|
|
Sortie(1);
|
|
}
|
|
// l'appel n'est pas possible si la matrice a des coefficients
|
|
if (nb_lign == 0)
|
|
{cout << "\n erreur la matrice n'a pas de coefficient, resolution impossible "
|
|
<< " nb_ligne= " << nb_lign << " nb_colonne= " << nb_col
|
|
<< "\n MatLapack::Inverse(..." << endl;
|
|
Sortie(1);
|
|
}
|
|
#endif
|
|
|
|
ENTIER_POUR_CLAPACK info = 0; // info de sortie en fait c'est le retour de la fonction
|
|
// init général
|
|
if (ipiv == NULL) {ipiv = new ENTIER_POUR_CLAPACK [nb_lign];};
|
|
MatLapack Msortie(*this);
|
|
// def de la taille adoc pour le vecteur de travail
|
|
// au premier appel on défini la taille
|
|
ENTIER_POUR_CLAPACK lwork = -1; // pour demander la taille a doc de la zone de travail
|
|
dgetri_(&Msortie.nb_lign,Msortie.mat.v,&lda,ipiv,travail.v,&lwork,&info);
|
|
GestionErreurInverse(info,1,"dgetri_"); // gestion d'erreurs éventuelles
|
|
int dim_optimal=travail(1); // sgetri sauvegarde la taille optimale dans v[0]
|
|
travail.Change_taille(dim_optimal);
|
|
// appel de la triangulation
|
|
dgetrf_(&Msortie.nb_lign,&Msortie.nb_col,Msortie.mat.v,&lda,ipiv,&info);
|
|
GestionErreurInverse(info,2,"dgetrf_"); // gestion d'erreurs éventuelles
|
|
// calcul de l'inverse
|
|
lwork=MaX(1,dim_optimal);
|
|
dgetri_(&Msortie.nb_lign,Msortie.mat.v,&lda,ipiv,travail.v,&lwork,&info);
|
|
GestionErreurInverse(info,3,"dgetri_"); // gestion d'erreurs éventuelles
|
|
// retour
|
|
return Msortie;
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n erreur , methode non implante pour la matrice lapack: "
|
|
<< Nom_resolution(type_resolution)
|
|
<< "\n MatLapack::Inverse()"<< endl;
|
|
Sortie(1);
|
|
};
|
|
};
|
|
// pour taire le warning
|
|
return MatLapack();
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// determine l'inverse d'une matrice carre
|
|
// mais en retournant l'inverse dans la matrice passée en paramètre
|
|
// qui doit être de même taille
|
|
MatLapack& MatLapack::Inverse(MatLapack& Msortie)
|
|
{switch (type_resolution)
|
|
{case GAUSS:
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
// l'appel n'est possible que si la matrice est carré
|
|
if (nb_lign != nb_col)
|
|
{cout << "\n erreur la matrice n'est pas carree, resolution impossible "
|
|
<< " nb_ligne= " << nb_lign << " nb_colonne= " << nb_col
|
|
<< "\n MatLapack::Inverse(..." << endl;
|
|
Sortie(1);
|
|
};
|
|
// l'appel n'est pas possible si la matrice a des coefficients
|
|
if (nb_lign == 0)
|
|
{cout << "\n erreur la matrice n'a pas de coefficient, resolution impossible "
|
|
<< " nb_ligne= " << nb_lign << " nb_colonne= " << nb_col
|
|
<< "\n MatLapack::Inverse(..." << endl;
|
|
Sortie(1);
|
|
};
|
|
if (nb_lign != Msortie.nb_lign)
|
|
{cout << "\n erreur la matrice resultat n'a pas le meme nombre de ligne ("<< Msortie.nb_lign
|
|
<< ") que this (" << nb_lign << "), resolution impossible "
|
|
<< "\n MatLapack::Inverse(..." << endl;
|
|
Sortie(1);
|
|
}
|
|
else if (nb_col != Msortie.nb_col)
|
|
{cout << "\n erreur la matrice resultat n'a pas le meme nombre de colonne ("<< Msortie.nb_col
|
|
<< ") que this (" << nb_col << "), resolution impossible "
|
|
<< "\n MatLapack::Inverse(..." << endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
|
|
ENTIER_POUR_CLAPACK info = 0; // info de sortie en fait c'est le retour de la fonction
|
|
// init général
|
|
if (ipiv == NULL) {ipiv = new ENTIER_POUR_CLAPACK [nb_lign];};
|
|
Msortie = (*this);
|
|
// def de la taille adoc pour le vecteur de travail
|
|
// au premier appel on défini la taille
|
|
ENTIER_POUR_CLAPACK lwork = -1; // pour demander la taille a doc de la zone de travail
|
|
dgetri_(&Msortie.nb_lign,Msortie.mat.v,&lda,ipiv,travail.v,&lwork,&info);
|
|
GestionErreurInverse(info,1,"dgetri_"); // gestion d'erreurs éventuelles
|
|
int dim_optimal=travail(1); // sgetri sauvegarde la taille optimale dans v[0]
|
|
travail.Change_taille(dim_optimal);
|
|
// appel de la triangulation
|
|
dgetrf_(&Msortie.nb_lign,&Msortie.nb_col,Msortie.mat.v,&lda,ipiv,&info);
|
|
GestionErreurInverse(info,2,"dgetrf_"); // gestion d'erreurs éventuelles
|
|
// calcul de l'inverse
|
|
lwork=MaX(1,dim_optimal);
|
|
dgetri_(&Msortie.nb_lign,Msortie.mat.v,&lda,ipiv,travail.v,&lwork,&info);
|
|
GestionErreurInverse(info,3,"dgetri_"); // gestion d'erreurs éventuelles
|
|
// retour
|
|
return Msortie;
|
|
break;
|
|
}
|
|
default:
|
|
{ cout << "\n erreur , methode non implante pour la matrice lapack: "
|
|
<< Nom_resolution(type_resolution)
|
|
<< "\n MatLapack::Inverse()"<< endl;
|
|
Sortie(1);
|
|
};
|
|
};
|
|
// pour taire le warning
|
|
return Msortie;
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Surcharge de l'operateur + : addition de deux matrices
|
|
// il y a création d'une matrice de sortie
|
|
MatLapack MatLapack::operator + (const MatLapack& mat_pl) const
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if ( ( nb_lign!=mat_pl.nb_lign )|| ( nb_col != mat_pl.nb_col ) )
|
|
{ cout << "Erreur : dimensions des matrices non egales\n"
|
|
<< " nb_ligne= " << nb_lign << " " << mat_pl.nb_lign
|
|
<< " nb_colonne= " << nb_col << " " << mat_pl.nb_col << " ";
|
|
cout << "MatLapack::OPERATOR+ (const MatLapack& )\n";
|
|
Sortie(1);
|
|
};
|
|
if (type_matrice != mat_pl.type_matrice)
|
|
{cout << " erreur les matrices sont de types différents "
|
|
<< Nom_matrice(type_matrice) << " " << Nom_matrice(mat_pl.type_matrice) << '\n';
|
|
cout << "MatLapack::operator+ (const MatLapack& mat_pl)" << endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
MatLapack result(*this);
|
|
result.mat +=mat_pl.mat;
|
|
return result;
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Surcharge de l'operateur - : soustraction entre deux matrices
|
|
MatLapack MatLapack::operator- (const MatLapack& mat_pl) const
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if ( ( nb_lign!=mat_pl.nb_lign )|| ( nb_col != mat_pl.nb_col ) )
|
|
{ cout << "Erreur : dimensions des matrices non egales\n"
|
|
<< " nb_ligne= " << nb_lign << " " << mat_pl.nb_lign
|
|
<< " nb_colonne= " << nb_col << " " << mat_pl.nb_col << " ";
|
|
cout << "MatLapack::OPERATOR- (const MatLapack& )\n";
|
|
Sortie(1);
|
|
};
|
|
if (type_matrice != mat_pl.type_matrice)
|
|
{cout << " erreur les matrices sont de types différents "
|
|
<< Nom_matrice(type_matrice) << " " << Nom_matrice(mat_pl.type_matrice) << '\n';
|
|
cout << "MatLapack::operator- (const MatLapack& mat_pl)" << endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
MatLapack result(*this);
|
|
// initialisation
|
|
result.mat -= mat_pl.mat;
|
|
return result;
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Surcharge de l'operateur - : oppose d'une matrice
|
|
MatLapack MatLapack::operator- () const
|
|
{ MatLapack result(*this);
|
|
// initialisation
|
|
result.mat = - mat;
|
|
return result;
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
MatLapack MatLapack::operator* (const MatLapack& mat_pl) const
|
|
// Surcharge de l'operateur * : multiplication de deux matrices
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if ( nb_col!= mat_pl.nb_lign )
|
|
{ cout << "Erreur : dimensions des matrices non valables, "
|
|
<< " pour la multiplication le nombre de colonne de la premiere matrice "
|
|
<< " doit etre egale au nombre de ligne de la seconde alors qu'ici on a respectivement: "
|
|
<< nb_col << " " << mat_pl.nb_lign << "\n";
|
|
cout << "MatLapack::OPERATOR* (const MatLapack& )\n";
|
|
Sortie(1);
|
|
};
|
|
if (type_matrice != mat_pl.type_matrice)
|
|
{ cout << " erreur les matrices sont de types différents "
|
|
<< Nom_matrice(type_matrice) << " " << Nom_matrice(mat_pl.type_matrice) << '\n';
|
|
cout << "MatLapack::operator* (const MatLapack& mat_pl)" << endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ // par défaut on initialise la matrice résultat en non symétrique (car le produit de deux matrices symétriques
|
|
// en en général non symétrique, et on ne renseigne pas la résolution
|
|
// la première matrice : n_lign x nb_col
|
|
// la deuxième matrice : mat_pl.n_lign x mat_pl.nb_col
|
|
// la matrice résultat : n_lign x mat_pl.nb_col
|
|
MatLapack result(nb_lign,mat_pl.nb_col,false,0.,RIEN_TYPE_RESOLUTION_MATRI,RIEN_PRECONDITIONNEMENT);
|
|
// utilisation de blas3
|
|
// void cblas_dgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
|
|
// const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
|
|
// const int K, const double alpha, const double *A,
|
|
// const int lda, const double *B, const int ldb,
|
|
// const double beta, double *C, const int ldc) __OSX_AVAILABLE_STARTING(__MAC_10_2,
|
|
// C := alpha*op( A )*op( B ) + beta*C,
|
|
// ici on prend alpha = 1. et beta = 0.
|
|
double alpha = 1.; double beta = 0.;
|
|
// MatLapack inter(nb_lign,mat_pl.nb_col,false,0.,RIEN_TYPE_RESOLUTION_MATRI,RIEN_PRECONDITIONNEMENT);
|
|
// le mat.v de la fin, ne sert à rien car beta = 0
|
|
cblas_dgemm(CblasRowMajor, CblasNoTrans,CblasNoTrans,nb_lign,mat_pl.nb_col,nb_col,alpha
|
|
,mat.v,nb_lign,mat_pl.mat.v,mat_pl.nb_lign,beta,result.mat.v,nb_col);
|
|
|
|
// for (int i=1;i<=nb_lign;i++)
|
|
// { for (int j=1;j<= nb_col;j++)
|
|
// for (int k=1;k<=nb_col;k++)
|
|
// (result.*Coor)(i,j) += (this->*Coor_const)(i,k) * (mat_pl.*Coor_const)(k,j);
|
|
// };
|
|
return result;
|
|
break;
|
|
}
|
|
// effectivement la suite n'est pas intéressante et ne sera pas utilisée, donc on laisse l'erreur s'afficher
|
|
// case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
// { // dans le cas du produit de deux matrices bandes, les bandes s'ajoutent
|
|
// MatLapack result(type_matrice,2*nb_col,nb_lign,false,0.,RIEN_TYPE_RESOLUTION_MATRI,RIEN_PRECONDITIONNEMENT);
|
|
// // mais pour l'instant je ne suis pas sûr que cela serve à quelque chose donc je met un message d'erreur
|
|
// cout << "\n erreur **** le produit de deux matrices bandes n'est pas implante "
|
|
// << "\n MatLapack::operator* ( ... " << endl;
|
|
// Sortie(1);
|
|
// return result;
|
|
// break;
|
|
// }
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::Colonne (.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
|
|
MatLapack toto;
|
|
return toto; // pour taire le compilo
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Surcharge de l'operateur * : multiplication d'une matrice par un scalaire
|
|
MatLapack MatLapack::operator* (const double& coeff) const
|
|
{ MatLapack result(*this);
|
|
// on n'utilise pas blas, car cela devra être fait dans la classe vecteur
|
|
// et sera donc directement répercuté ici
|
|
result *=coeff;
|
|
return result;
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Surcharge de l'operateur == : test d'egalite entre deux matrices
|
|
// Retourne 1 si les deux matrices sont egales
|
|
// 0 sinon
|
|
int MatLapack::operator== (const MatLapack& mat_pl) const
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
// dans le cas où les matrices ne sont pas de caractéristiques identiques
|
|
// il y a un message d'erreur
|
|
if (type_matrice != mat_pl.Type_matrice())
|
|
{cout << " erreur les matrices sont de types différents "
|
|
<< Nom_matrice(type_matrice) << " " << Nom_matrice(mat_pl.type_matrice) << '\n';
|
|
cout << "int MatLapack::operator== (..." << endl;
|
|
Sortie(1);
|
|
}
|
|
#endif
|
|
if ( ( nb_lign != mat_pl.nb_lign ) || ( nb_col != mat_pl.nb_col ) )
|
|
return 0;
|
|
if (mat != mat_pl.mat)
|
|
return 0;
|
|
return 1;
|
|
};
|
|
|
|
//met une valeur identique sur toute la ligne
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
void MatLapack::MetValLigne(int i,double x)
|
|
{ switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:case CARREE_SYMETRIQUE_LAPACK:
|
|
{ int le_max = nb_col+1;// +1 pour un test "inférieur pure" sur la ligne qui suit
|
|
for(int j=1; j < le_max; j++)
|
|
(this->*Coor)(i,j)=x;
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ int lb=(nb_col-1)/2; int dim = nb_lign;
|
|
int maxJ = MiN(dim,i+lb)+1;// +1 pour un test "inférieur pure" sur la ligne qui suit
|
|
for (int j= MaX(1,i-lb);j< maxJ;j++)
|
|
(this->*Coor)(i,j) = x;
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK:
|
|
// même symétrique il faut imposer sur toute la ligne car au-dessus et en-dessous de la diagonale,
|
|
// on n'adresse pas des grandeurs identiques
|
|
{int lb=nb_col-1; int dim = nb_lign;
|
|
int le_max = MiN(dim,i+lb)+1;// +1 pour un test "inférieur pure" sur la ligne qui suit
|
|
int le_min = MaX(1,i-lb); // on fait comme pour MatBand
|
|
// for (int j= MaX(1,i);j< le_max;j++)
|
|
for (int j= le_min;j< le_max;j++)
|
|
(this->*Coor)(i,j)=x;
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::MetValLigne (.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// Retourne la jeme colonne de la matrice
|
|
Vecteur MatLapack::Colonne(int j) const
|
|
// j est entre 1 et dim
|
|
{ Vecteur ret(nb_lign);
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:case CARREE_SYMETRIQUE_LAPACK:
|
|
{ int le_max = nb_lign+1;// +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
for (int i=1;i<le_max;i++)
|
|
ret(i)=(this->*Coor_const)(i,j);
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ int lb=(nb_col-1)/2; int dim = nb_lign;
|
|
int maxI = MiN(dim,j+lb)+1;// +1 pour un test "inférieur pure" sur la ligne qui suit
|
|
for (int i= MaX(1,j-lb);i< maxI;i++)
|
|
ret(i) = (this->*Coor_const)(i,j);
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK: // bande supérieure
|
|
{int lb=nb_col-1; int dim = nb_lign;
|
|
int le_max = MiN(dim,j+lb)+1;// +1 pour un test "inférieur pure" sur la boucle qui suit
|
|
int le_min = MaX(1,j-lb);
|
|
for (int i= le_min;i< le_max;i++)
|
|
ret(i) = (this->*Coor_const)(i,j);
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::Colonne (.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
return ret;
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// 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 MatLapack::ColonneSpe(int j) const
|
|
{ switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ Vecteur ret(nb_lign);
|
|
for (int i=1;i<=nb_lign;i++)
|
|
ret(i)=(this->*Coor_const)(i,j);
|
|
return ret;
|
|
break;
|
|
}
|
|
case CARREE_SYMETRIQUE_LAPACK:
|
|
// on fait pareil pour la matrice symétrique, car la colonne complète adresse
|
|
// deux parties différentes de symétrie, donc il faut récupérer toute la colonne
|
|
{ Vecteur ret(nb_lign);
|
|
for (int i=1;i<=nb_lign;i++)
|
|
ret(i)=(this->*Coor_const)(i,j);
|
|
return ret;
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ Vecteur col(nb_col); // mise a zero de ligne, nb_col = largeur totale de la bande
|
|
int lb=(nb_col-1)/2; int dim = nb_lign;
|
|
int maxI = MiN(dim,j+lb)+1; // +1 pour le test < sur le for qui suit
|
|
for (int i= MaX(1,j-lb);i< maxI;i++)
|
|
col(i-j+lb+1) = (this->*Coor_const)(i,j);
|
|
return col;
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK: // bande supérieure
|
|
{int lb=nb_col-1; int dim = nb_lign;
|
|
Vecteur col(1+2*lb); // mise a zero de ligne,
|
|
// méthode identique au cas de la matrice symétrique, car la colonne complète adresse
|
|
// deux parties différentes de symétrie, donc il faut récupérer toute la colonne
|
|
int le_max = MiN(dim,j+lb)+1;// +1 pour le test < sur le for qui suit
|
|
int le_min = MaX(1,j-lb);
|
|
for (int i= le_min;i < le_max;i++)
|
|
col(i-j+lb+1) = (this->*Coor_const)(i,j);
|
|
return col;
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::ColonneSpe (.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
Vecteur toto;
|
|
return toto; // pour taire le compilo
|
|
};
|
|
|
|
// remplace la Colonne de la matrice par la colonne fournie
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
void MatLapack::RemplaceColonneSpe(int j, const Vecteur & v)
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ if (nb_lign != v.Taille())
|
|
{ cout << " \nla taille de la colonne j= " << j <<" a remplacer n est pas bonne " ;
|
|
cout << " taille = " << v.Taille() << " tailleColMatrice = " << nb_lign;
|
|
cout << " void MatLapack::RemplaceColonne(int j,Vecteur & v) " << endl;
|
|
Sortie (1);
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK: case BANDE_SYMETRIQUE_LAPACK:
|
|
{ if (v.Taille() != nb_col)
|
|
{cout << "\nerreur d affectation pour le vecteur";
|
|
cout << " dim vecteur = " << v.Taille() << " au lieu de " << nb_col << '\n';
|
|
cout << "void MatLapack::RemplaceColonne(int j,Vecteur & v)" << endl;
|
|
};
|
|
Sortie (1);
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::RemplaceColonneSpe (.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:
|
|
{ for(int i=1; i <= nb_lign ; i++)
|
|
(this->*Coor)(i,j)=v(i);
|
|
break;
|
|
}
|
|
case CARREE_SYMETRIQUE_LAPACK:
|
|
// méthode identique au cas de la matrice non symétrique, car la colonne complète adresse
|
|
// deux parties différentes de symétrie, donc il faut récupérer toute la colonne
|
|
{ for(int i=1; i <= nb_lign ; i++)
|
|
(this->*Coor)(i,j)=v(i);
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{ int lb=(nb_col-1)/2; int dim = nb_lign;
|
|
int maxI = MiN(dim,j+lb)+1;// +1 pour inférieur pure sur la ligne qui suit
|
|
for (int i= MaX(1,j-lb);i< maxI;i++)
|
|
(this->*Coor)(i,j)=v(i);
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK: // bande supérieure
|
|
{int lb=nb_col-1; int dim = nb_lign;
|
|
Vecteur col(1+2*lb); // mise a zero de ligne,
|
|
// méthode identique au cas de la matrice non symétrique, car la colonne complète adresse
|
|
// deux parties différentes de symétrie, donc il faut récupérer toute la colonne
|
|
int le_max = MiN(dim,j+lb)+1;// +1 pour le test < sur le for qui suit
|
|
int le_min = MaX(1,j-lb);
|
|
for (int i= le_min;i<le_max;i++)
|
|
(this->*Coor)(i,j) = col(i-j+lb+1);
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::RemplaceColonneSpe (.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
};
|
|
//met une valeur identique sur toute la colonne
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
void MatLapack::MetValColonne(int j,double y)
|
|
{ switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:case CARREE_SYMETRIQUE_LAPACK:
|
|
// méthode identique entre symétrique et non symétrique, car pour le stockage symétrique
|
|
// la colonne complète adresse
|
|
// deux parties différentes de symétrie, donc il faut récupérer toute la colonne
|
|
{int le_maxi = nb_lign+1;
|
|
for(int i=1; i < le_maxi ; i++)
|
|
(this->*Coor)(i,j)=y;
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:
|
|
{int lb=(nb_col-1)/2; int dim = nb_lign;
|
|
int maxI = MiN(dim,j+lb)+1;// +1 pour le test < sur le for qui suit
|
|
for (int i= MaX(1,j-lb);i< maxI;i++)
|
|
(this->*Coor)(i,j)=y;
|
|
break;
|
|
}
|
|
case BANDE_SYMETRIQUE_LAPACK: // bande supérieure
|
|
// méthode identique au cas de la matrice non symétrique, car la colonne complète adresse
|
|
// deux parties différentes de symétrie, donc il faut récupérer toute la colonne
|
|
{int lb=nb_col-1; int dim = nb_lign;
|
|
int maxI = MiN(dim,j+lb)+1;// +1 pour le test < sur le for qui suit
|
|
int le_min = MaX(1,j-lb); // on fait comme MatBand
|
|
for (int i= le_min;i<maxI;i++)
|
|
(this->*Coor)(i,j) = y;
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::MetValColonne (.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
};
|
|
|
|
// Multiplication d'une ligne iligne de la matrice avec un vecteur de
|
|
// dimension = le nombre de colonne de la matrice
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
double MatLapack::Prod_Ligne_vec ( int iligne, const Vecteur& vec) const
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK: case CARREE_SYMETRIQUE_LAPACK:
|
|
{ if (nb_col != vec.Taille())
|
|
{ cout << " \n erreur la taille de la ligne i= " << iligne <<" n est pas correcte ";
|
|
cout << " tailleligneMatrice = " <<nb_col
|
|
<< " taillevecteur = " << vec.Taille();
|
|
cout << " double MatLapack::Prod_Ligne_vec ( int iligne,Vecteur& vec) " << endl;
|
|
Sortie (1);
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK: case BANDE_SYMETRIQUE_LAPACK:
|
|
{ if (vec.Taille() != nb_lign)
|
|
{cout << "\nerreur d affectation pour le vecteur";
|
|
cout << " dim vecteur = " << vec.Taille() << " au lieu de " << nb_lign << '\n';
|
|
cout << "void MatLapack::Prod_Ligne_vec(..." << endl;
|
|
};
|
|
Sortie (1);
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::Prod_Ligne_vec (.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
double res = 0;
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK: case CARREE_SYMETRIQUE_LAPACK:
|
|
{ for (int j=1;j <= nb_col;j++)
|
|
res += (this->*Coor_const)(iligne,j) * vec(j);
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:case BANDE_SYMETRIQUE_LAPACK:
|
|
{int lb=(nb_col-1)/2; int dim = nb_lign;
|
|
int maxJ = MiN(dim,iligne+lb)+1;// +1 pour le test < sur le for qui suit
|
|
for (int j=MaX(1,iligne-lb);j< maxJ;j++)
|
|
res += (this->*Coor_const)(iligne,j) * vec(j);
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::Prod_Ligne_vec (.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
return res;
|
|
};
|
|
|
|
// Multiplication d'un vecteur avec une colonne icol de la matrice
|
|
// dimension = le nombre de ligne de la matrice
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
double MatLapack::Prod_vec_col( int jcol, const Vecteur& vec) const
|
|
{
|
|
#ifdef MISE_AU_POINT
|
|
if (nb_lign != vec.Taille())
|
|
{ cout << " \n erreur la taille de la colonne j= " << jcol <<" n est pas correcte ";
|
|
cout << " tailleColMatrice = " <<nb_lign
|
|
<< " taillevecteur = " << vec.Taille();
|
|
cout << " double MatLapack::Prod_vec_col( int jcol,Vecteur& vec) " << endl;
|
|
Sortie (1);
|
|
};
|
|
#endif
|
|
double res = 0;
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK: case CARREE_SYMETRIQUE_LAPACK:
|
|
{for (int i=1;i <= nb_lign;i++)
|
|
res += vec(i) * (this->*Coor_const)(i,jcol) ;
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK: case BANDE_SYMETRIQUE_LAPACK:
|
|
{int lb=(nb_col-1)/2; int dim = nb_lign;
|
|
int maxI = MiN(jcol+lb,dim)+1;// +1 pour le test < sur le for qui suit
|
|
for (int i=MaX(1,jcol-lb);i< maxI;i++)
|
|
res += vec(i) * (this->*Coor_const)(i,jcol) ;
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::Prod_vec_col (.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
return res;
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// calcul du produit : (vec_1)^T * A * (vect_2)
|
|
double MatLapack::vectT_mat_vec(const Vecteur& vec1, const Vecteur& vec2) const
|
|
{ double resu=0.; // valeur de retour
|
|
#ifdef MISE_AU_POINT
|
|
if (vec1.Taille() != nb_lign)
|
|
{cout << "erreur de taille, la dimension de (vec1)^T= " << vec1.Taille() << " n'a pas la même dimension";
|
|
cout << " que le le nombre de ligne de la matrice= " << nb_lign
|
|
<< "\n MatLapack::vectT_mat_vec(const Vecteur& vec1, const Vecteur& vec2)" << endl;
|
|
Sortie(1);
|
|
}
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:case CARREE_SYMETRIQUE_LAPACK:
|
|
{ if (vec2.Taille() != nb_col)
|
|
{ cout << "erreur de taille, la dimension de (vec2)= " << vec2.Taille() << " n'a pas la même dimension";
|
|
cout << " que le le nombre de colonne de la matrice= " << nb_col
|
|
<< "\n MatLapack::vectT_mat_vec(const Vecteur& vec1, const Vecteur& vec2)" << endl;
|
|
Sortie(1);
|
|
};
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK:case BANDE_SYMETRIQUE_LAPACK:
|
|
{ if (vec2.Taille() != nb_lign)
|
|
{ cout << "erreur de taille, la dimension de (vec2)= " << vec2.Taille() << " n'a pas la même dimension";
|
|
cout << " que le le nombre de colonne de la matrice= " << nb_lign
|
|
<< "\n MatLapack::vectT_mat_vec(const Vecteur& vec1, const Vecteur& vec2)" << endl;
|
|
Sortie(1);
|
|
};
|
|
break;
|
|
};
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::vectT_mat_vec (.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
#endif
|
|
|
|
switch (type_matrice)
|
|
{ case RECTANGLE_LAPACK: case CARREE_LAPACK:case CARREE_SYMETRIQUE_LAPACK:
|
|
{ for (int i=1;i<=nb_lign;i++)
|
|
for (int j= 1;j<= nb_col;j++)
|
|
resu += vec1(i) * (this->*Coor_const)(i,j) * vec2(j);
|
|
break;
|
|
}
|
|
case BANDE_NON_SYMETRIQUE_LAPACK: case BANDE_SYMETRIQUE_LAPACK:
|
|
{ int lb=(nb_col-1)/2; int dim = nb_lign;
|
|
for (int i=1;i<=dim;i++)
|
|
{ int maxJ = MiN(dim,i+lb)+1;// +1 pour inférieur pure sur la ligne qui suit
|
|
for (int j= MaX(1,i-lb);j< maxJ;j++)
|
|
resu += vec1(i) * (this->*Coor_const)(i,j) * vec2(j);
|
|
};
|
|
break;
|
|
}
|
|
default:
|
|
cout << "\n erreur: cas de type de matrice: "<< Nom_matrice(type_matrice)
|
|
<< " non implante "
|
|
<< "\n MatLapack::vectT_mat_vec (.. "<<endl;
|
|
Sortie(1);
|
|
};
|
|
return resu;
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// retourne la place que prend la matrice en entier
|
|
int MatLapack::Place() const
|
|
{switch (type_matrice)
|
|
{ case CARREE_SYMETRIQUE_LAPACK: return nb_lign*(nb_col-1)/2+nb_col;
|
|
case CARREE_LAPACK: return nb_lign*nb_col;
|
|
case RECTANGLE_LAPACK: return (2*(nb_lign * nb_col));break;
|
|
case BANDE_SYMETRIQUE_LAPACK: return lda*ldb;
|
|
case BANDE_NON_SYMETRIQUE_LAPACK: return (2*((nb_lign+(nb_lign-1)/2) * nb_col));break;
|
|
default: return 0;
|
|
};
|
|
};
|
|
|
|
|
|
// =================== méthodes protégées =======================
|
|
/*
|
|
// 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.
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
void MatLapack::Triangulation (int N,MatLapack& B)
|
|
|
|
{ double somme;
|
|
|
|
// Determination des composantes de la matrice
|
|
// triangulaire inferieure B :
|
|
for (int k=1;k<=N;k++)
|
|
{
|
|
somme=0;
|
|
for (int j=1;j<=k-1;j++)
|
|
{
|
|
// somme=somme+ B(k,j)*B(k,j); // pow(B(k,j),2.0);
|
|
somme += B(k,j)*B(k,j); // pow(B(k,j),2.0);
|
|
};
|
|
B(k,k)=sqrt((*this)(k,k)-somme);
|
|
for (int i=k+1;i<=N;i++)
|
|
{
|
|
somme=0;
|
|
for (int j=1;j<=k-1;j++)
|
|
{
|
|
// somme=somme+(B(i,j)*B(k,j));
|
|
somme += (B(i,j)*B(k,j));
|
|
};
|
|
B(i,k)=((*this)(i,k)-somme)/B(k,k);
|
|
};
|
|
};
|
|
|
|
};
|
|
|
|
// résolution du problème triangularisé
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
void MatLapack::Resolution
|
|
(int N, const MatLapack& B,const Vecteur& b,Vecteur& resul) const
|
|
{ // Resolution du systeme B.y=b ( avec y=Bt.x )
|
|
// Les valeurs de y sont stockees dans le vecteur result
|
|
for (int i=1;i<=N;i++)
|
|
{
|
|
resul(i)=b(i);
|
|
// for (int j=1;j<=i-1;j++)
|
|
// resul(i)=resul(i)-B(i,j)*resul(j);
|
|
// resul(i)=resul(i)/B(i,i);
|
|
|
|
for (int j=1;j<=i-1;j++)
|
|
resul(i) -= B(i,j)*resul(j);
|
|
resul(i) /= B(i,i);
|
|
};
|
|
|
|
// Resolution du systeme Bt.x=y
|
|
// Les valeurs de x sont stockees dans le vecteur result
|
|
// B=B.Transpose();
|
|
// for (int i=N;i>=1;i--)
|
|
// {
|
|
// for (int j=i+1;j<=N;j++)
|
|
// result(i)=result(i)-B(j,i)*result(j);
|
|
// result(i)=result(i)/B(i,i);
|
|
// };
|
|
for (int i=N;i>=1;i--)
|
|
{
|
|
for (int j=i+1;j<=N;j++)
|
|
resul(i) -= B(j,i)*resul(j);
|
|
// resul(i)=resul(i)-B(j,i)*resul(j);
|
|
// resul(i)=resul(i)/B(i,i);
|
|
resul(i) /= B(i,i);
|
|
};
|
|
};
|
|
|
|
*/
|
|
// surcharge de l'operateur de lecture typée
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
istream & operator >> (istream & entree, MatLapack & mat_pl)
|
|
{ // vérification du type
|
|
string typa;
|
|
entree >> typa;
|
|
if (typa != "MatLapack_NxM")
|
|
{cout << "\n erreur en lecture d'une matrice type lapack on attendait la chaine: MatLapack_NxM"
|
|
<< " et on a lue: " << typa
|
|
<< "\n operator >> (istream & entree ... ";
|
|
Sortie (1);
|
|
};
|
|
// on lit le nombre de ligne et de colonne
|
|
int nb_l,nb_c;
|
|
entree >> nb_l >> nb_c;
|
|
// on lit le type et on vérifie que la matrice à le même type que this
|
|
entree >> typa;
|
|
if (Id_nom_matrice(typa.c_str()) != mat_pl.type_matrice)
|
|
{cout << "\n erreur en lecture du type de la matrice on attendait " << Nom_matrice(mat_pl.type_matrice)
|
|
<< " et on a lue: " << typa
|
|
<< "\n operator >> (istream & entree ... ";
|
|
Sortie (1);
|
|
};
|
|
// on redimentionne en fonction du type
|
|
|
|
if ((nb_l != mat_pl.nb_lign) || (nb_c != mat_pl.nb_col))
|
|
{ mat_pl.Change_taille(nb_l,nb_c);}
|
|
else
|
|
{ // même si la place est ok on met à jour le nombre de ligne et de colonne
|
|
mat_pl.nb_lign=nb_l;mat_pl.nb_col=nb_c;
|
|
};
|
|
|
|
// maintenant on peut lire les coordonnées
|
|
entree >> typa; // pour passer le mot clé valeurs=
|
|
int nb=mat_pl.mat.Taille();
|
|
for(int j=1; j <= nb; j++)
|
|
entree >> mat_pl.mat(j);
|
|
return entree;
|
|
};
|
|
|
|
// surcharge de l'operateur d'ecriture typée
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
ostream & operator << ( ostream & sort,const MatLapack & mat_pl)
|
|
{ // un indicateur donnant le type puis les datas
|
|
sort << "MatLapack_NxM " << mat_pl.nb_lign << " " << mat_pl.nb_col << " "
|
|
<< Nom_matrice(mat_pl.type_matrice) << " valeurs= " ;
|
|
int nb = mat_pl.mat.Taille();
|
|
for(int j=1; j <= nb; j++)
|
|
sort << setprecision(ParaGlob::NbdigdoCA()) << mat_pl.mat(j);
|
|
// retour
|
|
return sort;
|
|
};
|
|
|
|
// fonction interne de vérification des composantes
|
|
// arrêt si pb
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
void MatLapack::Verif_nb_composante(const Vecteur& b,string nom_routine) const
|
|
{
|
|
// l'appel n'est possible que si la matrice est carré
|
|
if (nb_lign != nb_col)
|
|
{cout << "\n erreur la matrice n'est pas carree, resolution impossible "
|
|
<< " nb_ligne= " << nb_lign << " nb_colonne= " << nb_col
|
|
<< "\n MatLapack::" << nom_routine << "(..." << endl;
|
|
Sortie(1);
|
|
}
|
|
// l'appel n'est pas possible si la matrice a des coefficients
|
|
if (nb_lign == 0)
|
|
{cout << "\n erreur la matrice n'a pas de coefficient, resolution impossible "
|
|
<< " nb_ligne= " << nb_lign << " nb_colonne= " << nb_col
|
|
<< "\n MatLapack::" << nom_routine << "(..." << endl;
|
|
Sortie(1);
|
|
}
|
|
// vérif du nombre de composante de b
|
|
if (b.Taille()!=nb_lign)
|
|
{cout << "\n erreur le nombre de composante du second membre: " << b.Taille()
|
|
<< " est differant du nombre de ligne de la matrice: " << nb_lign << " "
|
|
<< "\n MatLapack::" << nom_routine << "(..." << endl;
|
|
Sortie(1);
|
|
}
|
|
};
|
|
|
|
// fonction interne de vérification des composantes
|
|
// arrêt si pb
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
void MatLapack::Verif_nb_composante(const Tableau <Vecteur>& b,string nom_routine) const
|
|
{
|
|
// l'appel n'est possible que si la matrice est carré
|
|
if (nb_lign != nb_col)
|
|
{cout << "\n erreur la matrice n'est pas carree, resolution impossible "
|
|
<< " nb_ligne= " << nb_lign << " nb_colonne= " << nb_col
|
|
<< "\n MatLapack::" << nom_routine << "(..." << endl;
|
|
Sortie(1);
|
|
};
|
|
// l'appel n'est pas possible si la matrice a des coefficients
|
|
if (nb_lign == 0)
|
|
{cout << "\n erreur la matrice n'a pas de coefficient, resolution impossible "
|
|
<< " nb_ligne= " << nb_lign << " nb_colonne= " << nb_col
|
|
<< "\n MatLapack::" << nom_routine << "(..." << endl;
|
|
Sortie(1);
|
|
};
|
|
// vérif du nombre de composante des seconds membres b
|
|
int nb_SM = b.Taille();
|
|
for (int i=1;i<=nb_SM;i++)
|
|
if (b(i).Taille()!=nb_lign)
|
|
{cout << "\n erreur le nombre de composante: " << b(i).Taille() << " du second membre numero: " << i
|
|
<< " est differant du nombre de ligne de la matrice: " << nb_lign << " "
|
|
<< "\n MatLapack::Resol_syst(const Tableau <Vecteur>&..." << endl;
|
|
Sortie(1);
|
|
};
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// gestion erreur dgesv
|
|
void MatLapack::GestionErreurDg_sv(const int& info) const
|
|
{ if (info == 0)
|
|
{ return;}
|
|
else if (info < 0)
|
|
{ cout << "\n pb resolution Dg_sv, le " << -info << " argument de la fonction"
|
|
<< " a une valeur illegale "
|
|
<< "\n arguments: ( N, NRHS, A, LDA, IPIV, B, LDB, INFO ) "
|
|
<< " \n MatLapack::GestionErreurDg_sv(const int& info) ";
|
|
// on génère une exception
|
|
throw ErrResolve_system_lineaire(1);
|
|
Sortie(1);
|
|
}
|
|
else if (info > 0)
|
|
{ cout << "\n pb resolution Dg_sv, U(" << info << "," << info <<") est exactement null"
|
|
<< " (cf doc lapack) la factorisation est complete, mais U est exactement singuliere"
|
|
<< " donc la solution ne peut-etre calculee "
|
|
<< " \n MatLapack::GestionErreurDg_sv(const int& info) ";
|
|
// on génère une exception
|
|
throw ErrResolve_system_lineaire(1);
|
|
// cout << "\n corriger le debug de Matlapack "; // Sortie(1);
|
|
};
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// gestion erreur dpbsv
|
|
void MatLapack::GestionErreurDpbsv(const int& info) const
|
|
{ if (info == 0)
|
|
{ return;}
|
|
else if (info < 0)
|
|
{ cout << "\n pb resolution Dpbsv, le " << -info << " argument de la fonction"
|
|
<< " a une valeur illegale "
|
|
<< "\n arguments: ( UPLO, N, KD, NRHS, AB, LDAB, B, LDB, INFO ) "
|
|
<< " \n MatLapack::GestionErreurDpbsv(const int& info) ";
|
|
// on génère une exception
|
|
throw ErrResolve_system_lineaire(1);
|
|
Sortie(1);
|
|
}
|
|
else if (info > 0)
|
|
{ cout << "\n pb resolution Dpbsv, le mineur prindipal d'ordre (" << info << ") de la matrice est non defini positif"
|
|
<< " (cf doc lapack) la factorisation complete ne peut-etre effectuee, "
|
|
<< " donc la solution ne peut-etre calculee "
|
|
<< " \n MatLapack::GestionErreurDpbsv(const int& info) ";
|
|
// on génère une exception
|
|
throw ErrResolve_system_lineaire(1);
|
|
// cout << "\n corriger le debug de Matlapack "; // Sortie(1);
|
|
};
|
|
};
|
|
|
|
#ifndef MISE_AU_POINT
|
|
inline
|
|
#endif
|
|
// gestion erreur Inverse
|
|
void MatLapack::GestionErreurInverse(const int& info,int cas,const string nom_routine_lapack) const
|
|
{ if (info == 0)
|
|
{ return;}; // cas pas de pb
|
|
// cas de pb éventuel
|
|
if (info < 0)
|
|
{ cout << "\n pb resolution, le " << -info << " argument de la fonction"
|
|
<< " a une valeur illegale (cf doc lapack)";
|
|
}
|
|
else if (info > 0)
|
|
{ cout << "\n pb resolution, U(" << info << "," << info <<") est exactement null"
|
|
<< " (cf doc lapack) la factorisation est complete, mais U est exactement singuliere"
|
|
<< " donc la solution ne peut-etre calculee ";
|
|
};
|
|
switch (cas)
|
|
{ case 2:
|
|
{cout << "\n arguments: ( N, A, LDA, IPIV, WORK, LWORK, INFO ) "; break;}
|
|
case 1: case 3:
|
|
{cout << "\n arguments: ( M, N, A, LDA, IPIV, INFO ) "; break;}
|
|
default:
|
|
{ cout << "\n **** erreur: cas= " << cas <<" non pris en compte"
|
|
<< "\n MatLapack::GestionErreurInverse(...";
|
|
// on génère une exception
|
|
throw ErrResolve_system_lineaire(1);
|
|
Sortie(1);
|
|
};
|
|
};
|
|
cout << " \n routine : " << nom_routine_lapack << " cas = " << cas
|
|
<< "\n MatLapack::GestionErreurInverse(... ";
|
|
// on génère une exception
|
|
throw ErrResolve_system_lineaire(1);
|
|
Sortie(1);
|
|
};
|
|
|
|
|
|
|
|
#endif
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|