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

#ifndef  Algo_zero_deja_inclus
  #include "Algo_zero.h"
  #include "ParaGlob.h"
#endif

// contient que les méthodes nécessitant les templates, qui doivent donc être
//incluses dans le .h

#ifndef  Algo_zero_deja_inclus
  #include "MathUtil.h"
  #include "ConstMath.h"
  #include "ParaGlob.h"


    //recherche du zéro d'une fonction en utilisant la méthode de Newton-Raphson incrémentale
    // ici il s'agit d'une fonction à une variable
    //  *Ptfonc : le pointeur de la fonction dont il faut chercher le zéro
    //  *Ptder_fonc : pointeur de la dérivée de la fonction
    //               pour ces deux fonctions, la variable alpha est un facteur de charge qui varie de 
    //               0 à 1. Lorsque alpha=0, la racine vaut val_initiale, et ce que l'on 
    //               cherche c'est la racine pour alpha = 1. Lorsque alpha varie progressivement
    //               de 0 à 1, la racine est sensée varier progressivement de val_initiale à résidu
    //               l'utilisation d'alpha permet de faire du Newton incrémentale. 
    //   pour ces deux fonctions, l'argument test ramène
    //     . 1 si le calcul a été ok, -1 s'il y a eu un pb, mais on peut continuer, 0 s'il y a eu un pb
    //       fatal, qui invalide le calcul de la fonction et/ou de la dérivée
    //  val_initiale : une valeur initiale de x pour la recherche de zéro
    //  max_delta_x : si > 0 donne le norme maxi permise sur delta_x au cours d'une itération
    //               si ||delta_x|| est plus grand on fait delta_x = max_delta_x * delta_x / || delta_x||
    //  ramène en sortie:
    //     un booléen qui indique si la résolution est correcte ou non
    //     racine : la racine trouvée
    //     der_at_racine : contient en retour la valeur de la dérivée pour la racine trouvée
    //     nb_incr_total : le nombre total d'incrément qui a été nécessaire
    //     nb_iter_total : le nombre total d'itération, qui cumule les iter de tous les incréments
