// 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 "Debug.h"
#include "MatBand.h"
#include <iomanip>
#include "CharUtil.h"
#include "MatDiag.h"
#include "ParaGlob.h"


// par defaut
MatBand::MatBand () : 
   Mat_abstraite(BANDE_SYMETRIQUE,CHOLESKY,RIEN_PRECONDITIONNEMENT)
  { lb = 0;
    dim = 0;
    pt = NULL;
   };
   
// def d'une matrice bande
MatBand::MatBand (Enum_matrice type_mat,int lb1, int dim1 ) : 
   Mat_abstraite(type_mat,CHOLESKY,RIEN_PRECONDITIONNEMENT)
  {
    #ifdef MISE_AU_POINT 
     if ((lb1 < 0) || (dim1 < 0))
       {cout << "erreur de dimensionnement d une matrice bande";
        cout << " LB =  " << lb1 << " dim = " << dim1 << '\n';
        cout << "MatBand::MatBand (int lb, int dim )" << endl;
        Sortie(1);
      }
	  if (type_mat != BANDE_NON_SYMETRIQUE)	
       {cout << "erreur de type de matrice, demande de la construction d'une matrice " << Nom_matrice(type_mat) 
		       << " avec l'utilisation d'une MatBand: impossible, MatBand ne convient que pour les matrices symetriques "
				 << " essayer un autre type de matrice !!! ";
		  if (ParaGlob::NiveauImpression() > 2) 
		     cout << "\n MatBand::MatBand (Enum_matrice type_mat,int lb1, int dim1 )";		 
        cout <<  endl;
        Sortie(1);
       }
    #endif  
    lb = lb1;
    dim = dim1;
    try {pt = new double [lb*dim];}
    catch (ErrSortieFinale)
         // cas d'une direction voulue vers la sortie
         // on relance l'interuption pour le niveau supérieur
       { ErrSortieFinale toto;
         throw (toto);
       }
    catch (...)
     {cout << "\n Erreur d'allocation de memoire: taille de reel demandee: "
           << lb <<"*"<<dim <<"= "<< (lb*dim) ;
      if (ParaGlob::NiveauImpression() > 2)
          cout << "\n MatBand::MatBand (... ";
      Sortie(1);
     };
  };
    
// def d'une matrice bande avec 
// initialisation des composantes a la valeur a
MatBand::MatBand (Enum_matrice type_mat, int lb1, int dim1 , double a) :
   Mat_abstraite(type_mat,CHOLESKY,RIEN_PRECONDITIONNEMENT)
  {
    #ifdef MISE_AU_POINT 
    if ((lb1 < 0) || (dim1 < 0))
       {cout << "erreur de dimensionnement d une matrice bande";
        cout << " LB =  " << lb1 << " dim = " << dim1 << '\n';
        cout << "MatBand::MatBand (int lb1, int dim1 , double a)" << endl;
        Sortie(1);
       }
	   if (type_mat != BANDE_SYMETRIQUE)
       {cout << "erreur de type de matrice, demande de la construction d'une matrice " << Nom_matrice(type_mat) 
		           << " avec l'utilisation d'une MatBand: impossible, MatBand ne convient que pour les matrices symetriques "
				         << " essayer un autre type de matrice !!! ";
		      if (ParaGlob::NiveauImpression() > 2)
		        cout << "\n MatBand::MatBand (Enum_matrice type_mat, int lb1, int dim1 , double a) ";
        cout <<  endl;
        Sortie(1);
       }
    #endif    
    lb = lb1;
    dim = dim1;
    try {pt = new double [lb*dim];}
    catch (ErrSortieFinale)
         // cas d'une direction voulue vers la sortie
         // on relance l'interuption pour le niveau supérieur
       { ErrSortieFinale toto;
         throw (toto);
       }
    catch (...)
     {cout << "\n Erreur d'allocation de memoire: taille de reel demandee: "
           << lb <<"*"<<dim <<"= "<< (lb*dim) ;
      if (ParaGlob::NiveauImpression() > 2)
          cout << "\n MatBand::MatBand (... ";
      Sortie(1);
     };
    for (double * pti = pt; pti < pt+lb*dim; pti++)
      *pti = a;
  };
  
// de copie, il y a creation d'une deuxieme matrice     
MatBand::MatBand (const MatBand& m) : 
   Mat_abstraite(m)
  { lb = m.Nb_colonne();
    dim = m.Nb_ligne();
    try {pt = new double [lb*dim];}
    catch (ErrSortieFinale)
         // cas d'une direction voulue vers la sortie
         // on relance l'interuption pour le niveau supérieur
       { ErrSortieFinale toto;
         throw (toto);
       }
    catch (...)
     {cout << "\n Erreur d'allocation de memoire: taille de reel demandee: "
           << lb <<"*"<<dim <<"= "<< (lb*dim) ;
      if (ParaGlob::NiveauImpression() > 2)
          cout << "\n MatBand::MatBand (... ";
      Sortie(1);
     };
    for (double * pti = pt,* pt2i = m.pt; pti < pt+lb*dim; pti++,pt2i++)
      *pti = *pt2i;
  };  

// DESTRUCTEUR :
MatBand::~MatBand ()  
  { if ( pt != NULL) 
         delete [] pt;
   };           

// fonction permettant de creer une nouvelle instance d'element
Mat_abstraite * MatBand::NouvelElement() const
 { Mat_abstraite * a;
   a = new MatBand(*this);
   return a;
  };
        	
		
// surcharge de l'opérateur d'affectation, cas de matrices abstraites
Mat_abstraite & MatBand::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 << "MatBand::operator = ( const Mat_abstraite & b)" << endl;
       Sortie(1);
      }
   #endif

   const MatBand & a = *((MatBand*) & 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 ((lb!=a.lb)||(dim!=a.dim))
      {cout << "erreur dans l'operation d'affectation";
       cout << " les dimensions ne sont pas equivalentes dim =  " << dim << " , "  << a.dim 
            << " largeur de bande = " << lb  << " , " << a.lb << '\n';
       cout << "MatBand::operator = ( const Mat_abstraite & a)" << endl;
       Sortie(1);
      };
   #endif
   int total = lb*dim;
   for (int i=0;i< total; i++)
     pt[i] = a.pt[i];
   return *this;
  };
  
