// 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