template <class T>
bool  Algo_zero::Newton_raphson
              (T& instance,double (T::*Ptfonc) (double & alpha,double & x,int& test)
               ,double (T::*Ptder_fonc) (double & alpha,double & x,int& test)
               ,double val_initiale,double & racine,double & der_at_racine
               ,int& nb_incr_total, int& nb_iter_total,double max_delta_x) const
    { // initialisation des différentes variables intermédiaires
      nb_iter_total = nb_incr_total = 0;
      int iter = 0; // nombre d'itérations 
      int incre = 1;  // nombre d'incrément pour converger, au débu = 1
      double x = val_initiale; // initialisation de la valeur courante de x
      int phase_convergence = 0; // a priori on est en phase de convergence
          //  = -3 delta_x trop grand ;   = -2 delta_x trop petit ; = -1 trop d'itération
          //  = -4 : erreur fatale (ou inconnue) sur le calcul du résidu ou de la dérivée
          //  = 0  non convergence normal; = 1 convergence sur un incrément
          //  = 2  convergence finale
      double alpha = 1; // initialisation du facteur de charge
      double delta_alpha = 1; // incrément d'alpha, intit à 1
      double alpha_prec = 0; 
      double x_prec = x; // valeur précédente en début d'incrément
      double delta_alpha_prec = 1; // le delta alpha de l'incrément précédent
      double delta_x_prec = 0.; // le delta x de l'incrément précédent 
      int test =1; // init par défaut
      int glob_err_res = 0;// pour la gestion de test == -1
      int glob_err_der = 0;// pour la gestion de test == -1
      double dabs_residu=ConstMath::tresgrand; // init
      double delta_x_sauve=ConstMath::tresgrand; // ne sert que pour la sortie des messages d'erreur
        // valeur énorme volontairement

      do // boucle globale
       {      
        // calcul de la fonction et test de convergence
        double residu = (instance.*Ptfonc) (alpha,x,test);
        if (  (test == 1) // c-a-d que le calcul est correct
            ||(test == -1) // petit pb mais on continue
           )
         {
          #ifdef MISE_AU_POINT
           if ((!isfinite(residu)) || (isnan(residu)))
            if (((permet_affichage==0) && (ParaGlob::NiveauImpression() >1))|| (permet_affichage > 1))
              { cout << "\n Algo_zero::Newton_raphson: iteration: " << iter  ;
                cout << "\n *** pb resolution residu " << residu << " infini " ;
                cout << "\n Algo_zero::Newton_raphson( " << endl;
                throw ErrNonConvergence_Newton();
              };
          #endif
          if (test == -1) {glob_err_res++;} else {glob_err_res=0;}; // incrémentation ou remise à 0
          // -- manip pour le calcul des précisions
          dabs_residu = Dabs(residu);
          double prec_res_reelle = prec_res_abs + dabs_residu * prec_res_rel;
  //        if (prec_res < 0) // si l'indication est négative , cela veut dire que l'on veut une grandeur négative
  //          prec_res_reelle = MaX(-prec_res * dabs_residu,-prec_res);
          // --
          if (dabs_residu > prec_res_reelle)
            { // cas où l'on n'a pas encore convergé
              der_at_racine = (instance.*Ptder_fonc) (alpha,x,test); // calcul de la dérivée de la fonction
              if (test == -1) {glob_err_der++;} else {glob_err_der=0;}; // incrémentation ou remise à 0
              if ( (test == 1) // cas sans pb
                  || ((test == -1) && ((glob_err_res+glob_err_der) < maxi_test_moins1)) // test non fatal et on est < à maxi_test_moins1 séries de test mauvais
                 )                                   // en supposant que l'on a successivement le résidu et la dérivée en
                                                     // mauvais calcul, sinon ce sera l'un ou l'autre
                { // dans le cas où la dérivée est trop petite on limite
                  if (Dabs(der_at_racine) >= ConstMath::trespetit)
                       der_at_racine = Signe(der_at_racine,ConstMath::trespetit);
        //             {if (der_at_racine >= 0.) der_at_racine =  ConstMath::trespetit;
        //              else der_at_racine = - ConstMath::trespetit;};
                  // détermination du delta_x
                  double delta_x = delta_x_sauve = - residu/der_at_racine;
                  // on applique éventuellement une limitation
                  double val_abs_delta_x = Dabs(delta_x);
                  if ((max_delta_x > ConstMath::trespetit)&&(val_abs_delta_x > max_delta_x))
                   {double limitation = max_delta_x / val_abs_delta_x;
                    delta_x = delta_x_sauve = delta_x * limitation;
                    #ifdef MISE_AU_POINT
                     if (((permet_affichage==0) && (ParaGlob::NiveauImpression() > 4)) || (permet_affichage > 2))
                       { cout << "\n Algo_zero::Newton_raphson: limitation de delta_x: " << limitation;
                       };
                    #endif
                   };
                  
                  x = x + delta_x; // nouvelle valeur de x
                  iter++; // incrémentation du nombre d'itération
// non je ne comprend pas la partie que je commente car le résidu ne change pas donc le test
// sur la convergence ne sera jamais valide, donc cela ne sert à rien de ré-étudier la convergence
//                  // on étudie de nouveau la convergence
//                  dabs_residu = Dabs(residu);
//                  prec_res_reelle = prec_res_abs + dabs_residu * prec_res_rel;
//      //            if (prec_res < 0) // si l'indication est négative , cela veut dire que l'on veut une grandeur négative
//      //               prec_res_reelle = MaX(-prec_res * dabs_residu,-prec_res);
// ???              if (dabs_residu < prec_res_reelle)
//                    { phase_convergence = 1;
//                      // on sauvegarde les variables
//                      delta_x_prec = x - x_prec;
//                      x_prec = x;
//                      alpha_prec = alpha;
//                      delta_alpha_prec = delta_alpha;
//                    }
//                  else
                  // == étude des divergences éventuelles
                   {if (iter > iter_max) {phase_convergence=-1;} // trop d'itération
                    else if (Dabs(delta_x) < mini_delta_x)
                        { phase_convergence=-2; } // delta_x est rendu trop petit
                    else if (Dabs(delta_x) < coef_mini_delta_x * (MaX(Dabs(x-delta_x),Dabs(x))) )                  { phase_convergence=-4; } // delta_x est rendu trop petit
                    else if (Dabs(delta_x) > maxi_delta_x)
                      { phase_convergence = -3;} // delta_x est rendu trop grand
                    else // sinon c'est une non convergence normale, et on continue à itérer
                     { phase_convergence = 0;};
                   };
                 }
              else if ( (!test) || ((glob_err_res+glob_err_der+1) > maxi_test_moins1)) // si test = 0 -> !test= true
                {// que ce soit une erreur fatale ou une suite de mauvais calcul: on signale une divergence
                 phase_convergence = -4;
                }
              else // c'est une erreur non répertoriée
               {cout << "\n %% Algo_zero::Newton_raphson: cas non prevu de retour d'erreur sur la derivee =  "
                     << test << " %% arret algo Newton %% " << flush;
                Sortie(1);
               };
             }
          else // cas ou on a convergé
            { phase_convergence = 1;
              // on sauvegarde les variables
              delta_x_prec = x - x_prec;
              x_prec = x;
              alpha_prec = alpha;
              delta_alpha_prec = delta_alpha;
              #ifdef MISE_AU_POINT
               if (((permet_affichage==0) && (ParaGlob::NiveauImpression() > 6)) || (permet_affichage > 6))
                 { cout << "\n Algo_zero::Newton_raphson: convergence pour " << iter
                        << " iteration(s) " << incre << " increment(s) "
                        << " , alpha=  " << alpha << " et delta_alpha= " << delta_alpha;
                   cout << residu << " residuresidu "<< endl;
                 };
              #endif
             };
         }
        else if (!test)
         {// on signale que l'on est en divergence due à une erreur fatale sur le résidu
          phase_convergence = -4;
          if (x == val_initiale) // cela veut dire que dès le début il y a une erreur fatale
           {if ((permet_affichage > 1) || (ParaGlob::NiveauImpression() > 3))
              cout << "\n %% Algo_zero::Newton_raphson: des le debut, erreur fatale sur le calcul du residu "
                   << " %% arret algo Newton %% " << endl;
            break; // on ne peut pas continuer
           };
         }
        else // c'est une erreur non répertoriée
         {cout << "\n %% Algo_zero::Newton_raphson: cas non prevu de retour d'erreur sur le residu =  "
               << test << " %% arret algo Newton %% " << flush;
          Sortie(1);
         };

        // 2 ----- action en fonction de la convergence ou non
        // -- cas d'une non convergence particulière c-a-d en dehors du cas  phase_convergence == 0  
        // dans le cas où on est en divergence on cherche à augmenter le nombre d'incréments
        if (phase_convergence < 0)
          { if (incre <= nbMaxiIncre)
              {
                #ifdef MISE_AU_POINT
                 if (((permet_affichage==0) && (ParaGlob::NiveauImpression() > 5))|| (permet_affichage > 5))
                   { cout << "\n Algo_zero::Newton_raphson: non convergence pour alpha=  " << alpha
                          << " et delta_alpha= " << delta_alpha
                          << " ( incre= "<<incre << " iter= "<<iter
                          << " ||residu||= "<<dabs_residu
                          << ") " ;
                     if (((permet_affichage==0) && (ParaGlob::NiveauImpression() > 6))|| (permet_affichage > 6))
                      cout << "\n phase_convergence= "<< phase_convergence
                           << " Dabs(delta_x)= " << Dabs(delta_x_sauve)
                           << "(mini_delta_x= "<<mini_delta_x<<" maxi_delta_x= "<<maxi_delta_x
                           << " iter_max= "<<iter_max << ") "
                           << " prec_res_reelle=" << (prec_res_abs + dabs_residu * prec_res_rel);
                      cout << flush;
                   };
                #endif
                // on diminue l'incrément de charge par deux
                delta_alpha *= 0.5;
                nb_iter_total += iter; // stockage du nombre d'iter
                iter = 0; // on reinitialise le nombre d'itération
                // on repart avec une valeur d'alpha plus petite
                alpha = alpha_prec + delta_alpha;
                x = x_prec;
                // on incrémente le compteur d'incréments
                incre++; 
               }
            else // on a dépassé le nombre d'incrément permis, on arrête
              {
                #ifdef MISE_AU_POINT
                 if (((permet_affichage==0) && (ParaGlob::NiveauImpression() > 6))|| (permet_affichage > 6))
                   { cout << "\n %% Algo_zero::Newton_raphson: on atteind le nombre d'increment maxi  "
                          << nbMaxiIncre << " %% arret algo Newton %% " << endl;
                   };
                #endif
                break;
              };
           }
        else if (phase_convergence == 1)
         { // on vient de converger, on regarde si l'on a terminé les incréments de charge
           if (alpha >= 1.) 
             { phase_convergence = 2;} // on a terminé
           else
             { // sinon il faut incrémenter alpha
               alpha += delta_alpha; if (alpha > 1.) alpha = 1.; // pb d'erreur d'arrondi
               // on initialise x en tenant compte du delta_x de l'incrément précédent (dons on extrapole)
               // cela peut-être une bonne idée comme une mauvaise, car cela peut donner
               // des valeurs abérantes, mais normalement on ne fera qu'une passe car on aura une valeur de test
               // de retour des résidus et dérivées en 0 donc avec une erreur fatale
               if (Dabs(delta_x_prec) >= ConstMath::trespetit) 
                   x += delta_x_prec * delta_alpha / delta_alpha_prec;
              }   
          };     
        // -- dans le cas où phase_convergence==0 c-a-d non convergence normal, on refait une boucle
        //    ce qui a pour conséquence d'augmenter iter, et de préciser la racine, c'est tout       
        } while (phase_convergence < 2) ;
   
     racine = x;
     // mise à jour de l'indicateur du nombre d'incréments
     nb_incr_total = incre;
     nb_iter_total += iter; // idem pour les itérations
     // retour de la validitée de la recherche du zéro
     if (phase_convergence == 2)
        return true;
     else
        return false;   
     }
 
    // idem précédemment, mais pour une fonction à valeur vectorielle
    // les matrices sont telles que: der_at_racine(i,j) = df(i)/d(x(j), soit: df = der_at_racine * dx
    // recherche du zéro d'une fonction en utilisant la méthode de Newton-Raphson incrémentale
    // ici il s'agit d'une fonction à n variables
    //  *Ptfonc : le pointeur de la fonction dont il faut chercher le zéro
    //  *Ptder_fonc : pointeur de la dérivée de la fonction
    //               pour ces deux fonctions, la variable alpha est un facteur de charge qui varie de 
    //               0 à 1. Lorsque alpha=0, la racine vaut val_initiale, et ce que l'on 
    //               cherche c'est la racine pour alpha = 1. Lorsque alpha varie progressivement
    //               de 0 à 1, la racine est sensée varier progressivement de val_initiale à résidu
    //               l'utilisation d'alpha permet de faire du Newton incrémentale. 
    //   pour ces deux fonctions, l'argument test ramène
    //     . 1 si le calcul a été ok, -1 s'il y a eu un pb, mais on peut continuer, 0 s'il y a eu un pb
    //       fatal, qui invalide le calcul de la fonction et/ou de la dérivée
    //  val_initiale : une valeur initiale de x pour la recherche de zéro
    //  max_delta_x : si > 0 donne le norme maxi permise sur delta_x au cours d'une itération
    //               si ||delta_x|| est plus grand on fait delta_x = max_delta_x * delta_x / || delta_x||
    //  ramène en sortie:
    //     un booléen qui indique si la résolution est correcte ou non
    //     racine : la racine trouvée
    //     der_at_racine : contient en retour la valeur de la dérivée pour la racine trouvée
    //     nb_incr_total : le nombre total d'incrément qui a été nécessaire
    //     nb_iter_total : le nombre total d'itération, qui cumule les iter de tous les incréments
 template <class T>
 bool  Algo_zero::Newton_raphson
              (T& instance,Vecteur& (T::*Ptfonc) (const double & alpha,const Vecteur & x,int& test)
               ,Mat_abstraite& (T::*Ptder_fonc) (const double & alpha,const Vecteur & x,int& test)
               ,const Vecteur& val_initiale,Vecteur & racine,Mat_abstraite & der_at_racine
               ,int& nb_incr_total, int& nb_iter_total,double max_delta_x) const
    { int dime = val_initiale.Taille();    
      // initialisation des différentes variables intermédiaires
      nb_iter_total = nb_incr_total = 0;
      int iter = 0; // nombre d'itérations
      int incre = 1;  // nombre d'incrément pour converger, au débu = 1
      Vecteur x = val_initiale; // initialisation de la valeur courante de x
      int phase_convergence = 0; // a priori on est en phase de convergence
          //  = -3 delta_x trop grand ;   = -2 delta_x trop petit ; = -1 trop d'itération
          //  = 0  non convergence normal; = 1 convergence sur un incrément
          //  = 2  convergence finale
      double alpha = 1; // initialisation du facteur de charge
      double delta_alpha = 1; // incrément d'alpha, intit à 1
      double alpha_prec = 0; 
      Vecteur x_prec = x; // valeur précédente en début d'incrément
      double delta_alpha_prec = 1; // le delta alpha de l'incrément précédent
      Vecteur delta_x_prec(dime); // le delta x de l'incrément précédent 
      // grandeurs de calcul
      Vecteur residu(dime),delta_x(dime),xmoinsdeltax(dime);
      int test =1; // init par défaut
      int glob_err_res = 0;// pour la gestion de test == -1
      int glob_err_der = 0;// pour la gestion de test == -1
      double residu_Max_val_abs=ConstMath::tresgrand; // init
//      double delta_x_sauve=ConstMath::tresgrand; // ne sert que pour la sortie des messages d'erreur
        // valeur énorme volontairement

      do // boucle globale
       {      
        // calcul de la fonction et test de convergence
        residu = (instance.*Ptfonc) (alpha,x,test);
        if (  (test == 1) // c-a-d que le calcul est correct
            ||(test == -1) // petit pb mais on continue
           )
         {
          residu_Max_val_abs = residu.Max_val_abs();
          #ifdef MISE_AU_POINT
           if ((!isfinite(residu_Max_val_abs)) || (isnan(residu_Max_val_abs)))
            if (((permet_affichage==0) && (ParaGlob::NiveauImpression() >1))|| (permet_affichage > 1))
              { cout << "\n Algo_zero::Newton_raphson: iteration: " << iter  ;
                cout << "\n residu: "; residu.Affiche();
                cout << "\n *** pb resolution residu infini " ;
                cout << "\n Algo_zero::Newton_raphson( " << flush;
                throw ErrNonConvergence_Newton();
              };
          #endif
          if (test == -1) {glob_err_res++;} else {glob_err_res=0;}; // incrémentation ou remise à 0
          // -- manip sur les précisions
          // si prec_res est négatif, il s'agit d'une précision relative
  //        double prec_res_reelle = prec_res;
          double prec_res_reelle = prec_res_abs + residu_Max_val_abs * prec_res_rel;
  //        if (prec_res < 0) // si l'indication est négative , cela veut dire que l'on veut une grandeur négative
  //          prec_res_reelle = MaX(-prec_res * residu_Max_val_abs,-prec_res);
          if (residu_Max_val_abs > prec_res_reelle)
            { // cas où l'on n'a pas encore convergé
              der_at_racine = (instance.*Ptder_fonc) (alpha,x,test); // calcul de la dérivée de la fonction
              #ifdef MISE_AU_POINT
               if (((permet_affichage==0) && (ParaGlob::NiveauImpression() >=8))|| (permet_affichage >= 8))
                 { cout << "\n Algo_zero::Newton_raphson: iteration: " << iter  ;
                   cout << "\n residu: "; residu.Affiche();
                   cout << "\n matrice tangente: " ; der_at_racine.Affiche();
                   cout << flush;
                 };
              #endif
              if (test == -1) {glob_err_der++;} else {glob_err_der=0;}; // incrémentation ou remise à 0
              if ( (test == 1) // cas sans pb
                  || ((test == -1) && ((glob_err_res+glob_err_der) < maxi_test_moins1)) // test non fatal et on est < à maxi_test_moins1 séries de test mauvais
                 )                                   // en supposant que l'on a successivement le résidu et la dérivée en
                                                     // mauvais calcul, sinon ce sera l'un ou l'autre
                { delta_x = -der_at_racine.Resol_systID_2 (residu,delta_x);
                  // on applique éventuellement une limitation
                  double val_abs_delta_x = delta_x.Norme();
                  if ((max_delta_x > ConstMath::trespetit)&&(val_abs_delta_x > max_delta_x))
                   {double limitation = max_delta_x / val_abs_delta_x;
                    delta_x *=  limitation;
                    #ifdef MISE_AU_POINT
                     if (((permet_affichage==0) && (ParaGlob::NiveauImpression() > 4)) || (permet_affichage > 2))
                       { cout << "\n Algo_zero::Newton_raphson: limitation de delta_x: " << limitation;
                       };
                    #endif
                   };
                  
                  x = x + delta_x; // nouvelle valeur de x
                  #ifdef MISE_AU_POINT
                   if (((permet_affichage==0) && (ParaGlob::NiveauImpression() > 7))|| (permet_affichage > 7))
                     { cout << "\n  delta_X :"; delta_x.Affiche();
                       cout << "\n x :"  ;x.Affiche();
                       cout << flush ;
                     };
                  #endif
                  iter++; // incrémentation du nombre d'itération
// non je ne comprend pas la partie que je commente car le résidu ne change pas donc le test
// sur la convergence ne sera jamais valide, donc cela ne sert à rien de ré-étudier la convergence
//              // on étudie de nouveau la convergence
//              residu_Max_val_abs = residu.Max_val_abs();
//              prec_res_reelle = prec_res_abs + residu_Max_val_abs * prec_res_rel;
//  //            if (prec_res < 0) // si l'indication est négative , cela veut dire que l'on veut une grandeur négative
//  //               prec_res_reelle = MaX(-prec_res * residu_Max_val_abs,-prec_res);
//              if (residu_Max_val_abs < prec_res_reelle)
//                { phase_convergence = 1;
//                  // on sauvegarde les variables
//                  delta_x_prec = x - x_prec;
//                  x_prec = x;
//                  alpha_prec = alpha;
//                  delta_alpha_prec = delta_alpha;
//                }
//              else
                 // == étude des divergences éventuelles
                  {xmoinsdeltax = x-delta_x;
                   if (iter > iter_max) {phase_convergence=-1;} // trop d'itération
                   else if (delta_x.Max_val_abs() < mini_delta_x)
                       { phase_convergence=-2; } // delta_x est rendu trop petit
                   else if ((delta_x.Max_val_abs()
                             < coef_mini_delta_x * (MaX( xmoinsdeltax.Max_val_abs(),x.Max_val_abs())) ) )
                       { phase_convergence=-4; } // delta_x est rendu trop petit
                   else if (delta_x.Max_val_abs() > maxi_delta_x)
                     { phase_convergence = -3;} // delta_x est rendu trop grand
                   else // sinon c'est une non convergence normale, et on continue à itérer
                    { phase_convergence = 0;};
                  };
                }
              else if ( (!test) || ((glob_err_res+glob_err_der+1) > maxi_test_moins1)) // si test ==0 -> !test = true
                {// que ce soit une erreur fatale ou une suite de mauvais calcul: on signale une divergence
                 phase_convergence = -4;
                }
              else // c'est une erreur non répertoriée
               {cout << "\n %% Algo_zero::Newton_raphson: cas non prevu de retour d'erreur sur la derivee =  "
                     << test << " %% arret algo Newton %% " << flush;
                Sortie(1);
               };
             }
          else // cas ou on a convergé
            { phase_convergence = 1;
              // on sauvegarde les variables
              delta_x_prec = x - x_prec;
              x_prec = x;
              alpha_prec = alpha;
              delta_alpha_prec = delta_alpha;
              #ifdef MISE_AU_POINT
               if (((permet_affichage==0) && (ParaGlob::NiveauImpression() > 6))|| (permet_affichage > 6))
                 { cout << "\n Algo_zero::Newton_raphson: convergence pour " << iter
                        << " iteration(s) " << incre << " increment(s) "
                        << " , alpha=  " << alpha << " et delta_alpha= " << delta_alpha;
                   residu.Affiche();cout << " residu "<< flush;
                 };
              #endif
             };
         }
        else if (!test)
          {// on signale que l'on est en divergence due à une erreur fatale sur le résidu
           phase_convergence = -4;
           if (x == val_initiale) // cela veut dire que dès le début il y a une erreur fatale
            {if ((permet_affichage > 1) || (ParaGlob::NiveauImpression() > 3))
               cout << "\n %% Algo_zero::Newton_raphson: des le debut, erreur fatale sur le calcul du residu "
                    << " %% arret algo Newton %% " << endl;
             break; // on ne peut pas continuer
            };
          }
        else // c'est une erreur non répertoriée
         {cout << "\n %% Algo_zero::Newton_raphson: cas non prevu de retour d'erreur sur le residu =  "
               << test << " %% arret algo Newton %% " << flush;
          Sortie(1);
         };

        // 2 ----- action en fonction de la convergence ou non
        // -- cas d'une non convergence particulière c-a-d en dehors du cas  phase_convergence == 0  
        // dans le cas où on est en divergence on cherche à augmenter le nombre d'incréments
        if (phase_convergence < 0)
          { if (incre <= nbMaxiIncre)
              {
                #ifdef MISE_AU_POINT
                 if (((permet_affichage==0) && (ParaGlob::NiveauImpression() > 5))|| (permet_affichage > 5))
                   { cout << "\n Algo_zero::Newton_raphson: non convergence pour alpha=  " << alpha
                          << " et delta_alpha= " << delta_alpha
                          << " ( incre= "<<incre << " iter= "<<iter
                          << " ||residu||= "<<residu_Max_val_abs
                          << ") ";
                     if (((permet_affichage==0) && (ParaGlob::NiveauImpression() > 6))|| (permet_affichage > 6))
                      cout << "\n phase_convergence= "<< phase_convergence
                           << " delta_x.Max_val_abs()= " << delta_x.Max_val_abs()
                           << "(mini_delta_x= "<<mini_delta_x<<" maxi_delta_x= "<<maxi_delta_x
                           << " iter_max= "<<iter_max << ") "
                           << " prec_res_reelle=" << (prec_res_abs + residu_Max_val_abs * prec_res_rel);
                     cout << flush;
                   };
                #endif
                // on diminue l'incrément de charge par deux
                delta_alpha *= 0.5;
                nb_iter_total += iter; // stockage du nombre d'iter
                iter = 0; // on reinitialise le nombre d'itération
                // on repart avec une valeur d'alpha plus petite
                alpha = alpha_prec + delta_alpha;
                x = x_prec;
                // on incrémente le compteur d'incréments
                incre++; 
               }
            else // on a dépassé le nombre d'incrément permis, on arrête
              {
                #ifdef MISE_AU_POINT
                if (((permet_affichage==0) && (ParaGlob::NiveauImpression() > 6))|| (permet_affichage > 6))
                   { cout << "\n %% Algo_zero::Newton_raphson: on atteind le nombre d'increment maxi  "
                          << nbMaxiIncre << " %% arret algo Newton %% " << flush;
                   };
                #endif
                break;
              };
           }
        else if (phase_convergence == 1)
         { // on vient de converger, on regarde si l'on a terminé les incréments de charge
           if (alpha >= 1.) 
             { phase_convergence = 2;} // on a terminé
           else
             { // sinon il faut incrémenter alpha
               alpha += delta_alpha; if (alpha > 1.) alpha = 1.; // pb d'erreur d'arrondi
               // on initialise x en tenant compte du delta_x de l'incrément précédent (dons on extrapole)
               // cela peut-être une bonne idée comme une mauvaise, car cela peut donner
               // des valeurs abérantes, mais normalement on ne fera qu'une passe car on aura une valeur de test
               // de retour des résidus et dérivées en 0 donc avec une erreur fatale
               if (delta_x_prec.Max_val_abs() >= ConstMath::trespetit)
                   x += delta_x_prec * (delta_alpha / delta_alpha_prec);
              }   
          };     
        // -- dans le cas où phase_convergence==0 c-a-d non convergence normal, on refait une boucle
        //    ce qui a pour conséquence d'augmenter iter, et de préciser la racine, c'est tout       
        } while (phase_convergence < 2) ;
   
     racine = x;
     // mise à jour de l'indicateur du nombre d'incréments
     nb_incr_total = incre;
     nb_iter_total += iter; // idem pour les itérations
     // retour de la validitée de la recherche du zéro
     if (phase_convergence == 2)
        return true;
     else
        return false;   
     }
 
    // idem précédemment, mais pour une fonction à valeur vectorielle et ..
    //  la méthode externe calcul la valeur de la fonction et de la dérivée en même temps
    // les matrices sont telles que: der_at_racine(i,j) = df(i)/d(x(j), soit: df = der_at_racine * dx
    // recherche du zéro d'une fonction en utilisant la méthode de Newton-Raphson incrémentale
    // ici il s'agit d'une fonction à n variables
    //  *Ptfonc : le pointeur de la fonction dont il faut chercher le zéro
    //  *Ptder_fonc : pointeur de la dérivée de la fonction, calcul également en même temps le vecteur résidu 
    //                c'est à dire la fonction que l'on veut annuler
    //               pour ces deux fonctions, la variable alpha est un facteur de charge qui varie de 
    //               0 à 1. Lorsque alpha=0, la racine vaut val_initiale, et ce que l'on 
    //               cherche c'est la racine pour alpha = 1. Lorsque alpha varie progressivement
    //               de 0 à 1, la racine est sensée varier progressivement de val_initiale à résidu
    //               l'utilisation d'alpha permet de faire du Newton incrémentale. 
    //   pour ces deux fonctions, l'argument test ramène
    //     . 1 si le calcul a été ok, -1 s'il y a eu un pb, mais on peut continuer, 0 s'il y a eu un pb
    //       fatal, qui invalide le calcul de la fonction et/ou de la dérivée
    //  val_initiale : une valeur initiale de x pour la recherche de zéro
    //  max_delta_x : si > 0 donne le norme maxi permise sur delta_x au cours d'une itération
    //               si ||delta_x|| est plus grand on fait delta_x = max_delta_x * delta_x / || delta_x||
    //  ramène en sortie:
    //     un booléen qui indique si la résolution est correcte ou non
    //     racine : la racine trouvée
    //     der_at_racine : contient en retour la valeur de la dérivée pour la racine trouvée
    //     nb_incr_total : le nombre total d'incrément qui a été nécessaire
    //     nb_iter_total : le nombre total d'itération, qui cumule les iter de tous les incréments
   template <class T>
   bool  Algo_zero::Newton_raphson
               (T& instance, Vecteur& (T::*Ptfonc) (const double & alpha,const Vecteur & x,int& test)
                ,Mat_abstraite& (T::*Ptder_fonc) (const double & alpha,const Vecteur & x,Vecteur& res,int& test)
                ,const Vecteur& val_initiale,Vecteur & racine,Mat_abstraite & der_at_racine
                ,int& nb_incr_total, int& nb_iter_total,double max_delta_x) const
    { int dime = val_initiale.Taille();    
      // initialisation des différentes variables intermédiaires
      nb_iter_total = nb_incr_total = 0;
      int iter = 0; // nombre d'itérations
      int incre = 1;  // nombre d'incrément pour converger, au débu = 1
      Vecteur x = val_initiale; // initialisation de la valeur courante de x
      int phase_convergence = 0; // a priori on est en phase de convergence
          //  = -3 delta_x trop grand ;   = -2 delta_x trop petit ; = -1 trop d'itération
          //  = 0  non convergence normal; = 1 convergence sur un incrément
          //  = 2  convergence finale
      double alpha = 1; // initialisation du facteur de charge
      double delta_alpha = 1; // incrément d'alpha, intit à 1
      double alpha_prec = 0; 
      Vecteur x_prec = x; // valeur précédente en début d'incrément
      double delta_alpha_prec = 1; // le delta alpha de l'incrément précédent
      Vecteur delta_x_prec(dime); // le delta x de l'incrément précédent 
      // grandeurs de calcul
      Vecteur residu(dime),delta_x(dime),xmoinsdeltax(dime);
      int test =1; // init par défaut
      int glob_err_res = 0;// pour la gestion de test == -1
      int glob_err_der = 0;// pour la gestion de test == -1
      double residu_Max_val_abs=ConstMath::tresgrand; // init

      
      do // boucle globale
       {      
        // 1  ----- calcul de la fonction et test de convergence
        residu = (instance.*Ptfonc) (alpha,x,test);
        if (  (test == 1) // c-a-d que le calcul est correct
            ||(test == -1) // petit pb mais on continue
           )
         {
          residu_Max_val_abs = residu.Max_val_abs();
          #ifdef MISE_AU_POINT
           if ((!isfinite(residu_Max_val_abs)) || (isnan(residu_Max_val_abs)))
            if (((permet_affichage==0) && (ParaGlob::NiveauImpression() >1))|| (permet_affichage > 1))
              { cout << "\n Algo_zero::Newton_raphson: iteration: " << iter  ;
                cout << "\n residu: "; residu.Affiche();
                cout << "\n *** pb resolution residu infini " ;
                cout << "\n Algo_zero::Newton_raphson( " << flush;
                throw ErrNonConvergence_Newton();
              };
          #endif
          if (test == -1) {glob_err_res++;} else {glob_err_res=0;}; // incrémentation ou remise à 0
          // -- manip sur les précisions
          // si prec_res est négatif, il s'agit d'une précision relative
   //        double prec_res_reelle = prec_res;
          double prec_res_reelle = prec_res_abs + residu_Max_val_abs * prec_res_rel;
   //        if (prec_res < 0) // si l'indication est négative , cela veut dire que l'on veut une grandeur négative
   //          prec_res_reelle = MaX(-prec_res * residu_Max_val_abs,-prec_res);
          if (residu_Max_val_abs > prec_res_reelle)
             { // cas où l'on n'a pas encore convergé
               der_at_racine = (instance.*Ptder_fonc) (alpha,x,residu,test); // calcul de la dérivée de la fonction
                                                                             // + la fonction elle-même donc le nouveau résidu
               if (test == -1) {glob_err_der++;} else {glob_err_der=0;}; // incrémentation ou remise à 0
               #ifdef MISE_AU_POINT
                if (((permet_affichage==0) && (ParaGlob::NiveauImpression() >=8))|| (permet_affichage >= 8))
                  { cout << "\n Algo_zero::Newton_raphson: iteration: " << iter  ;
                    cout << "\n residu: "; residu.Affiche();
                    cout << "\n matrice tangente: " ; der_at_racine.Affiche();
                    cout << flush;
                  };
               #endif
       
              if ( (test == 1) // cas sans pb
                  || ((test == -1) && ((glob_err_res+glob_err_der) < maxi_test_moins1)) // test non fatal et on est < à maxi_test_moins1 séries de test mauvais
                 )                                   // en supposant que l'on a successivement le résidu et la dérivée en
                                                     // mauvais calcul, sinon ce sera l'un ou l'autre
                { // détermination du delta_x
                  delta_x = - der_at_racine.Resol_systID_2 (residu,delta_x);
                  if ((!isfinite(delta_x(1))) || (isnan(delta_x(1))))
                     { cout << "\n *** pb resolution delta_x est infini " << delta_x  ;
                       cout << "\n Algo_zero::Newton_raphson( " << flush;
      //                 //--- debug
      //                 cout << "\n debug Algo_zero::Newton_raphson( ";
      //                 residu = (instance.*Ptfonc) (alpha,x);
      //                 der_at_racine = (instance.*Ptder_fonc) (alpha,x,residu); // calcul de la dérivée de la fonction
      //                 cout << "\n iteration: " << iter  ;
      //                 cout << "\n residu: "; residu.Affiche();
      //                 cout << "\n matrice tangente: " ; der_at_racine.Affiche();
      //                 cout << flush;
      //                 // --- fin debug

                       throw ErrNonConvergence_Newton();
                     };
                  // on applique éventuellement une limitation
                  double val_abs_delta_x = delta_x.Norme();
                  if ((max_delta_x > ConstMath::trespetit)&&(val_abs_delta_x > max_delta_x))
                   {double limitation = max_delta_x / val_abs_delta_x;
                    delta_x *=  limitation;
                    #ifdef MISE_AU_POINT
                     if (((permet_affichage==0) && (ParaGlob::NiveauImpression() > 4)) || (permet_affichage > 2))
                       { cout << "\n Algo_zero::Newton_raphson: limitation de delta_x: " << limitation;
                       };
                    #endif
                   };

                  x = x + delta_x; // nouvelle valeur de x
                  #ifdef MISE_AU_POINT
                   if (((permet_affichage==0) && (ParaGlob::NiveauImpression() >=8))|| (permet_affichage >= 8))
                     { cout << "\n  delta_X :"; delta_x.Affiche();
                       cout << "\n x :"  ;x.Affiche();
                       cout << flush ;
                     };
                  #endif
                  iter++; // incrémentation du nombre d'itération
                  // on étudie de nouveau la convergence, car ici le résidu a été recalculé
                  // en même temps que le calcul de l'opérateur tangent
                  residu_Max_val_abs = residu.Max_val_abs();
                  prec_res_reelle = prec_res_abs + residu_Max_val_abs * prec_res_rel;
      //            if (prec_res < 0) // si l'indication est négative , cela veut dire que l'on veut une grandeur négative
      //               prec_res_reelle = MaX(-prec_res * residu_Max_val_abs,-prec_res);
                  if (residu_Max_val_abs < prec_res_reelle)
                    { phase_convergence = 1;
                      // on sauvegarde les variables
                      delta_x_prec = x - x_prec;
                      x_prec = x;
                      alpha_prec = alpha;
                      delta_alpha_prec = delta_alpha;
                    }
                  // == étude des divergences éventuelles
                  else
                   {xmoinsdeltax = x-delta_x;
                    if (iter > iter_max) {phase_convergence=-1;} // trop d'itération
                    else if (delta_x.Max_val_abs() < mini_delta_x)
                        { phase_convergence=-2; } // delta_x est rendu trop petit
                    else if ((delta_x.Max_val_abs()
                              < coef_mini_delta_x * (MaX( xmoinsdeltax.Max_val_abs(),x.Max_val_abs())) ) )
                        { phase_convergence=-4; } // delta_x est rendu trop petit
                    else if (delta_x.Max_val_abs() > maxi_delta_x)
                      { phase_convergence = -3;} // delta_x est rendu trop grand
                    else // sinon c'est une non convergence normale, et on continue à itérer
                     { phase_convergence = 0;};
                   };
                }
              else if ( (!test) || ((glob_err_res+glob_err_der+1) > maxi_test_moins1)) // si test == 0 -> !test = true
                {// que ce soit une erreur fatale ou une suite de mauvais calcul: on signale une divergence
                 phase_convergence = -4;
                }
              else // c'est une erreur non répertoriée
               {cout << "\n %% Algo_zero::Newton_raphson: cas non prevu de retour d'erreur sur la derivee =  "
                     << test << " %% arret algo Newton %% "
                     << " glob_err_res+glob_err_der+1= "<<glob_err_res+glob_err_der+1
                     << " maxi_test_moins1= "<<maxi_test_moins1
                     << flush;
                Sortie(1);
               };
             }
           else // cas ou on a convergé
             { phase_convergence = 1;
               // on sauvegarde les variables
               delta_x_prec = x - x_prec;
               x_prec = x;
               alpha_prec = alpha;
               delta_alpha_prec = delta_alpha;
               #ifdef MISE_AU_POINT
                if (((permet_affichage==0) && (ParaGlob::NiveauImpression() > 5))|| (permet_affichage > 5))
                  { cout << "\n Algo_zero::Newton_raphson: convergence pour " << iter
                         << " iteration(s) " << incre << " increment(s) "
                         << " , alpha=  " << alpha << " et delta_alpha= " << delta_alpha;
                    cout << "\n  residu: "; residu.Affiche();cout << flush;
                  };
               #endif
             };
          }
        else if (!test)
          {// on signale que l'on est en divergence due à une erreur fatale sur le résidu
           phase_convergence = -4;
           if (x == val_initiale) // cela veut dire que dès le début il y a une erreur fatale
            {if ((permet_affichage > 1) || (ParaGlob::NiveauImpression() > 3))
               cout << "\n %% Algo_zero::Newton_raphson: des le debut, erreur fatale sur le calcul du residu "
                    << " %% arret algo Newton %% " << endl;
             break; // on ne peut pas continuer
            };
          }
        else // c'est une erreur non répertoriée
         {cout << "\n %% Algo_zero::Newton_raphson: cas non prevu de retour d'erreur sur le residu =  "
               << test << " %% arret algo Newton %% " << flush;
          Sortie(1);
         };

        // 2 ----- action en fonction de la convergence ou non
        // -- cas d'une non convergence particulière c-a-d en dehors du cas  phase_convergence == 0  
        // dans le cas où on est en divergence on cherche à augmenter le nombre d'incréments
        if (phase_convergence < 0)
          { if (incre <= nbMaxiIncre)
              {
                #ifdef MISE_AU_POINT
                 if (((permet_affichage==0) && (ParaGlob::NiveauImpression() > 5))|| (permet_affichage > 5))
                   { cout << "\n Algo_zero::Newton_raphson: non convergence pour alpha=  " << alpha
                          << " et delta_alpha= " << delta_alpha
                          << " ( incre= "<<incre << " iter= "<<iter
                          << " ||residu||= "<<residu_Max_val_abs
                          << ") " ;
                     if (((permet_affichage==0) && (ParaGlob::NiveauImpression() > 6))|| (permet_affichage > 6))
                      cout << "\n phase_convergence= "<< phase_convergence
                           << " delta_x.Max_val_abs()= " << delta_x.Max_val_abs()
                           << "(mini_delta_x= "<<mini_delta_x<<" maxi_delta_x= "<<maxi_delta_x
                           << " iter_max= "<<iter_max << ") "
                           << " prec_res_reelle=" << (prec_res_abs + residu_Max_val_abs * prec_res_rel);
                     cout << flush;
                   };
                #endif
                // on diminue l'incrément de charge par deux
                delta_alpha *= 0.5;
                nb_iter_total += iter; // stockage du nombre d'iter
                iter = 0; // on reinitialise le nombre d'itération
                // on repart avec une valeur d'alpha plus petite
                alpha = alpha_prec + delta_alpha;
                x = x_prec;
                // on incrémente le compteur d'incréments
                incre++; 
              }
            else // on a dépassé le nombre d'incrément permis, on arrête
              {
                #ifdef MISE_AU_POINT
                if (((permet_affichage==0) && (ParaGlob::NiveauImpression() > 6))|| (permet_affichage > 6))
                   { cout << "\n Algo_zero::Newton_raphson: %% on atteind le nombre d'increment maxi  "
                          << nbMaxiIncre << " %% arret algo Newton %% " << flush;
                   };
                #endif
                break;
              };
          }
        else if (phase_convergence == 1)
         { // on vient de converger, on regarde si l'on a terminé les incréments de charge
           if (alpha >= 1.) 
             { phase_convergence = 2;} // on a terminé
           else
             { // sinon il faut incrémenter alpha
               alpha += delta_alpha; if (alpha > 1.) alpha = 1.; // pb d'erreur d'arrondi
               // on initialise x en tenant compte du delta_x de l'incrément précédent (donc on extrapole)
               // cela peut-être une bonne idée comme une mauvaise, car cela peut donner
               // des valeurs abérantes, mais normalement on ne fera qu'une passe car on aura une valeur de test
               // de retour des résidus et dérivées en 0 donc avec une erreur fatale
               x += delta_x_prec * (delta_alpha / delta_alpha_prec);
             };
         };
        // -- dans le cas où phase_convergence==0 c-a-d non convergence normal, on refait une boucle
        //    ce qui a pour conséquence d'augmenter iter, et de préciser la racine, c'est tout       
       } while (phase_convergence <= 1) ;
   
     racine = x;
     // mise à jour de l'indicateur du nombre d'incréments
     nb_incr_total = incre;
     nb_iter_total += iter; // idem pour les itérations
     // retour de la validitée de la recherche du zéro
     if (phase_convergence == 2)
        return true;
     else
        return false;   
     }
     
 #endif