// surcharge de l'opérateur d'affectation, cas de matrices bandes
MatBand & MatBand::operator = ( const MatBand & a)
 { 
   #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 ((lb!=a.lb)||(dim!=a.dim))
      {cout << "erreur dans l'operation d'affectation";
       cout << " les dimensions ne sont pas equivalentes dim =  " << dim << a.dim 
            << " largeur de bande = " << lb << a.lb << '\n';
       cout << "MatBand::operator = ( const MatBand & a)" << endl;
       Sortie(1);
      }
   #endif
   int total = lb*dim;
   for (int i=0;i< total; i++)
     pt[i] = a.pt[i];
   return *this;
  };
        
	// Surcharge de l'operateur += : addition d'une matrice a la matrice courante
void MatBand::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 << "MatBand::operator+= (const Mat_abstraite& b)" << endl;
       Sortie(1);
      }
    // si les tailles ne sont pas identique pb
    if (this->Nb_ligne() != b.Nb_ligne())
      {cout << " les matrices ont un nombre de ligne différent " 
            << this->Nb_ligne() << "  " << b.Nb_ligne() << '\n';
       cout << "MatBand::operator+= (const Mat_abstraite& b)" << endl;
       Sortie(1);
      }
   #endif
   
   if (b.Type_matrice() != DIAGONALE)
    {  // cas de deux matrices bandes
       const MatBand & mat_b = *((MatBand*) & b);
       #ifdef MISE_AU_POINT 
        if ((lb!=mat_b.lb)||(dim!=mat_b.dim))
         {cout << "erreur dans l'operation d'affectation";
          cout << " les dimensions ne sont pas equivalentes dim =  " << dim << mat_b.dim 
               << " largeur de bande = " << lb << mat_b.lb << '\n';
          cout << "MatBand::operator+= (const Mat_abstraite& b)" << endl;
          Sortie(1);
          }
       #endif
       int total = lb*dim;
       for (int i=0;i< total; i++)
          pt[i] += mat_b.pt[i];
      } 
   else // cas d'une matrice argument diagonale
    {  int nb_lign = this->Nb_ligne();
       const MatDiag & mat_d = *((MatDiag*) & b);
	   for (int i=1;i<=nb_lign;i++)
	    (*this)(i,i) += mat_d(i,i);
     }
  };

	// Surcharge de l'operateur -= : soustraction d'une matrice a la matrice courante
void MatBand::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 << "MatBand::operator-= (const Mat_abstraite& b)" << endl;
       Sortie(1);
      }
    // si les tailles ne sont pas identique pb
    if (this->Nb_ligne() != b.Nb_ligne())
      {cout << " les matrices ont un nombre de ligne différent " 
            << this->Nb_ligne() << "  " << b.Nb_ligne() << '\n';
       cout << "MatBand::operator-= (const Mat_abstraite& b)" << endl;
       Sortie(1);
      }
   #endif
   
   if (b.Type_matrice() != DIAGONALE)
    {  // cas de deux matrices bandes
       const MatBand & mat_b = *((MatBand*) & b);
       #ifdef MISE_AU_POINT 
        if ((lb!=mat_b.lb)||(dim!=mat_b.dim))
         {cout << "erreur dans l'operation d'affectation";
          cout << " les dimensions ne sont pas equivalentes dim =  " << dim << mat_b.dim 
               << " largeur de bande = " << lb << mat_b.lb << '\n';
          cout << "MatBand::operator-= (const Mat_abstraite& b)" << endl;
          Sortie(1);
          }
       #endif
       int total = lb*dim;
       for (int i=0;i< total; i++)
          pt[i] -= mat_b.pt[i];
      } 
   else // cas d'une matrice argument diagonale
    {  int nb_lign = this->Nb_ligne();
       const MatDiag & mat_d = *((MatDiag*) & b);
	   for (int i=1;i<=nb_lign;i++)
	    (*this)(i,i) -= mat_d(i,i);
     }
  };
		
	// Surcharge de l'operateur *= : multiplication de la  matrice courante par un scalaire
void MatBand::operator*= (const double r)
 { 
   int total = lb*dim;
   for (int i=0;i< total; i++)
     pt[i] *= r;
  };

	// Surcharge de l'operateur == : test d'egalite entre deux matrices
int MatBand::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 << "int MatBand::operator== (const Mat_abstraite& mat_pl)" << endl;
       Sortie(1);
      }
   #endif
   // si les tailles ne sont pas identique => différentes
   if (this->Nb_ligne() != b.Nb_ligne()) 
       return 0;
   const MatBand & mat_b = *((MatBand*) & b);
   if ((lb!=mat_b.lb)||(dim!=mat_b.dim))
       return 0;
   int total = lb*dim;
   for (int i=0;i< total; i++)
     if ( pt[i] != mat_b.pt[i])
       return 0;
   return 1;
  };
 
 // transfert des informations de *this dans la matrice passée en paramètre
 // la matrice paramètre est au préalable, mise à 0.
 void MatBand::Transfert_vers_mat( Mat_abstraite & b )
{b.Initialise(0.); // init de la matrice
 // puis on transfert
 for (int i=1;i<=dim;i++)
   { for (int j=1;j<= lb; j++)
       b(i,j)=(*this)(i,j);
   };
};

