// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
// AUTHOR : Gérard Rio
// E-MAIL  : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.

/************************************************************************
 *     DATE:        23/01/97                                            *
 *                                                                $     *
 *     AUTEUR:      G RIO   (mailto:gerardrio56@free.fr)                *
 *                                                                $     *
 *     PROJET:      Herezh++                                            *
 *                                                                $     *
 ************************************************************************
 *     BUT: definition de matrices creuses a valeurs reelles.           *
 *     Le stockage est quelconque (pas symétrique par exemple).         *
 *     Ici le stockage est de type en colonne compressée.               *
 *                                                                $     *
 *     ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''     *                                                                      *
 *     VERIFICATION:                                                    *
 *                                                                      *
 *     !  date  !   auteur   !       but                          !     *
 *     ------------------------------------------------------------     *
 *     !        !            !                                    !     *
 *                                                                $     *
 *     ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''     *
 *     MODIFICATIONS:                                                   *
 *     !  date  !   auteur   !       but                          !     *
 *     ------------------------------------------------------------     *
 *                                                                $     *
 ************************************************************************/
#ifndef MAT_CREUSE_COMPCOL_H
#define MAT_CREUSE_COMPCOL_H

//#include "Debug.h"
#include <iostream>
#include "Mat_abstraite.h"
#include "MathUtil.h"

#include "compcol_double.h"
#include "comprow_double.h"
#include "coord_double.h"

#include "Vecteur.h"
#include "Tableau2_T.h"

/// @addtogroup Les_classes_Matrices
///  @{
///


class Mat_creuse_CompCol : public Mat_abstraite  , public CompCol_Mat_double
{   // surcharge de l'operator de lecture typée
    friend istream & operator >> (istream &, Mat_creuse_CompCol &);
    // surcharge de l'operator d'ecriture typée
    friend ostream & operator << (ostream &, const Mat_creuse_CompCol &);
  public :
    // CONSTRUCTEURS :
    Mat_creuse_CompCol ();  // par defaut
    
    // constructeur parmettant la création d'une matrice creuse 
    // ici aucun élément de la matrice n'est créé, il faut ensuite utiliser l'initialisation
    // pour créer les éléments non nulles dans la matrice
    // M : nombre de lignes
    // N : nombre de colonnes
    //   
    Mat_creuse_CompCol(int M, int N);
    
    // constructeur parmettant la création d'une matrice creuse complete
    // à partir de la donnée d'un ensemble de sous matrices carrées
    // les petites matrices sont carrées, ils sont définit par un tableau ligne qui contiend
    // la position d'un terme de la petite matrice dans la grande matrice
    // exemple : pet_mat(k)(i), pet_mat(k)(j) : indique le numéro de ligne  et le numéro de colonne 
    //           d'un élément non nul dans la grande matrice
    // M : nombre de lignes
    // N : nombre de colonnes
    //   
    Mat_creuse_CompCol(int M, int N,const Tableau < Tableau <int> >& petites_matrices); 

    // de copie à partir d'une même instance
    Mat_creuse_CompCol(const Mat_creuse_CompCol &S) ;
    
    // de copie à partir d'une sparse matrice compressée par ligne
 //   Mat_creuse_CompCol(const Mat_creuse_CompRow &R) ;
    // de copie à partir d'une sparse matrice non compressée
 //   Mat_creuse_CompCol(const Mat_creuse_Coord &CO) ;
        
    // DESTRUCTEUR :
    ~Mat_creuse_CompCol () {};
    
		
    // fonction permettant de creer une nouvelle instance d'element du même type
    Mat_abstraite * NouvelElement() const;
        	
    // METHODES PUBLIQUES :
		
        // surcharge de l'opérateur d'affectation : cas de matrices abstraites
    Mat_abstraite & operator = ( const Mat_abstraite &);
 
    // transfert des informations de *this dans la matrice passée en paramètre
    // la matrice paramètre est au préalable, mise à 0.
    void Transfert_vers_mat( Mat_abstraite & A );

        // surcharge de l'opérateur d'affectation : cas de matrices creuses
    Mat_creuse_CompCol & operator = ( const Mat_creuse_CompCol &);
        