// sortie d'une valeur acces au coordonnees en matrice carree
double& MatBand::operator () (int i, int j )
 { 
   #ifdef MISE_AU_POINT 
	   bool erreur = false;
    if (lb==0)
      {cout << "\n erreur d acces aux composantes d une matrice bande";
       cout << " la largeur de bande est nulle " << '\n';
       cout << "double MatBand::operator () (int i, int j )" << endl;
		     erreur = true;
      };
    if (i<1 || abs(i-j) >= lb )
      {cout << "\n erreur d acces aux composantes d une matrice bande";
       cout << " indice hors borne (" << i << "," << j <<"), lb = " << lb; 
       cout << ", condition inacceptable: i<1 || abs(i-j) >= lb " << '\n';
       cout << "double MatBand::operator () (int i, int j )" << endl;
		     erreur = true;
      };
    if (j<1 || j>dim)
      {cout << "\n erreur d acces aux composantes d une matrice bande";
       cout << " second indice hors borne j = " << j <<" dim = " << dim << '\n';
       cout << "double MatBand::operator () (int i, int j )" << endl;
		     erreur = true;
      };
	   if (erreur)
	     Sortie(1);
   #endif
   if (i>=j)   // acces a la partie inferieure
       return pt[i*(lb-1)+j-1];
   else       // acces theorique a la partie sup qui est egale a la partie inf
       return pt[j*(lb-1)+i-1];     
 };
 
// sortie en lecture seule, d'une valeur acces au coordonnees en matrice carree
double MatBand::operator () (int i, int j ) const 
 { 
   #ifdef MISE_AU_POINT  
    bool erreur = false;
    if (lb==0)
      {cout << "erreur d acces aux composantes d une matrice bande";
       cout << " la largeur de bande est nulle " << '\n';
       cout << "double MatBand::operator () (int i, int j )" << endl;
		     erreur = true;
      };
//    if (i<1 || abs(i-j) >= lb )
    if (i<1 )
      {cout << "erreur d acces aux composantes d une matrice bande";
       cout << " indice hors borne (" << i << "," << j <<"), lb = " << lb; 
       cout << ", condition acceptable: i<1 || abs(i-j) >= lb " << '\n';
       cout << "double MatBand::operator () (int i, int j )" << endl;
 		    erreur = true;
     };
    if (j<1 || j>dim)
      {cout << "erreur d acces aux composantes d une matrice bande";
       cout << " second indice hors borne j = " << j <<" dim = " << dim << '\n';
       cout << "double MatBand::operator () (int i, int j )" << endl;
		     erreur = true;
      };
	   if (erreur)
	     Sortie(1);
   #endif
   // dans le cas où les indices sont hors bornes, on sort 0
   if (abs(i-j) >= lb)
     return 0;
   if (i>=j)   // acces a la partie inferieure
       return pt[i*(lb-1)+j-1];
   else       // acces theorique a la partie sup qui est egale a la partie inf
       return pt[j*(lb-1)+i-1];     
 };
 
// Retourne la ieme ligne de la matrice
Vecteur MatBand::Ligne(int i) const 
  {  Vecteur ligne(dim); // mise a zero de ligne
     int maxJ=MiN(dim,i+lb-1);
     for (int j= MaX(1,i-lb+1);j<= maxJ;j++)
           ligne(j) = (*this)(i,j);
     return ligne;  
   };
 
// Retourne la ieme ligne de la matrice
// sous le format de stokage propre a la matrice
// donc a n'utiliser que comme sauvegarde en parralele
// avec la fonction RemplaceLigne
Vecteur MatBand::LigneSpe(int i) const 
  {  Vecteur ligne(2*lb-1); // mise a zero de ligne
     int le_min = MaX(1,i-lb+1);
     int le_max = MiN(dim,i+lb-1);
     for (int j= le_min;j<= le_max;j++)
           ligne(j-i+lb) = (*this)(i,j);
     return ligne;  
   };
// remplace la ligne de la matrice par la ligne fournie
void MatBand::RemplaceLigneSpe(int i,const Vecteur & v)
  { 
     #ifdef MISE_AU_POINT 
      if (v.Taille() != 2*lb-1)
        {cout << "\nerreur d affectation pour le vecteur";
         cout << " dim vecteur = " << v.Taille() << " au lieu de " << 2*lb-1 <<'\n';
         cout << "void MatBand::RemplaceLigne(int i,Vecteur & v)" << endl;
        }
     #endif
     for (int j= MaX(1,i-lb+1);j<= MiN(dim,i+lb-1);j++)
           (*this)(i,j)=v(j-i+lb);
   };
//met une valeur identique sur toute la ligne
void MatBand::MetValLigne(int i,double x)
   { 
     for (int j= MaX(1,i-lb+1);j<= MiN(dim,i+lb-1);j++)
          (*this)(i,j) = x;
    };
// Retourne la jeme colonne de la matrice
Vecteur MatBand::Colonne(int j) const 
  {  Vecteur col(dim); // mise a zero de colonne
     // bien se rappeler qu'ici lb contient la diagonale
     int maxI = MiN(dim,j+lb-1);
     int le_mini = MaX(1,j-lb+1);
     for (int i= le_mini;i<= maxI;i++)
           col(i) = (*this)(i,j);
     return col;  
   };
		
// 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 MatBand::ColonneSpe(int j) const 
  {  Vecteur col(2*lb-1); // mise a zero de colonne
     for (int i= MaX(1,j-lb+1);i<= MiN(dim,j+lb-1);i++)
           col(i-j+lb) = (*this)(i,j);
     return col;  
   };
// remplace la Colonne de la matrice par la colonne fournie
void MatBand::RemplaceColonneSpe(int j,const Vecteur & v)
  { 
    #ifdef MISE_AU_POINT 
      if (v.Taille() != 2*lb-1)
        {cout << "\nerreur d affectation pour le vecteur";
         cout << " dim vecteur = " << v.Taille() << " au lieu de " << 2*lb-1 << '\n';
         cout << "void MatBand::RemplaceColonne(int j,Vecteur & v)" << endl;
        };
		Sortie (1);
     #endif
     for (int i= MaX(1,j-lb+1);i<= MiN(dim,j+lb-1);i++)
           (*this)(i,j) = v(i-j+lb) ;
   };
//met une valeur identique sur toute la colonne
void MatBand::MetValColonne(int j,double y)
  { 
     for (int i= MaX(1,j-lb+1);i<= MiN(dim,j+lb-1);i++)
         (*this)(i,j)=y;
   };

// Affichage des valeurs de la matrice
// uniquement le valeurs de la bande inferieur
void MatBand::Affiche () const 
  { cout << "\n affichage d une matrice bande de largeur= " << 
         lb << " dimension= " << dim << '\n';
    int i,j;
    for (i=1;i<=dim;i++)
      { for (j=1;j<=i;j++)
          if(abs(i-j)<lb)
             { cout << setw(4) << " i= "  << i << " j= " << j << " *  "
                    << setw(14) << setprecision(7) << (*this)(i,j) << '\n';
             }  
      }
    cout << endl ;                
  };		
// 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
void MatBand::Affiche1(int min_i,int max_i,int pas_i,int min_j,int max_j,int pas_j)
 const 
  { cout << "\n affichage d une matrice bande de largeur= " << 
         lb << " dimension= " << dim << '\n';
    int i,j;
    cout << "           ";
    for (j=min_j;j<=max_j;j+=pas_j) cout << "col  " << setw(4) << j << "  " ;
    cout << '\n';
    for (i=min_i;i<=max_i;i+=pas_i)
      { cout << "lig " << setw(4) << i << " * ";
        for (j=min_j;j<=max_j;j+=pas_j)
          if(abs(i-j)<lb)
             cout << setw(11) << setprecision(7) << (*this)(i,j);
          else
             cout << "...........";   
        cout << '\n';
      }
    cout <<  endl;                
  };   
// 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
void MatBand::Affiche2
   (int min_i,int max_i,int pas_i,int min_j,int max_j,int pas_j,int nd)
   const 
  { cout << "\n affichage d une matrice bande de largeur= " << 
         lb << " dimension= " << dim << '\n';
    if (nd < 7) 
      { cout << " dimension des colonnes trop faible (<7 ! ) " << endl;
        return;
      }       
    int i,j;
    SortBlanc(7+9);

    for (j=min_j;j<=max_j;j+=pas_j) 
      { cout << "col " << setw(2) << j ;
        SortBlanc(nd-6);
       }  
    cout << '\n';
    for (i=min_i;i<=max_i;i+=pas_i)
      { cout << "li " << setw(3) << i << "* ";
        for (j=min_j;j<=max_j;j+=pas_j)
          if(abs(i-j)<=lb-1)
             cout << setw(nd) << setprecision(7) << (*this)(i,j);
          else
             Sortcar(nd,'.');   
        cout << '\n';
      }
    cout <<  endl;                
  };   

void MatBand::Change_taille(int lb1, int dim1 ) // changement de la taille de la matrice 
  { 
    #ifdef MISE_AU_POINT 
     if ((lb < 0) || (dim1 < 0))
       {cout << "erreur de dimensionnement d une matrice bande";
        cout << " LB =  " << lb << " dim = " << dim << '\n';
        cout << "void MatBand::Change_taille(int lb, int dim )" << endl;
      }
    #endif
  
    if ( (lb != lb1) || (dim != dim1))
      { delete [] pt;
        lb = lb1;
        dim = dim1;
        try {pt = new double [lb*dim];}
        catch (ErrSortieFinale)
             // cas d'une direction voulue vers la sortie
             // on relance l'interuption pour le niveau supérieur
           { ErrSortieFinale toto;
             throw (toto);
           }
        catch (...)
         {cout << "\n Erreur d'allocation de memoire: taille de reel demandee: "
               << lb <<"*"<<dim <<"= "<< (lb*dim) ;
          if (ParaGlob::NiveauImpression() > 2)
              cout << "\n MatBand::Change_taille (... ";
          Sortie(1);
         };
        for (double * pti = pt; pti <= pt+lb*dim-1; pti++)
         *pti = 0.; 
       }
  };       
void MatBand::Initialise (double a) // initialisation de la matrice a la valeur "a"
  { for (double * pti = pt; pti <= pt+lb*dim-1; pti++)
         *pti = a;
  }; 
void MatBand::Libere ()   	// Liberation de la place memoire 
  { if ( pt != NULL) 
         delete [] pt;
    lb = 0;
    dim = 0;
    pt = NULL;     
  };
  
// ---------------------------------------------------------------  
//        Resolution du systeme Ax=b
// ---------------------------------------------------------------

			//1) avec en sortie un new vecteur
Vecteur MatBand::Resol_syst 
   (const Vecteur& b,const double &tole, const int maxi, const int rest)
   { // sauvegarde du vecteur b
     Vecteur bb(b);
     if (type_resolution==CHOLESKY)
       ResoBand (*this,lb,dim,bb);
     else
       Mat_abstraite::Resolution_syst(b,bb,tole,maxi,rest);
     return bb;
   };
			 //2) avec en sortie le vecteur d'entree 
Vecteur& MatBand::Resol_systID 
    (Vecteur& b,const double &tole, const int maxi, const int rest)
   { // résolution
     if (type_resolution==CHOLESKY)
       ResoBand (*this,lb,dim,b);
     else
      {Vecteur bb(b);
       Mat_abstraite::Resolution_syst(b,bb,tole,maxi,rest);
       b=bb;
       }
     return b;
   };
		    //3) avec en entrée un tableau de vecteur second membre et
		    //    en sortie un nouveau tableau de vecteurs  
Tableau <Vecteur> MatBand::Resol_syst 
 (const Tableau <Vecteur>& b,const double &tole, const int maxi, const int rest) 
   { // sauvegarde de b
     Tableau <Vecteur> bb(b);
     if (type_resolution==CHOLESKY)
        ResoBand (*this,lb,dim,bb);
     else
       Mat_abstraite::Resolution_syst(b,bb,tole,maxi,rest);
     return bb;
   };
		    //4) avec en entrée un tableau de vecteur second membre et
		    //    en sortie le tableau de vecteurs d'entree 