    //-----------------------------------------------------------------------------------
    // --- plusieurs fonctions virtuelles qui agissent en général sur la matrice ---------
    //-----------------------------------------------------------------------------------
    // ici l'idée est d'éviter de construire une nouvelle matrice pour des questions
    // d'encombrement, c'est surtout des méthodes utiles pour le stockage de matrice 
    // de raideur. On met le nombre mini de fonction, pour pouvoir faire des expressions.
    // Cependant aucune de ces fonctions n'est en désaccord avec les fonctions membres 
    // qui crée une matrice en résultat (cf. Mat_pleine)
		
    // Surcharge de l'operateur += : addition d'une matrice a la matrice courante
    // cette méthode n'est possible si et seulement les termes existant de mat_pl
    void operator+= (const Mat_abstraite& mat_pl);

    // Surcharge de l'operateur -= : soustraction d'une matrice a la matrice courante
    // cette méthode n'est possible si et seulement les termes existant de mat_pl
    void operator-= (const Mat_abstraite& mat_pl);
 
    // Surcharge de l'operateur *= : multiplication de la  matrice courante par un scalaire
    void operator*= (const double r);
		
    //------------------------------------------------------------------------------------
    // --- fin de plusieurs fonctions virtuelles qui agissent en général sur la matrice --
    //------------------------------------------------------------------------------------
        
    // Surcharge de l'operateur == : test d'egalite entre deux matrices
    int operator== (const Mat_abstraite& mat_pl) const ;
		
    // acces aux composantes comme pour une matrice carree
    // il y a un message en cas de mauvaises valeurs
    double& operator () (int i, int j )
      { return CompCol_Mat_double::set(i-1,j-1);};
		
		
    // 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
    inline bool Existe(int i, int j) const 
	     {return Exista(i,j);};

     
    // acces en lecture seul, aux composantes comme pour une matrice carree
    // il y a un message en cas de mauvais indice
    // si les indices sont possibles mais que l'élément n'existe pas retour 0
    double operator () (int i, int j ) const 
     {return CompCol_Mat_double::operator () (i-1,j-1);};
    
    // Retourne la ieme ligne de la matrice en lecture / écriture
    // pas util dans le cas des matrices creuses donc
    // non implemente
	   Vecteur& Ligne_set (int i)
     { cout <<" fonction non implementee : "
	           << " Vecteur& Mat_creuse_CompCol::operator() (int i) " << endl;
	      Sortie(1);     
	      i = i; // pour ne pas avoir de message a la compilation !!
	      Vecteur* toto = new Vecteur();     	         
	      return *toto;     
     };
	   
    // Retourne la ieme ligne de la matrice uniquement en lecture 
    // pas util dans le cas des matrices creuses  donc
    // non implemente
	   Vecteur operator() (int i) const 
	    { cout  << " fonction non implementee : "
	            << " Vecteur& Mat_creuse_CompCol::operator() (int i) " << endl;
	      i = i; // pour ne pas avoir de message a la compilation !!     	         
	      return Vecteur(0);     
     };
	         
    // Retourne la ieme ligne de la matrice (lent !!!)
		  Vecteur Ligne(int i) const ;

    // Retourne la ieme ligne de la matrice (lent !!!)
    // sous le format de stokage propre a la matrice
    // donc a n'utiliser que comme sauvegarde en parralele
    // avec la fonction RemplaceLigne
    Vecteur LigneSpe(int i) const ;
    // remplace la ligne de la matrice par la ligne fournie (lent !!!)
    void RemplaceLigneSpe(int i,const Vecteur & v);
    
    //met une valeur identique sur toute la ligne (lent !!!)
    void MetValLigne(int i,double val);
    
    // Retourne la jeme colonne de la matrice (rapide)
    Vecteur Colonne(int j) const ;
    
    // Retourne la jeme colonne de la matrice (rapide)
    // sous le format de stokage propre a la matrice
    // donc a n'utiliser que comme sauvegarde en parralele
    // avec la fonction RemplaceColonne
    Vecteur ColonneSpe(int j) const ;
    // remplace la Colonne de la matrice par la ligne fournie (rapide)
    void RemplaceColonneSpe(int j, const Vecteur & v);
    //met une valeur identique sur toute la colonne (rapide)
    void MetValColonne(int j,double val);