Tableau <Vecteur>& MatBand::Resol_systID 
 (Tableau <Vecteur>& b,const double &tole, const int maxi, const int rest) 
   { // résolution
     if (type_resolution==CHOLESKY)
        ResoBand (*this,lb,dim,b);
     else
      {Tableau <Vecteur> bb(b);
       Mat_abstraite::Resolution_syst(b,bb,tole,maxi,rest);
       b=bb;
       }
     return b;
   };

// ++++ méthodes  spécifique à la classe
	// Resolution du systeme Ax=b, sans triangulation, c'est-à-dire que l'on
	// considère que la matrice stocke actuellement la triangulation
Vecteur& MatBand::Resol_systID_sans_triangul ( Vecteur& b)
   {
      #ifdef MISE_AU_POINT 
         if (dim != b.Taille())
           {cout << "erreur taille matrice < taille vecteur";
            cout << "dimMat =  " << dim << " dimVect = " << b.Taille() << '\n';
            cout << "MatBand::Resol_systID_sans_triangul ( Vecteur& b)" << endl;
            Sortie(1); 
            }
      #endif
      RESOLT(*this,b,dim,lb);
      return b;
    };
	
//5) avec en sortie le dernier vecteur d'entree, le premier étant le second membre
//   et restant inchangé, en sortie c'est donc soit le retour ou soit vortie, les
//   deux étant identiques
Vecteur& MatBand::Resol_systID_2 (const Vecteur& b,Vecteur& vortie
		                     , const double &tole,const int maxi,const int rest)
   { if (type_resolution==CHOLESKY)
       {vortie = b;
        ResoBand (*this,lb,dim,vortie);
        }
     else
       Mat_abstraite::Resolution_syst(b,vortie,tole,maxi,rest);
     return vortie;
   };
   
// ===== RÉSOLUTION EN DEUX TEMPS ================ :
		//    1) préparation de la matrice donc modification de la matrice éventuellement
		//       par exemple pour les matrices bandes avec cholesky : triangulation
void MatBand::Preparation_resol()
{ // dans le cas de cholesky on triangule
  if (type_resolution==CHOLESKY)
   {
      double PMIN= 1.E31;  // maxi a priori
      double PMAX= 0.;   // mini a priori
      double PABS; // variable intermediaire
      for (int I=1;I<=dim;I++)
         { PABS= Dabs(this->c(lb,I));
           PMIN= MiN(PMIN,PABS);
           PMAX= MaX(PMAX,PABS);
         }
//    controle avant reduction
      if (ParaGlob::NiveauImpression() >= 4 ) 
          cout << "\nMatBand avant triangulation " << "PMIN= " << PMIN << " PMAX= " << PMAX << '\n';
      Vecteur CC(dim);  // vecteur de travail      
      REDUCT(*this,lb,dim,CC);
		
//    controle apres reduction
      PMIN= 1.E31;
      PMAX= 0.;
      for (int I=1;I<=dim;I++)
         { PABS= Dabs(this->c(lb,I));
           PMIN= MiN(PMIN,PABS);
           PMAX= MaX(PMAX,PABS);
         }
      if (ParaGlob::NiveauImpression() >= 4 ) 
       { cout << "MatBand après triangulation " << "PMIN= " << PMIN << " PMAX= " << PMAX;
         cout << endl; }     
    };
};
		//    2) *** résolution sans modification de la matrice DOIT ÊTRE PRÉCÉDÉ DE L'APPEL DE 
		//           Preparation_resol     
		      // a) avec en sortie un new vecteur
Vecteur MatBand::Simple_Resol_syst (const Vecteur& b,const double &tol
		       ,const int maxi,const int rest) const 
   { // sauvegarde du vecteur b
     Vecteur bb(b);
     if (type_resolution==CHOLESKY)
       RESOLT(*this,bb,dim,lb);
     else
       Mat_abstraite::Resolution_syst(b,bb,tol,maxi,rest);
     return bb;
   };
		      // b) avec en sortie le vecteur d'entree 
Vecteur& MatBand::Simple_Resol_systID (Vecteur& b,const double &tol
		       ,const int maxi,const int restart ) const 
   { // résolution
     if (type_resolution==CHOLESKY)
       RESOLT(*this,b,dim,lb);
     else
      {Vecteur bb(b);
       Mat_abstraite::Resolution_syst(b,bb,tol,maxi,restart);
       b=bb;
       }
     return b;
   };
		      // c) avec en entrée un tableau de vecteur second membre et
		      //    en sortie un nouveau tableau de vecteurs  
Tableau <Vecteur> MatBand::Simple_Resol_syst 
		                     (const Tableau <Vecteur>& b,const double &tol
		       ,const int maxi,const int restart ) const 
   { // sauvegarde de b
     Tableau <Vecteur> bb(b);
     if (type_resolution==CHOLESKY)
      { // boucle pour résoudre les n seconds membres
        int nbSM = b.Taille();
        for (int i=1;i<=nbSM;i++)
          RESOLT(*this,bb(i),dim,lb);
       }   
     else
       Mat_abstraite::Resolution_syst(b,bb,tol,maxi,restart);
     return bb;
   };
		      // d) avec en entrée un tableau de vecteur second membre et
		      //    en sortie le tableau de vecteurs d'entree 