    // Affichage des valeurs de la matrice
    void Affiche () const ;
	   // Affiche une partie de la matrice (util pour le debug)
	   // min_ et max_ sont les bornes de la sous_matrice
	   // pas_ indique le pas en i et j pour les indices
	   void Affiche1(int min_i,int max_i,int pas_i,int min_j,int max_j,int pas_j) const ;   
 	  // Affiche une partie de la matrice idem si dessus
 	  // mais avec un nombre de digit (>7) = nd
 	  // si < 7 ne fait rien
	   void Affiche2(int min_i,int max_i,int pas_i,int min_j,int max_j,int pas_j,int nd) const ;   
   
    // changement de la taille de la matrice : c-a-d permet de redimentionner 
    // la matrice creuse
    // 1) cas d'un utilisateur qui connait un ensemble de petite matrice qui doivent être contenue
    // dans la matrice creuse, la matrice initiale est effacée
    // les petites matrices sont carrées, ils sont définit par un tableau ligne qui contiend
    // la position d'un terme de la petite matrice dans la grande matrice
    // exemple : pet_mat(k)(i), pet_mat(k)(j) : indique le numéro de ligne  et le numéro de colonne 
    //           d'un élément non nul dans la grande matrice
    // M et N : nb ligne et nb colonne, 
    void Change_taille(const int M, const int N ,const Tableau < Tableau <int> >& petites_matrices); 
      
    // 2) cas d'un utilisateur qui connait le stockage par colonne
    // M et N : nb ligne et nb colonne, 
    // pointeur_colonne : donne dans le vecteur de stockage global le début de chaque colonne
    // nz     : nb de termes non nulles dans la matrice
    void Change_taille(const int M, const int N ,const  Tableau<int>& pointeur_colonne,const int nz = 0);
     
    // ramène le nombre de colonne
    inline int Nb_colonne ()  const{ return CompCol_Mat_double::dim(1);} ; 
    // ramène le nombre de ligne
    inline int Nb_ligne()  const { return  CompCol_Mat_double::dim(0);};
    // initialisation toutes les valeurs de la matrice a la valeur "a"
    void Initialise (double a); 
    // Liberation de la place memoire
    void Libere () ; 
    // test si la matrice est symétrique ou pas 
    // ici pour l'instant on considère qu'elle ne l'est pas 		
    int Symetrie ()  const { return 0; } 		
	 
			 // Resolution du systeme Ax=b
			 // la verification de taille n'est faite que pour le debug
				//1) avec en sortie un new vecteur
	   Vecteur Resol_syst ( const Vecteur& b,const double &tol = tol_defaut
		      ,const int maxit = maxit_defaut,const int restart = restart_defaut)  ;
			 //2) avec en sortie le vecteur d'entree 
	   Vecteur& Resol_systID ( Vecteur& b,const double &tol = tol_defaut
		      ,const int maxit = maxit_defaut,const int restart = restart_defaut);
    //3) avec en entrée un tableau de vecteur second membre et
    //    en sortie un nouveau tableau de vecteurs  
	   Tableau <Vecteur> Resol_syst 
	          (const Tableau <Vecteur>& b,const double &tol = tol_defaut
		      ,const int maxit = maxit_defaut,const int restart = restart_defaut) ;
    //4) avec en entrée un tableau de vecteur second membre et
    //    en sortie le tableau de vecteurs d'entree
	   Tableau <Vecteur>& Resol_systID 
	          (Tableau <Vecteur>& b,const double &tol = tol_defaut
		      ,const int maxit = maxit_defaut,const int restart = restart_defaut) ;
		  //5) avec en sortie le dernier vecteur d'entree, le premier étant le second membre
		  //   et restant inchangé, en sortie c'est donc soit le retour ou soit vortie, les
		  //   deux étant identiques
	   Vecteur& Resol_systID_2 (const Vecteur& b,Vecteur& vortie
		                     , const double &tol = tol_defaut,const int maxit = maxit_defaut
		                     ,const int restart = restart_defaut);
		                     
	   // ===== RÉSOLUTION EN DEUX TEMPS ================ :
	   //    1) préparation de la matrice donc modification de la matrice éventuellement
	   //       par exemple pour les matrices bandes avec cholesky : triangulation
	   void Preparation_resol();
	   //    2) *** résolution sans modification de la matrice DOIT ÊTRE PRÉCÉDÉ DE L'APPEL DE
	   //           Preparation_resol
	         // a) avec en sortie un new vecteur
	   Vecteur Simple_Resol_syst (const Vecteur& b,const double &tol = tol_defaut
	          ,const int maxit = maxit_defaut,const int restart = restart_defaut) const ;
	         // b) avec en sortie le vecteur d'entree
	   Vecteur& Simple_Resol_systID (Vecteur& b,const double &tol = tol_defaut
	          ,const int maxit = maxit_defaut,const int restart = restart_defaut) const ;
	         // c) avec en entrée un tableau de vecteur second membre et
	         //    en sortie un nouveau tableau de vecteurs
	   Tableau <Vecteur> Simple_Resol_syst
	                        (const Tableau <Vecteur>& b,const double &tol = tol_defaut
	          ,const int maxit = maxit_defaut,const int restart = restart_defaut) const ;
	         // d) avec en entrée un tableau de vecteur second membre et
	         //    en sortie le tableau de vecteurs d'entree
	   Tableau <Vecteur>& Simple_Resol_systID
	                           (Tableau <Vecteur>& b,const double &tol = tol_defaut
	          ,const int maxit = maxit_defaut,const int restart = restart_defaut) const ;
	         // e) avec en sortie le dernier vecteur d'entree, le premier étant le second membre
	         //   et restant inchangé, en sortie c'est donc soit le retour ou soit vortie, les
	         //   deux étant identiques
	   Vecteur& Simple_Resol_systID_2 (const Vecteur& b,Vecteur& vortie
	                                    , const double &tol = tol_defaut
	                                    ,const int maxit = maxit_defaut
	                                    ,const int restart = restart_defaut) const ;
	   // ===== FIN RÉSOLUTION EN DEUX TEMPS ================ :
	   
	   // Multiplication d'un vecteur par une matrice ( (vec)t * A )
    Vecteur Prod_vec_mat ( const  Vecteur& vec)  const ;
	   // idem mais on utilise la place du second vecteur pour le résultat
	   Vecteur& Prod_vec_mat ( const Vecteur& vec, Vecteur & resul) const ;
	   // Multiplication d'une matrice par un vecteur ( A * vec )
    Vecteur Prod_mat_vec ( const  Vecteur& vec) const ;
	   // idem mais on utilise la place du second vecteur pour le résultat
    Vecteur& Prod_mat_vec ( const Vecteur& vec, Vecteur & resul) const;
    // Multiplication d'une ligne iligne de la matrice avec un vecteur de
    // dimension = le nombre de colonne de la matrice
    double Prod_Ligne_vec ( int iligne, const Vecteur& vec) const ;
    // Multiplication d'un vecteur avec une colonne icol de la matrice 
    // dimension = le nombre de ligne de la matrice
    double Prod_vec_col( int icol, const Vecteur& vec) const ; 
		
	   // calcul du produit : (vec_1)^T * A * (vect_2)
	   double vectT_mat_vec(const Vecteur& vec1, const Vecteur& vec2) const ;

    // retourne la place que prend la matrice en entier
    int Place() const { return (NumNonzeros() * 3 + dim(1) + 5);};

//  ====== quelques méthodes spécifiques à la classe MV_Vector<double>==============
    // définie la surcharge de multiplication d'une matrice par un MV_Vector
    // ici de type constant, ramène un nouveau MV_Vector<double>
    inline virtual MV_Vector<double> operator * (const MV_Vector<double> & vec) const
      { return this->CompCol_Mat_double::operator *(vec); };
    // multiplication matrice transposée fois vecteur ce qui est équivalent à 
    // la transposée du résultat de : multiplication vecteur transposé fois matrice
    inline virtual MV_Vector<double> trans_mult(const MV_Vector<double> &x) const 
      { return this->CompCol_Mat_double::trans_mult(x); };  
//      ====== fin des quelques méthodes spécifiques à la classe MV_Vector<double>=======
       
  protected :
    
    // VARIABLES PROTEGEES :
 };
 /// @}  // end of group

#endif