Tableau <Vecteur>& MatBand::Simple_Resol_systID 
		                        (Tableau <Vecteur>& b,const double &tol
		       ,const int maxi,const int restart ) const 
   { // résolution
     if (type_resolution==CHOLESKY)
      { // boucle pour résoudre les n seconds membres
        int nbSM = b.Taille();
        for (int i=1;i<=nbSM;i++)
          RESOLT(*this,b(i),dim,lb);
       }   
     else
      {Tableau <Vecteur> bb(b);
       Mat_abstraite::Resolution_syst(b,bb,tol,maxi,restart);
       b=bb;
       }
     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
Vecteur& MatBand::Simple_Resol_systID_2 (const Vecteur& b,Vecteur& vortie
		                                 , const double &tol,const int maxi,const int restart) const 
   { if (type_resolution==CHOLESKY)
      { vortie = b;
        RESOLT(*this,vortie,dim,lb);
       } 
     else
       Mat_abstraite::Resolution_syst(b,vortie,tol,maxi,restart);
     return vortie;
   };
		// ===== FIN RÉSOLUTION EN DEUX TEMPS ================ :
         // enchainement de la resolution avec controle des pivots
void MatBand::ResoBand (MatBand & TRAID,int LB,int NBDL,Vecteur& SG)
   {
      #ifdef MISE_AU_POINT 
         if (dim > SG.Taille())
           {cout << "erreur taille matrice < taille vecteur";
            cout << "dimMat =  " << dim << " dimVect = " << SG.Taille() << '\n';
            cout << "void MatBand::ResoBand "
                 << "(MatBand & TRAID,int LB,int NBDL,Vecteur& SG)" << endl;
           }
      #endif

      double PMIN= 1.E31;  // maxi a priori
      double PMAX= 0.;   // mini a priori
      double PABS; // variable intermediaire
      for (int I=1;I<=NBDL;I++)
         { PABS= Dabs(TRAID.c(LB,I));
           PMIN= MiN(PMIN,PABS);
           PMAX= MaX(PMAX,PABS);
         }
      if (ParaGlob::NiveauImpression() >= 4 ) 
        cout << "\n MatBand resolution " << "PMIN= " << PMIN << " PMAX= " << PMAX ;
      Vecteur CC(NBDL);  // vecteur de travail      
     
     // --- essai de traitement particulier de ré-équilibrage de pivots, basé sur la diagonale
/*     double rapport_maxi = 10000; // maxi entre le pivot max et le pivot considéré
     double coef_mult = 100; // doit être = à la racine de rapport_maxi
     double min_pivot = PMAX/rapport_maxi;
     Tableau <bool> pivot_modifier(NBDL,false); // le tableau des pivots modifiés
     bool regularisation_pivot = false; // indicateur global
     if (PMIN < min_pivot) // sinon tout est ok, ce n'est pas la peine de tout re-regarder
      {regularisation_pivot=true;
       cout << "\n ** regularisation pivot ";
       //numLC: numéro de la ligne ou de la colonne
       for (int numLC=1;numLC<=NBDL;numLC++)
       { double val_diag = Dabs(TRAID.c(LB,numLC));
         if (val_diag < min_pivot)
           {pivot_modifier(numLC)=true; // on mémorise
            // on change la ligne de numéro numLC
            int maxJ=MiN(dim,numLC+lb-1);
            for (int j= MaX(1,numLC-lb+1);j<= maxJ;j++)
               (*this)(numLC,j) *= coef_mult;
            // comme uniquement la moitié de la matrice est stocké, cela signifie que la colonne a également
            // été modifiée, sauf le pivot sur la diagonale qui doit être remultiplié
            (*this)(numLC,numLC) *= coef_mult;  
               
 //           // on change la colonne, l'indice de ligne s'appelle I
 //           int maxI = MiN(dim,numLC+lb-1);
 //           for (int I= MaX(1,numLC-lb+1);I<= maxI;I++)
 //              (*this)(I,numLC) *= coef_mult; 
            // on change le second membre
            SG(numLC) *= coef_mult;   
           };
       };
      }; */
//---------- essai 1 première partie
//      double pmaxAvantReduc = PMAX;
//---------- fin essai		
             
      REDUCT(TRAID,LB,NBDL,CC);
      
//    controle apres reduction
      PMIN= 1.E31;
      PMAX= 0.;
      for (int I=1;I<=NBDL;I++)
         { PABS= Dabs(TRAID.c(LB,I));
           PMIN= MiN(PMIN,PABS);
           PMAX= MaX(PMAX,PABS);
//---------- essai 1 deuxième partie
//      double fact = 0.000001;
//      if (PMAX > fact * pmaxAvantReduc)
//		  {this->c(lb,I) = PMAX = fact * pmaxAvantReduc ; };//fact * pmaxAvantReduc;
//---------- fin essai		
         }
      if (ParaGlob::NiveauImpression() >= 4 ) 
       { cout << "\n MatBand resolution " << "PMIN= " << PMIN << " PMAX= " << PMAX;
         cout << endl;  }

      RESOLT(TRAID,SG,NBDL,LB);
      
/*      // s'il y a eu une régularisation il faut en tenir compte sur le résultat
      if (regularisation_pivot)
       for (int i=1;i<=NBDL;i++)
        if (pivot_modifier(i))
          SG(i) *= coef_mult; */
      
    };
         // idem précédent mais avec plusieurs seconds membres
void MatBand::ResoBand (MatBand & TRAID,int LB,int NBDL,Tableau <Vecteur>& SG)
   {
      int nbSM = SG.Taille();
      #ifdef MISE_AU_POINT 
        for (int i=1;i<=nbSM;i++)
         if (dim > SG(i).Taille())
           {cout << "erreur taille matrice < taille vecteur";
            cout << "dimMat =  " << dim << " dimVect = " << SG(i).Taille() << '\n';
            cout << "void MatBand::ResoBand (pour plusieurs second membres)"
                 << "(MatBand & TRAID,int LB,int NBDL,Tableau <Vecteur>& SG)" << endl;
           }
      #endif

      double PMIN= 1.E31;  // maxi a priori
      double PMAX= 0.;   // mini a priori
      double PABS; // variable intermediaire
      for (int I=1;I<=NBDL;I++)
         { PABS= Dabs(TRAID.c(LB,I));
           PMIN= MiN(PMIN,PABS);
           PMAX= MaX(PMAX,PABS);
         }
      if (ParaGlob::NiveauImpression() >= 4 ) 
         cout << "\n MatBand resolution " << "PMIN= " << PMIN << " PMAX= " << PMAX ;
      Vecteur CC(NBDL);  // vecteur de travail      
      REDUCT(TRAID,LB,NBDL,CC);
      
//    controle apres reduction
      PMIN= 1.E31;
      PMAX= 0.;
      for (int I=1;I<=NBDL;I++)
         { PABS= Dabs(TRAID.c(LB,I));
           PMIN= MiN(PMIN,PABS);
           PMAX= MaX(PMAX,PABS);
         }
      if (ParaGlob::NiveauImpression() >= 4 ) 
        { cout << "\n MatBand resolution " << "PMIN= " << PMIN << " PMAX= " << PMAX;
          cout << endl; }

      // boucle pour résoudre les n seconds membres
      for (int i=1;i<=nbSM;i++)
         RESOLT(TRAID,SG(i),NBDL,LB);
    };
    
			// Multiplication d'un vecteur par une matrice ( (vec)t * A )
Vecteur MatBand::Prod_vec_mat ( const Vecteur& vec) const 
    {
      #ifdef MISE_AU_POINT 
         if (dim != vec.Taille())
           {cout << "erreur taille matrice < taille vecteur";
            cout << "dimMat =  " << dim << " dimVect = " << vec.Taille() << '\n';
            cout << "Vecteur& MatBand::Prod_vec_mat ( Vecteur& vec,Vecteur & res)" << endl;
           }
      #endif

      Vecteur res(dim);
      for (int j=1;j<=dim;j++)
        for (int i= MaX(1,j-lb+1);i<= MiN(dim,j+lb-1);i++)
           res(j) = res (j) + vec(i) * (*this)(i,j) ;
      return res;

     };
    
			// Multiplication d'un vecteur par une matrice ( (vec)t * A )
			// ici on utilise la place du second vecteur pour le résultat
Vecteur& MatBand::Prod_vec_mat ( const Vecteur& vec,Vecteur & res) const 
    {
      #ifdef MISE_AU_POINT 
         if (dim != vec.Taille())
           {cout << "erreur taille matrice < taille vecteur";
            cout << "dimMat =  " << dim << " dimVect = " << vec.Taille() << '\n';
            cout << "Vecteur MatBand::Prod_vec_mat ( Vecteur& vec)" << endl;
           }
      #endif
      res.Change_taille(dim);
      for (int j=1;j<=dim;j++)
       {res(j) = 0.;
        for (int i= MaX(1,j-lb+1);i<= MiN(dim,j+lb-1);i++)
           res(j) += vec(i) * (*this)(i,j) ;
        }   
      return res;

     };

		// Multiplication d'une matrice par un vecteur ( A * vec )
Vecteur MatBand::Prod_mat_vec (const Vecteur& vec) const 
    {
      #ifdef MISE_AU_POINT 
         if (dim != vec.Taille())
           {cout << "erreur taille matrice < taille vecteur";
            cout << "dimMat =  " << dim << " dimVect = " << vec.Taille() << '\n';
            cout << "Vecteur MatBand::Prod_mat_vec (const Vecteur& vec)" << endl;
           }
      #endif

      Vecteur res(dim);
      for (int i=1;i<=dim;i++)
        for (int j= MaX(1,i-lb+1);j<= MiN(dim,i+lb-1);j++)
           res(i) = res (i) + (*this)(i,j) * vec(j);
      return res;
    };

		// Multiplication d'une matrice par un vecteur ( A * vec )
	    // ici on utilise la place du second vecteur pour le résultat
Vecteur& MatBand::Prod_mat_vec (const Vecteur& vec,Vecteur & res) const 
    {
      #ifdef MISE_AU_POINT 
         if (dim != vec.Taille())
           {cout << "erreur taille matrice < taille vecteur";
            cout << "dimMat =  " << dim << " dimVect = " << vec.Taille() << '\n';
            cout << "Vecteur& MatBand::Prod_mat_vec (const Vecteur& vec,Vecteur & res)" << endl;
           }
      #endif

      res.Change_taille(dim);
      for (int i=1;i<=dim;i++)
        { res (i) = 0.;
        for (int j= MaX(1,i-lb+1);j<= MiN(dim,i+lb-1);j++)
           res(i) +=  (*this)(i,j) * vec(j);
        }   
      return res;
    };

// Multiplication d'une ligne iligne de la matrice avec un vecteur de
// dimension = le nombre de colonne de la matrice
double MatBand::Prod_Ligne_vec ( int iligne,const Vecteur& vec) const 
   {
     #ifdef MISE_AU_POINT
      if (dim != vec.Taille())
       { cout << " \n erreur la taille du vecteur = " << vec.Taille() <<" n est pas correcte ";
         cout << " dimMatrice = " << dim;
         cout << " double MatBand::Prod_Ligne_vec ( int iligne,Vecteur& vec)" << endl;
         Sortie (1);
        }
     #endif 
     double res = 0;     
     for (int j=MaX(1,iligne+1-lb);j<=MiN(iligne-1+lb,dim);j++)
		res += (*this)(iligne,j) * vec(j);  
     return res;   
    };
// Multiplication d'un vecteur avec une colonne icol de la matrice 
// dimension = le nombre de ligne de la matrice
double MatBand::Prod_vec_col( int jcol,const Vecteur& vec) const 
   {
     #ifdef MISE_AU_POINT
      if (dim != vec.Taille())
       { cout << " \n erreur la taille du vecteur = " << vec.Taille() <<" n est pas correcte ";
         cout << " dimMatrice = " << dim;
         cout << " double MatBand::Prod_vec_col( int icol,Vecteur& vec)" << endl;
         Sortie (1);
        }
     #endif 
     double res = 0;     
     for (int i=MaX(1,jcol+1-lb);i<=MiN(jcol-1+lb,dim);i++)
		res +=  vec(i) * (*this)(i,jcol) ;  
     return res;   
    };
		
// calcul du produit : (vec_1)^T * A * (vect_2)
double MatBand::vectT_mat_vec(const Vecteur& vec1, const Vecteur& vec2)  const 
{ double resu=0.; // valeur de retour
  #ifdef MISE_AU_POINT 
  if (vec1.Taille() != dim)
      {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= " << dim 
            << "\n MatBand::vectT_mat_vec(const Vecteur& vec1, const Vecteur& vec2)"  << endl;
       Sortie(1);
      }
  if (vec2.Taille() != dim)
      {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= " << dim 
            << "\n MatBand::vectT_mat_vec(const Vecteur& vec1, const Vecteur& vec2)"  << endl;
       Sortie(1);
      }
  #endif

  for (int i=1;i<=dim;i++)
    for (int j= MaX(1,i-lb+1);j<= MiN(dim,i+lb-1);j++)
	    resu += vec1(i) * (*this)(i,j) * vec2(j);
  return resu;	
};

      
//========================================================             
//                     METHODES PROTEGEES
//========================================================             

 void MatBand::REDUCT (MatBand & TRAID,int LB,int NDDL,Vecteur& CC)
//C     reduction de la matrice de raideur avant resolution

//      REAL*8     TRAID(LB,*)
//      DIMENSION CC(NDDL)
  
  {   int I,NI,JJ,KK;
      double C,D;
      
      TRAID.c(LB,1) = 1. / TRAID.c (LB,1);
      for ( NI=2; NI<=NDDL ;NI++) 
       {if (NI <= LB) 
         {CC(1)=TRAID.c(LB-NI+1,NI);
          if (NI> 2)
           { for (JJ=2; JJ<=NI-1; JJ++) 
             { C=TRAID.c(LB-NI+JJ,NI);
               for (KK=1; KK<= JJ-1; KK++)  
                   C=C-CC(KK)*TRAID.c(LB-JJ+KK,JJ);
               CC(JJ)=C;
              }
            }  
          I=1;
          for (JJ=LB-NI+1;JJ<=LB-1; JJ++)
            { TRAID.c(JJ,NI)=CC(I)*TRAID.c(LB,I);
              I=I+1;
             } 
          D=TRAID.c(LB,NI);
          for(KK=1;KK<=NI-1;KK++)
            D=D-CC(KK)*TRAID.c(LB-NI+KK,NI);
         }
        else 
         {      
          CC(1)=TRAID.c(1,NI);
          for (JJ=2;JJ<=LB-1;JJ++)
            { C=TRAID.c(JJ,NI);
              for (KK=1;KK<=JJ-1;KK++)
                C=C-CC(KK)*TRAID.c(KK+LB-JJ,NI-LB+JJ);
              CC(JJ)=C;
             } 
          I=1;
          for (JJ=1;JJ<=LB-1;JJ++)
            { TRAID.c(JJ,NI)=CC(I)*TRAID.c(LB,NI-LB+JJ);
              I=I+1;
             }
          D=TRAID.c(LB,NI);
          for (KK=1;KK<=LB-1;KK++)
            D=D-CC(KK)*TRAID.c(KK,NI);
         }   
       TRAID.c(LB,NI)=1./D;
     }  
 };
  
  
  
  

// ************************* RESOLUTION ****************************

void MatBand::RESOLT(const MatBand & TRAID,Vecteur & TSM,int NDDL,int LB) const 

 {

   //    DIMENSION TRAID(LB,1),TSM(1)
  
//---------------------------execution ---------------------------
  
      int NI,JJ,NII;
      double R;
      for (NI=2;NI<=NDDL;NI++)
        { if (NI<=LB)
           { R=TSM(NI);
             for (JJ=LB-NI+1;JJ<=LB-1;JJ++)
                R=R-TRAID.cConst(JJ,NI)*TSM(NI-LB+JJ);
             TSM(NI)=R;
            }
          else
           {   
             R=TSM(NI);
             for (JJ=1;JJ<=LB-1;JJ++)
               R=R-TRAID.cConst(JJ,NI)*TSM(NI-LB+JJ);
             TSM(NI)=R;
            } 
        }

//
      for( NI=1;NI<=NDDL;NI++)
        TSM(NI)=TSM(NI)*TRAID.cConst(LB,NI);
//
     for (NI=1;NI<=NDDL-1;NI++)
       { NII=NDDL-NI;
         if (NI+1<=LB)
          { R=TSM(NII);
            for (JJ=1;JJ<=NI;JJ++)
              R=R-TRAID.cConst(LB-JJ,NII+JJ)*TSM(NII+JJ);
            TSM(NII)=R;
          }
         else
          {   
            R=TSM(NII);
            for (JJ=1;JJ<=LB-1;JJ++)
              R=R-TRAID.cConst(LB-JJ,NII+JJ)*TSM(NII+JJ);
            TSM(NII)=R;
          }  
        }  
  };


// surcharge de l'operateur de lecture typée
istream & operator >> (istream & entree, MatBand & mat)
  { // vérification du type
    string type;
    entree >> type;
    if (type != "MatBand")
      {Sortie (1);
       return entree;
       }
    //   les dimensions
    entree >> mat.lb >> mat.dim ; 
    // les datas
    for (double * pti = mat.pt; pti <= mat.pt+mat.lb * mat.dim-1; pti++)
      entree >> *pti ;
    return entree;      
  };
  
// surcharge de l'operateur d'ecriture typée
ostream & operator << ( ostream & sort,const  MatBand & mat)
  { // un indicateur donnant le type puis les dimensions
    sort << "MatBand " << mat.lb << " " << mat.dim << " " ; 
    // les datas
    for (double * pti = mat.pt; pti <= mat.pt+mat.lb * mat.dim-1; pti++)
      sort << *pti << " ";
    return sort;      
  };