/*
 * Copyright (C) 2007-2008 Dynare Team
 *
 * This file is part of Dynare.
 *
 * Dynare 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.
 *
 * Dynare 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 Dynare.  If not, see <http://www.gnu.org/licenses/>.
 */

#include <cstdlib>
#include <cstdio>
#include <iostream>
#include <string>
#include <ctime>
#include <stack>
#include <cmath>
#include <fstream>

#include "ModelTree.hh"
//#include "pctimer_h.hh"
#include "Model_Graph.hh"
#include "SymbolGaussElim.hh"
//#define SIMPLIFY

using namespace std;
int max_nb_table_y, max_nb_in_degree_edges=0, max_nb_out_degree_edges=0;

SymbolicGaussElimination::SymbolicGaussElimination()
{
#ifdef SAVE
  file_open = false;
#endif
}

#ifdef SIMPLIFY
void
SymbolicGaussElimination::list_table_u(int pos)
{
  int i;
  t_table_u *table_u;
  table_u = First_table_u->pNext;
  i = 0;
  while(table_u)
    {
      if((table_u->type > 7) || (table_u->type < 1))
        {
          cout << "Error : table_u->type=" << int(table_u->type) << " i=" << i << " pos=" << pos << "\n";
          system("pause");
          exit( -1);
        }
      i++;
      table_u = table_u->pNext;
    }
}

//==================================================================================
void
SymbolicGaussElimination::print_free_u_list()
{
  int i;
  cout << "nb_free_u_list : " << nb_free_u_list << " \n";
  cout << "-----------------------------------------------\n";
  for(i = 0;i < nb_free_u_list;i++)
    cout << i << " " << free_u_list[i] << "\n";
}

//==================================================================================
void
SymbolicGaussElimination::set_free_u_list(int index)
{
#ifndef UNSORTED
  int i = 0, j;
  if(nb_free_u_list > 0)
    {
      while((free_u_list[i] > index) && (i < nb_free_u_list))
        i++;
      for(j = nb_free_u_list;j > i;j--)
        {
          free_u_list[j] = free_u_list[j - 1];
          if(j >= MAX_FREE_U_LIST)
            {
              cout << "Error to much j=" << j << " Maximum allowed=" << MAX_FREE_U_LIST << "\n";
              system("pause");
              exit( -1);
            }
        }
    }
  if(i >= MAX_FREE_U_LIST)
    {
      cout << "Error to much i=" << i << " Maximum allowed=" << MAX_FREE_U_LIST << "\n";
      system("pause");
      exit( -1);
    }
  free_u_list[i] = index;
  nb_free_u_list++;
  if(nb_free_u_list >= MAX_FREE_U_LIST)
    {
      cout << "Error to much free_u_list=" << nb_free_u_list << " Maximum allowed=" << MAX_FREE_U_LIST << "\n";
      system("pause");
      exit( -1);
    }
  if(nb_free_u_list > max_nb_free_u_list)
    max_nb_free_u_list = nb_free_u_list;
#else
  free_u_list[nb_free_u_list++] = index;
#endif
}

int
SymbolicGaussElimination::get_free_u_list(bool dynamic)
{
  int i;
  if((nb_free_u_list > 0) && (simplification_allowed))
    {
      i = free_u_list[--nb_free_u_list];
      return (i);
    }
  else
    if(dynamic)
      return (u_count++);
    else
      return (u_count++);
}
#endif /**SIMPLIFY**/

//==================================================================================
void
SymbolicGaussElimination::init_glb()
{
  tol = TOL;
  vertex_count = 0;
  nb_table_u = 0;
  nb_prologue_save_table_y = 0;
  nb_first_save_table_y = 0;
  nb_middle_save_table_y = 0;
  nb_last_save_table_y = 0;
  nb_prologue_save_table_u = 0;
  nb_first_save_table_u = 0;
  nb_middle_save_table_u = 0;
  nb_last_save_table_u = 0;
  middle_count_loop = 0;
  u_count = 0;
  y_kmin = y_kmax = 0;
  prologue_save_table_u = NULL;
  first_save_table_u = NULL;
  first_save_i_table_u = NULL;
  middle_save_table_u = NULL;
  middle_save_i_table_u = NULL;
  last_save_table_u = NULL;
  prologue_save_table_y = NULL;
  first_save_table_y = NULL;
  first_save_i_table_y = NULL;
  middle_save_table_y = NULL;
  middle_save_i_table_y = NULL;
  last_save_table_y = NULL;
  pos_nb_first_save_table_u = 0;
#ifdef SIMPLIFY
  nb_free_u_list = 0;
  max_nb_free_u_list1 = 0;
  max_nb_free_u_list = 0;
#endif
  save_direct=true;
}


void
SymbolicGaussElimination::write_to_file_table_y( t_table_y *save_table_y, t_table_y *save_i_table_y, int nb_save_table_y, int nb_save_i_table_y)
{
  int i, k;
#ifdef PRINT_OUT
  cout << "nb_save_table_y=" << nb_save_table_y << "\n";
#endif
  SaveCode.write(reinterpret_cast<char *> (&nb_save_table_y), sizeof(nb_save_table_y));
  for(i = 0;i < nb_save_table_y;i++)
    {
      SaveCode.write(reinterpret_cast<char *>(&save_table_y[i].index), sizeof(save_table_y[i].index));
#ifdef PRINT_OUT
      cout << "+y[" << save_table_y[i].index << "]=";
#endif
      SaveCode.write(reinterpret_cast<char *>(&save_table_y[i].nb), sizeof(save_table_y[i].nb));
      for(k = 0;k < save_table_y[i].nb;k++)
        {
          SaveCode.write(reinterpret_cast<char *>(&save_table_y[i].u_index[k]), sizeof(save_table_y[i].u_index[k]));
          SaveCode.write(reinterpret_cast<char *>(&save_table_y[i].y_index[k]), sizeof(save_table_y[i].y_index[k]));
#ifdef PRINT_OUT
          cout << "u[" << save_table_y[i].u_index[k] << "]*y[" << save_table_y[i].y_index[k] << "]";
          if(k + 1 < save_table_y[i].nb)
            cout << "+";
          else
            cout << "\n";
#endif
        }
    }
  for(i = 0;i < nb_save_i_table_y;i++)
    {
      SaveCode.write(reinterpret_cast<char *>(&save_i_table_y[i].index), sizeof(save_i_table_y[i].index));
#ifdef PRINT_OUT
      cout << "+i y[" << save_i_table_y[i].index << "]=";
#endif
      SaveCode.write(reinterpret_cast<char *>(&save_i_table_y[i].nb), sizeof(save_i_table_y[i].nb));
      for(k = 0;k < save_i_table_y[i].nb;k++)
        {
          SaveCode.write(reinterpret_cast<char *>(&save_i_table_y[i].u_index[k]), sizeof(save_i_table_y[i].u_index[k]));
          SaveCode.write(reinterpret_cast<char *>(&save_i_table_y[i].y_index[k]), sizeof(save_i_table_y[i].y_index[k]));
#ifdef PRINT_OUT
          cout << "u[" << save_i_table_y[i].u_index[k] << "]*y[" << save_i_table_y[i].y_index[k] << "]";
          if(k + 1 < save_i_table_y[i].nb)
            cout << "+";
          else
            cout << "\n";
#endif
        }
    }
}

void
SymbolicGaussElimination::write_to_file_table_u(t_table_u *save_table_u, t_table_u *save_i_table_u, int nb_save_table_u)
{
  t_table_u *table_u;
  SaveCode.write(reinterpret_cast<char *>(&nb_save_table_u), sizeof(nb_save_table_u));
#ifdef PRINT_OUT
  cout << "**nb_save_table_u=" << nb_save_table_u << "\n";
#endif
#ifdef PRINT_OUT
  cout << "**save_table_u=" << save_table_u << "\n";
#endif
  while(save_table_u)
    {
#ifdef PRINT_OUT
      cout << "**save_table_u->type=" << int(save_table_u->type) << "\n";
#endif
      switch (save_table_u->type)
        {
        case 3:
        case 7:
          SaveCode.write(reinterpret_cast<char *>(&save_table_u->type), sizeof(save_table_u->type));
          SaveCode.write(reinterpret_cast<char *>(&save_table_u->index), sizeof(save_table_u->index));
          SaveCode.write(reinterpret_cast<char *>(&save_table_u->op1), sizeof(save_table_u->op1));
#ifdef PRINT_OUT
          if(save_table_u->type == 3)
            cout << "+u[" << save_table_u->index << "]=1/(1-u[" << save_table_u->op1 << "])\n";
          else
            cout << "+u[" << save_table_u->index << "]*=u[" << save_table_u->op1 << "]\n";
#endif /**PRINT_OUT**/
          break;
        case 1:
        case 2:
        case 6:
          SaveCode.write(reinterpret_cast<char *>(&save_table_u->type), sizeof(save_table_u->type));
          SaveCode.write(reinterpret_cast<char *>(&save_table_u->index), sizeof(save_table_u->index));
          SaveCode.write(reinterpret_cast<char *>(&save_table_u->op1), sizeof(save_table_u->op1));
          SaveCode.write(reinterpret_cast<char *>(&save_table_u->op2), sizeof(save_table_u->op2));
#ifdef PRINT_OUT
          if(save_table_u->type == 1)
            cout << "+u[" << save_table_u->index << "]=" << "u[" << save_table_u->op1 << "]*u[" << save_table_u->op2 << "]\n";
          else if(save_table_u->type == 2)
            cout << "+u[" << save_table_u->index << "]+=u[" << save_table_u->op1 << "]*u[" << save_table_u->op2 << "]\n";
          else
            cout << "+u[" << save_table_u->index << "]=1/(1-u[" << save_table_u->op1 << "]*u[" << save_table_u->op2 << "])\n";
#endif /**PRINT_OUT**/
          break;
        case 5:
          SaveCode.write(reinterpret_cast<char *>(&save_table_u->type), sizeof(save_table_u->type));
          SaveCode.write(reinterpret_cast<char *>(&save_table_u->index), sizeof(save_table_u->index));
#ifdef PRINT_OUT
          cout << "+push(u[" << save_table_u->index << "])\n";
#endif /**PRINT_OUT**/
          break;
        }
      table_u = save_table_u->pNext;
      free(save_table_u);
      save_table_u = table_u;
    }
#ifdef PRINT_OUT
  cout << "save_i_table_u=" << save_i_table_u << "\n";
#endif
  while(save_i_table_u)
    {
      switch (save_i_table_u->type)
        {
        case 3:
        case 7:
          SaveCode.write(reinterpret_cast<char *>(&save_i_table_u->type), sizeof(save_i_table_u->type));
          SaveCode.write(reinterpret_cast<char *>(&save_i_table_u->index), sizeof(save_i_table_u->index));
          SaveCode.write(reinterpret_cast<char *>(&save_i_table_u->op1), sizeof(save_i_table_u->op1));
          break;
        case 1:
        case 2:
        case 6:
          SaveCode.write(reinterpret_cast<char *>(&save_i_table_u->type), sizeof(save_i_table_u->type));
          SaveCode.write(reinterpret_cast<char *>(&save_i_table_u->index), sizeof(save_i_table_u->index));
          SaveCode.write(reinterpret_cast<char *>(&save_i_table_u->op1), sizeof(save_i_table_u->op1));
          SaveCode.write(reinterpret_cast<char *>(&save_i_table_u->op2), sizeof(save_i_table_u->op2));
          break;
        case 5:
          SaveCode.write(reinterpret_cast<char *>(&save_i_table_u->type), sizeof(save_i_table_u->type));
          SaveCode.write(reinterpret_cast<char *>(&save_i_table_u->index), sizeof(save_i_table_u->index));
          break;
        }
      table_u = save_i_table_u->pNext;
      free(save_i_table_u);
      save_i_table_u = table_u;
    }
#ifdef PRINT_OUT
  cout << "nb_save_table_u=" << nb_save_table_u << "\n";
#endif
}

void
SymbolicGaussElimination::write_to_file_table_u_b(t_table_u *save_table_u, t_table_u *last_table_u, int *nb_save_table_u, bool chk)
{
  t_table_u *table_u;
  bool OK=true;
  while(/*save_table_u!=last_table_u*/OK)
    {
#ifdef PRINT_OUT
      cout << "**save_table_u->type=" << int(save_table_u->type) << "\n";
#endif
      switch (save_table_u->type)
        {
        case 3:
        case 7:
          (*nb_save_table_u)++;
          //SaveCode << save_table_u->type << save_table_u->index << save_table_u->op1;
          SaveCode.write(reinterpret_cast<char *>(&save_table_u->type), sizeof(save_table_u->type));
          SaveCode.write(reinterpret_cast<char *>(&save_table_u->index), sizeof(save_table_u->index));
          SaveCode.write(reinterpret_cast<char *>(&save_table_u->op1), sizeof(save_table_u->op1));
#ifdef PRINT_OUT
          if(save_table_u->type == 3)
            cout << "+u[" << save_table_u->index << "]=1/(1-u[" << save_table_u->op1 << "])\n";
          else
            cout << "+u[" << save_table_u->index << "]*=u[" << save_table_u->op1 << "]\n";
#endif /**PRINT_OUT**/
          break;
        case 1:
        case 2:
        case 6:
          (*nb_save_table_u)++;
          //SaveCode << save_table_u->type << save_table_u->index << save_table_u->op1 << save_table_u->op2;
          SaveCode.write(reinterpret_cast<char *>(&save_table_u->type), sizeof(save_table_u->type));
          SaveCode.write(reinterpret_cast<char *>(&save_table_u->index), sizeof(save_table_u->index));
          SaveCode.write(reinterpret_cast<char *>(&save_table_u->op1), sizeof(save_table_u->op1));
          SaveCode.write(reinterpret_cast<char *>(&save_table_u->op2), sizeof(save_table_u->op2));
#ifdef PRINT_OUT
          if(save_table_u->type == 1)
            cout << "+u[" << save_table_u->index << "]=" << "u[" << save_table_u->op1 << "]*u[" << save_table_u->op2 << "]\n";
          else if(save_table_u->type == 2)
            cout << "+u[" << save_table_u->index << "]+=u[" << save_table_u->op1 << "]*u[" << save_table_u->op2 << "]\n";
          else
            cout << "+u[" << save_table_u->index << "]=1/(1-u[" << save_table_u->op1 << "]*u[" << save_table_u->op2 << "])\n";
#endif /**PRINT_OUT**/
          break;
        case 5:
          (*nb_save_table_u)++;
          //SaveCode << save_table_u->type << save_table_u->index;
          SaveCode.write(reinterpret_cast<char *>(&save_table_u->type), sizeof(save_table_u->type));
          SaveCode.write(reinterpret_cast<char *>(&save_table_u->index), sizeof(save_table_u->index));
#ifdef PRINT_OUT
          cout << "+push(u[" << save_table_u->index << "])\n";
#endif /**PRINT_OUT**/
          break;
        }
      if(chk)
        {
          OK=(save_table_u!=last_table_u);
          if(OK)
            {
              table_u = save_table_u->pNext;
              free(save_table_u);
              save_table_u = table_u;
            }
        }
      else
        {
          table_u = save_table_u->pNext;
          free(save_table_u);
          save_table_u = table_u;
          OK=(save_table_u!=last_table_u);
        }
    }
}


void
store_code(t_table_u **save_table_u, t_table_u *First_table_u, t_table_u *stop_table_u, t_table_y *table_y, t_table_y **s_table_y, int *nb_save_table_u, int *nb_save_table_y, int vertex_count)
{
  int i, k;
  t_table_u *s_table_u, *table_u;
  table_u = First_table_u->pNext;
  s_table_u = (t_table_u*)malloc(sizeof(t_table_u) - 3 * sizeof(int));
  (*save_table_u) = s_table_u;
  (*nb_save_table_u) = 0;
  while(table_u != stop_table_u)
    {
      if((table_u->type > 7) || (table_u->type < 1))
        {
          cout << "Error : table_u->type=" << int(table_u->type) << " *nb_save_table_u=" << *nb_save_table_u << "\n";
          exit( -1);
        }
      else
        {
          (*nb_save_table_u)++;
#ifdef PRINT_OUT
          cout << "--table_u->type=" << int(table_u->type) << "\n";
#endif
          switch (table_u->type)
            {
            case 3:
            case 7:
              s_table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u) - sizeof(int));
              s_table_u = s_table_u->pNext;
              s_table_u->type = table_u->type;
              s_table_u->index = table_u->index;
              s_table_u->op1 = table_u->op1;
#ifdef PRINT_OUT
              if(s_table_u->type == 3)
                cout << "st u[" << s_table_u->index << "]=-1/u[" << s_table_u->op1 << "]\n";
              else
                cout << "st u[" << s_table_u->index << "]*=u[" << s_table_u->op1 << "]\n";
#endif
              break;
            case 1:
            case 2:
            case 6:
              s_table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u) );
              s_table_u = s_table_u->pNext;
              s_table_u->type = table_u->type;
              s_table_u->index = table_u->index;
              s_table_u->op1 = table_u->op1;
              s_table_u->op2 = table_u->op2;
#ifdef PRINT_OUT
              if(s_table_u->type == 1)
                cout << "st u[" << s_table_u->index << "]=" << "u[" << s_table_u->op1 << "]*u[" << s_table_u->op2 << "]\n";
              else if(s_table_u->type == 2)
                cout << "st u[" << s_table_u->index << "]+=u[" << s_table_u->op1 << "]*u[" << s_table_u->op2 << "]\n";
              else
                cout << "st u[" << s_table_u->index << "]=1/(1-u[" << s_table_u->op1 << "]*u[" << s_table_u->op2 << "])\n";
#endif /**PRINT_OUT**/
              break;
            case 5:
              s_table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u) - 2 * sizeof(int));
              s_table_u = s_table_u->pNext;
              s_table_u->type = table_u->type;
              s_table_u->index = table_u->index;
#ifdef PRINT_OUT
              cout << "st push(u[" << s_table_u->index << "])\n";
#endif /**PRINT_OUT**/
              break;
            }
        }
      table_u = table_u->pNext;
    }
  s_table_u->pNext = NULL;
#ifdef PRINT_OUT
  cout << "*nb_save_table_u=" << *nb_save_table_u << "\n";
  cout << "vertex_count=" << vertex_count << " nb_save_table_y=" << nb_save_table_y << "\n";
#endif
  (*s_table_y) = (t_table_y*)malloc((vertex_count - *nb_save_table_y) * sizeof(t_table_y));
  for(i = 0;i < vertex_count - *nb_save_table_y;i++)
    {
      (*s_table_y)[i].u_index = (int*)malloc(table_y[i + *nb_save_table_y].nb * sizeof(int));
      (*s_table_y)[i].y_index = (int*)malloc(table_y[i + *nb_save_table_y].nb * sizeof(int));
      /*(*s_table_y)[i].y_lead_lag = (int*)malloc(table_y[i + *nb_save_table_y].nb * sizeof(int));*/
      (*s_table_y)[i].index = table_y[i + *nb_save_table_y].index;
      (*s_table_y)[i].nb = table_y[i + *nb_save_table_y].nb;
      for(k = 0;k < table_y[i + *nb_save_table_y].nb;k++)
        {
          (*s_table_y)[i].u_index[k] = table_y[i + *nb_save_table_y].u_index[k];
          (*s_table_y)[i].y_index[k] = table_y[i + *nb_save_table_y].y_index[k];
        }
    }
  *nb_save_table_y = vertex_count - *nb_save_table_y;
}



bool
SymbolicGaussElimination::Check_Regularity(t_table_u *first_u_blck, t_table_u *second_u_blck, t_table_u *third_u_blck)
{
  t_table_u *c_first_table_u, *c_second_table_u, *c_third_table_u;
  bool OK, Regular;
  int i=0;
  c_first_table_u = first_u_blck->pNext;
  c_second_table_u = second_u_blck->pNext;
  c_third_table_u = third_u_blck->pNext;
  OK = true;
  Regular = true;
  while(OK)
    {
      i++;
#ifdef PRINT_OUT
      cout << "c_first_table_u=" << c_first_table_u << " c_second_table_u=" << c_second_table_u << " c_third_table_u=" << c_third_table_u << "\n";
#endif /**PRINT_OUT**/
      if((c_first_table_u->type != c_second_table_u->type) ||
         (c_second_table_u->index + (c_second_table_u->index - c_first_table_u->index) != c_third_table_u->index))
        {
          Regular = false;
#ifdef PRINT_OUT
          cout << "c_first_table_u->type=" << (int)c_first_table_u->type << "?=c_second_table_u->type=" << (int)c_second_table_u->type << "\n";
          cout << "c_second_table_u->index+(c_second_table_u->index-c_first_table_u->index)=" << c_second_table_u->index + (c_second_table_u->index - c_first_table_u->index)
               << "?=c_third_table_u->index=" << c_third_table_u->index << "\n";
#endif
          break;
        }
      if(c_first_table_u->type != 5)
        {
          if(c_second_table_u->op1 + (c_second_table_u->op1 - c_first_table_u->op1) != c_third_table_u->op1)
            {
#ifdef PRINT_OUT
              cout << "c_second_table_u->op1+(c_second_table_u->op1-c_first_table_u->op1)=" << c_second_table_u->op1 + (c_second_table_u->op1 - c_first_table_u->op1)
                   << "?=c_third_table_u->op1=" << c_third_table_u->op1 << "\n";
#endif
              Regular = false;
              break;
            }
          if((c_first_table_u->type != 3) && (c_first_table_u->type != 7))
            {
              if(c_second_table_u->op2 + (c_second_table_u->op2 - c_first_table_u->op2) != c_third_table_u->op2)
                {
#ifdef PRINT_OUT
                  cout << "c_second_table_u->op2+(c_second_table_u->op2-c_first_table_u->op2)=" << c_second_table_u->op2 + (c_second_table_u->op2 - c_first_table_u->op2)
                       << "?=c_third_table_u->op2=" << c_third_table_u->op2 << "\n";
#endif
                  Regular = false;
                  break;
                }
            }
        }
      if(c_first_table_u->pNext != second_u_blck->pNext)
        {
          c_second_table_u = c_second_table_u->pNext;
          c_first_table_u = c_first_table_u->pNext;
          c_third_table_u = c_third_table_u->pNext;
        }
      else
        OK = 0;
    }
#ifdef PRINT_OUT
  cout << "Regular=" << Regular << " OK=" << OK << " i=" << i << "\n";
  cout << "Check_Regularity Regular=" << Regular << "\n";
#endif
  return (Regular);
}


//==================================================================================
t_table_u*
SymbolicGaussElimination::interpolation(t_model_graph* model_graph, t_table_y* table_y, int to_add, bool middle, int per)
{
  t_table_u *tmp_table_u=NULL, *old_table_u=NULL;
  t_table_u *c_first_table_u, *c_second_table_u, *c_third_table_u;
  int c_first_y_blck, c_second_y_blck;
  int i, j, k, up_to=0, op_count, k1, k2;
  bool OK, modify_u_count;
  int cur_pos, nb_table_u=0;
#ifdef PRINT_OUT
  cout << "in interpolation\n";
  cout << "--- u_count=" << u_count << "\n";
  print_Graph(model_graph);
#endif /**PRINT_OUT**/
  cur_pos=SaveCode.tellp();
  SaveCode.write(reinterpret_cast<char *>(&up_to), sizeof(up_to));
  if(middle)
    {
      up_to = 2;
      middle_count_loop = 1 + per;
      //middle_count_loop = 4;
    }
  else
    up_to = 2;
#ifdef SAVE
  if(middle)
    {
      middle_save_table_u = NULL;
#ifdef PRINT_OUT
      cout << "up_to=" << up_to << " middle_count_loop=" << middle_count_loop << "\n";
#endif
    }
#endif /**SAVE**/
  /*Adding the Gaussian Elimination coefficients for the uncomputed periods*/
  /*table_u = tmp_table_u = (t_table_u*)malloc(sizeof(t_table_u) - 3 * sizeof(int));*/
  c_third_table_u = third_u_blck;
  c_first_y_blck = first_y_blck;
  c_second_y_blck = second_y_blck;
  modify_u_count = true;
  for(i = 1;i < up_to;i++)
    {
      c_first_table_u = first_u_blck->pNext;
      c_second_table_u = second_u_blck->pNext;
#ifdef PRINT_OUT
      cout << "c_first_table_u=" << c_first_table_u << " c_second_table_u=" << c_second_table_u << "\n";
#endif
      old_table_u = tmp_table_u;
      op_count = 0;
      OK = true;
      while(OK)
        {
          if(c_first_table_u->type != c_second_table_u->type)
            {
              cout << "Error : Interpolation the two lists are not synchronized (middle=" << middle << ", op_count=" << op_count << ")\n";
              cout << "c_first_table_u=" << c_first_table_u << " c_second_table_u=" << c_second_table_u << "\n";
              cout << "y_kmin=" << y_kmin << " y_kmax=" << y_kmax << "\n";
              cout << "first_u_blck=" << first_u_blck << " second_u_blck=" << second_u_blck << "\n";
              cout << "c_first_table_u->type=" << int(c_first_table_u->type) << " c_second_table_u->type=" << int(c_second_table_u->type) << "\n";
              system("pause");
              exit( -1);
            }
          switch (c_first_table_u->type)
            {
            case 3:
            case 7:
              nb_table_u++;
#ifdef SAVE
              if(i == 1)
                {
                  if(middle)
                    {
                      nb_middle_save_table_u++;
                    }
                  else
                    {
#ifdef PRINT_OUT
                      cout << "!!!!!nb_first_save_table_u=" << nb_first_save_table_u << "\n";
#endif
                      nb_first_save_table_u++;
                    }
                  SaveCode.write(reinterpret_cast<char *>(&c_first_table_u->type), sizeof(c_first_table_u->type));
                  SaveCode.write(reinterpret_cast<char *>(&c_first_table_u->index), sizeof(c_first_table_u->index));
                  SaveCode.write(reinterpret_cast<char *>(&c_first_table_u->op1), sizeof(c_first_table_u->op1));
#ifdef PRINT_OUT
                  if(c_first_table_u->type == 3)
                    cout << ">u[" << c_first_table_u->index << "]=1/(1-u[" << c_first_table_u->op1 << "])\n"; //"]) (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n";
                  else
                    cout << ">u[" << c_first_table_u->index << "]*=u[" << c_first_table_u->op1 << "]\n"; //"] (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n";
#endif
                  k=c_second_table_u->index - c_first_table_u->index;
                  SaveCode.write(reinterpret_cast<char *>(&k), sizeof(k));
                  k1=c_second_table_u->op1- c_first_table_u->op1;
                  SaveCode.write(reinterpret_cast<char *>(&k1), sizeof(k1));
#ifdef PRINT_OUT
                  if(c_first_table_u->type == 3)
                    cout << ">i u[" << k << "]=1/(1-u[" << k1 << "])\n"; //"]) (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n";
                  else
                    cout << ">i u[" << k << "]*=u[" << k1 << "]\n"; //"] (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n";
#endif
                }
#endif /**SAVE**/
              break;
            case 1:
            case 2:
            case 6:
              nb_table_u++;
#ifdef SAVE
              if((i == 1) && (middle))
                {
                  nb_middle_save_table_u++;
                  SaveCode.write(reinterpret_cast<char *>(&c_first_table_u->type), sizeof(c_first_table_u->type));
                  SaveCode.write(reinterpret_cast<char *>(&c_first_table_u->index), sizeof(c_first_table_u->index));
                  SaveCode.write(reinterpret_cast<char *>(&c_first_table_u->op1), sizeof(c_first_table_u->op1));
                  SaveCode.write(reinterpret_cast<char *>(&c_first_table_u->op2), sizeof(c_first_table_u->op2));
#ifdef PRINT_OUT
                  if((c_first_table_u->type == 1) && (middle))
                    cout << ">u[" << c_first_table_u->index << "]=u[" << c_first_table_u->op1 << "]*u[" << c_first_table_u->op2 << "]\n"; //"] (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n";
                  else if(c_first_table_u->type == 2)
                    cout << ">u[" << c_first_table_u->index << "]+=u[" << c_first_table_u->op1 << "]*u[" << c_first_table_u->op2 << "]\n"; //"] (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n";
                  else
                    cout << ">u[" << c_first_table_u->index << "]=1/(1-u[" << c_first_table_u->op1 << "]*u[" << c_first_table_u->op2 << "])\n"; //"]) (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n";
#endif //PRINT_OUT

                  k=c_second_table_u->index - c_first_table_u->index;
                  SaveCode.write(reinterpret_cast<char *>(&k), sizeof(c_first_table_u->index));
                  k1=c_second_table_u->op1- c_first_table_u->op1;
                  SaveCode.write(reinterpret_cast<char *>(&k1), sizeof(c_first_table_u->op1));
                  k2=c_second_table_u->op2- c_first_table_u->op2;
                  SaveCode.write(reinterpret_cast<char *>(&k2), sizeof(c_first_table_u->op2));
#ifdef PRINT_OUT
                  if((c_first_table_u->type == 1) && (middle))
                    cout << ">u[" << k << "]=u[" << k1 << "]*u[" << k2 << "]\n"; //"] (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n";
                  else if(c_first_table_u->type == 2)
                    cout << ">u[" << k << "]+=u[" << k1 << "]*u[" << k2 << "]\n"; //"] (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n";
                  else
                    cout << ">u[" << k << "]=1/(1-u[" << k1 << "]*u[" << k2 << "])\n"; //"]) (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n";
#endif //PRINT_OUT
                }
#endif /**SAVE**/
              break;
            case 5:
              nb_table_u++;
#ifdef SAVE
              if((i == 1) && (middle))
                {
                  nb_middle_save_table_u++;
                  SaveCode.write(reinterpret_cast<char *>(&c_first_table_u->type), sizeof(c_first_table_u->type));
                  SaveCode.write(reinterpret_cast<char *>(&c_first_table_u->index), sizeof(c_first_table_u->index));
#ifdef PRINT_OUT
                  cout << ">push(u[" << c_first_table_u->index << "])\n"; // "]) (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n";
#endif //PRINT_OUT
                  k=c_second_table_u->index - c_first_table_u->index;
                  SaveCode.write(reinterpret_cast<char *>(&k), sizeof(c_first_table_u->index));
#ifdef PRINT_OUT
                  cout << ">push(u[" << k << "])\n"; // "]) (" << tmp_table_u << ", " << tmp_table_u->pNext << ")\n";
#endif //PRINT_OUT

                }
#endif /**SAVE**/
              break;
            }
          if((c_first_table_u->pNext != second_u_blck->pNext) && (c_first_table_u->pNext != NULL) && (c_second_table_u->pNext != NULL))
            {
              tmp_table_u=c_first_table_u ;
              c_first_table_u = c_first_table_u->pNext;
              free(tmp_table_u);
              tmp_table_u=c_second_table_u ;
              c_second_table_u = c_second_table_u->pNext;
              free(tmp_table_u);
            }
          else
            {
              if(c_first_table_u->pNext != second_u_blck->pNext) /*||(second_u_blck->pNext!=NULL))*/
                {
                  cout << "c_first_table_u->pNext=" << c_first_table_u->pNext << " second_u_blck->pNext=" << second_u_blck->pNext << "\n";
                  cout << "Error not synchronize graph interpolation\n";
                  exit( -1);
                }
              OK = 0;
            }
          op_count++;
        }
      SaveCode.flush();
      k=SaveCode.tellp();
      SaveCode.seekp(cur_pos);
      SaveCode.write(reinterpret_cast<char *>(&nb_table_u), sizeof(nb_table_u));
      SaveCode.seekp(k);
      if(middle)
        {
          first_y_blck = c_first_y_blck;
          second_y_blck = c_second_y_blck;
#ifdef SAVE
          if((i == 1) && (middle))
            {
              middle_save_table_y = (t_table_y*)malloc((nb_endo) * sizeof(t_table_y));
              middle_save_i_table_y = (t_table_y*)malloc((nb_endo) * sizeof(t_table_y));
              nb_middle_save_table_y = nb_endo;
            }
#endif /**SAVE**/
          for(j = 0;j < nb_endo;j++)
            {
#ifdef SAVE
              if(i == 1)
                {
                  if(middle)
                    {
                      middle_save_i_table_y[j].u_index = (int*)malloc(table_y[first_y_blck].nb * sizeof(int));
                      middle_save_i_table_y[j].y_index = (int*)malloc(table_y[first_y_blck].nb * sizeof(int));
                      middle_save_i_table_y[j].index = table_y[second_y_blck].index - table_y[first_y_blck].index;
                      middle_save_i_table_y[j].nb = table_y[first_y_blck].nb;
                      for(k = 0;k < table_y[first_y_blck].nb;k++)
                        {
                          middle_save_i_table_y[j].u_index[k] = table_y[second_y_blck].u_index[k] - table_y[first_y_blck].u_index[k];
                          middle_save_i_table_y[j].y_index[k] = table_y[second_y_blck].y_index[k] - table_y[first_y_blck].y_index[k];
                        }
                      middle_save_table_y[j].u_index = (int*)malloc(table_y[first_y_blck].nb * sizeof(int));
                      middle_save_table_y[j].y_index = (int*)malloc(table_y[first_y_blck].nb * sizeof(int));
                      /*middle_save_table_y[j].y_lead_lag = (int*)malloc(table_y[first_y_blck].nb * sizeof(int));*/
                      middle_save_table_y[j].index = table_y[first_y_blck].index;
                      middle_save_table_y[j].nb = table_y[first_y_blck].nb;
                      for(k = 0;k < table_y[first_y_blck].nb;k++)
                        {
                          middle_save_table_y[j].u_index[k] = table_y[first_y_blck].u_index[k];
                          middle_save_table_y[j].y_index[k] = table_y[first_y_blck].y_index[k];
                        }
                    }
                  else
                    {
                      first_save_i_table_y[j].u_index = (int*)malloc(table_y[first_y_blck].nb * sizeof(int));
                      first_save_i_table_y[j].y_index = (int*)malloc(table_y[first_y_blck].nb * sizeof(int));
                      first_save_i_table_y[j].index = table_y[second_y_blck].index - table_y[first_y_blck].index;
                      first_save_i_table_y[j].nb = table_y[first_y_blck].nb;
                      for(k = 0;k < table_y[first_y_blck].nb;k++)
                        {
                          first_save_i_table_y[j].u_index[k] = table_y[second_y_blck].u_index[k] - table_y[first_y_blck].u_index[k];
                          first_save_i_table_y[j].y_index[k] = table_y[second_y_blck].y_index[k] - table_y[first_y_blck].y_index[k];
                        }
                      first_save_table_y[j].u_index = (int*)malloc(table_y[first_y_blck].nb * sizeof(int));
                      first_save_table_y[j].y_index = (int*)malloc(table_y[first_y_blck].nb * sizeof(int));
                      first_save_table_y[j].index = table_y[first_y_blck].index;
                      first_save_table_y[j].nb = table_y[first_y_blck].nb;
                      for(k = 0;k < table_y[first_y_blck].nb;k++)
                        {
                          first_save_table_y[j].u_index[k] = table_y[first_y_blck].u_index[k];
                          first_save_table_y[j].y_index[k] = table_y[first_y_blck].y_index[k];
                        }
                    }
                }
#endif /**SAVE**/
              table_y[vertex_count].nb = table_y[first_y_blck].nb;
              table_y[vertex_count].index = table_y[first_y_blck].index + (i - 1) * (table_y[second_y_blck].index - table_y[first_y_blck].index);
#ifdef PRINT_OUT
              cout << ">y[" << table_y[vertex_count].index << "]=";
#endif /**PRINT_OUT**/
              for(k = 0;k < table_y[vertex_count].nb;k++)
                {
                  table_y[vertex_count].u_index[k] = table_y[first_y_blck].u_index[k] + (i - 1) * (table_y[second_y_blck].u_index[k] - table_y[first_y_blck].u_index[k]);
                  table_y[vertex_count].y_index[k] = table_y[first_y_blck].y_index[k] + (i - 1) * (table_y[second_y_blck].y_index[k] - table_y[first_y_blck].y_index[k]);
#ifdef PRINT_OUT
                  cout << "u[" << table_y[vertex_count].u_index[k] << "]*y[" << table_y[vertex_count].y_index[k] << "]";
                  if(k + 1 == table_y[vertex_count].nb)
                    cout << "\n";
                  else
                    cout << "+";
#endif /**PRINT_OUT**/
                }
              vertex_count++;
              first_y_blck++;
              second_y_blck++;
            }
        }
    }
#ifdef PRINT_OUT
  cout << "+++ u_count=" << u_count << "\n";
  cout << "out interpolation\n";
#endif

  return (old_table_u);
}


//==================================================================================
bool
SymbolicGaussElimination::Loop_Elimination(t_model_graph* model_graph)
{
  int i, j, pos, i1, i_per;
  bool there_is_a_loop, go_on = true;
  /* destroy the loop*/
  int per = 0;
  there_is_a_loop = 0;
  i_per = 0;
  for(i = 0;i < model_graph->nb_vertices;i++)
    {
      if(nstacked)
        {
          if(model_graph->vertex[i].nb_in_degree_edges > 0)
            {
              if(!(i_per % nb_endo))
                {
                  per++;
                  if(per == 1)
                    {
                      first_u_blck = table_u;
#ifdef PRINT_OUT
                      cout << "first_u_blck=" << first_u_blck << "\n";
#endif
                    }
                  if(per == 2)
                    {
                      pos_nb_first_save_table_u = nb_table_u;
                      second_u_blck = table_u;
#ifdef PRINT_OUT
                      cout << "second_u_blck=" << second_u_blck << "\n";
#endif
                    }
                  if(per == 3)
                    {
                      go_on = false;
                      third_u_blck = table_u;
#ifdef PRINT_OUT
                      cout << "third_u_blck\n";
#endif
                    }
                }
              i_per++;
            }
        }
      pos = -1;
      for(j = 0;((j < model_graph->vertex[i].nb_in_degree_edges) && (pos < 0));j++)
        if(model_graph->vertex[i].in_degree_edge[j].index == i)
          {
            // at least one loop
            pos = j;
          }
      if(pos >= 0)
        {
          there_is_a_loop = 1;
#ifdef DEBUGR
          cout << "--------------------------------------------------------------- \n";
          cout << "Loop elimination on vertex " << model_graph->vertex[i].index << "\n";
#endif /**DEBUGR**/
          if(nstacked)
            {
              /*Store the u_count associated to the loop to re-use it in case of
                loop creation during vertex supression step*/
              loop_table_u_count[i] = model_graph->vertex[i].in_degree_edge[pos].u_count;
              loop_table_vertex_index[i] = model_graph->vertex[i].index;
            }
          if(go_on)
            {
              nb_table_u++;
              table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u) - sizeof(int));
              table_u = table_u->pNext;
              table_u->pNext = NULL;
              table_u->type = 3;
              table_u->index = model_graph->vertex[i].in_degree_edge[pos].u_count;
              table_u->op1 = model_graph->vertex[i].in_degree_edge[pos].u_count;
#ifdef PRINT_OUT
              cout << "u[" << table_u->op1 << "]=1/(-u[" << table_u->op1 << "]) " << table_u << "\n"; // "])  (nb_free_u_list : " << nb_free_u_list  << " u_count=" << u_count << ")\n";
#endif /**PRINT_OUT**/
              for(j = 0;j < model_graph->vertex[i].nb_in_degree_edges;j++)
                {
                  if(j != pos)
                    {
                      i1 = model_graph->vertex[i].in_degree_edge[j].index;
#ifdef DIRECT_COMPUTE
                      nb_table_u++;
                      table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u) - sizeof(int));
                      table_u = table_u->pNext;
                      table_u->pNext = NULL;
                      table_u->type = 7;
                      table_u->index = model_graph->vertex[i].in_degree_edge[j].u_count;
                      table_u->op1 = model_graph->vertex[i].in_degree_edge[pos].u_count;
#ifdef PRINT_OUT
                      cout << "u[" << model_graph->vertex[i].in_degree_edge[j].u_count << "]*=u[" << table_u->op1 << "] " << table_u << "\n"; // "] (nb_free_u_list : " << nb_free_u_list  << " u_count=" << u_count << ")\n";
#endif /**PRINT_OUT**/
#endif /**DIRECT_COMPUTE**/
                    }
                }
            }
          /*elimination of the loop*/
          /*in the out_degree*/
          pos = -1;
          for(j = 0;((j < model_graph->vertex[i].nb_out_degree_edges) && (pos < 0));j++)
            if(model_graph->vertex[i].out_degree_edge[j].index == i)
              pos = j;
          if(pos < 0)
            {
              cout << "Error: not symetric on a loop on vertex " << model_graph->vertex[i].index << "\n";
              print_Graph(model_graph);
              system("pause");
              exit( -1);
            }
          for(j = pos + 1;j < model_graph->vertex[i].nb_out_degree_edges;j++)
            {
              model_graph->vertex[i].out_degree_edge[j - 1].index = model_graph->vertex[i].out_degree_edge[j].index;
              model_graph->vertex[i].out_degree_edge[j - 1].u_count = model_graph->vertex[i].out_degree_edge[j].u_count;
            }
          model_graph->vertex[i].nb_out_degree_edges--;
          /*in the in_degree*/
          pos = -1;
          for(j = 0;((j < model_graph->vertex[i].nb_in_degree_edges) && (pos < 0));j++)
            if(model_graph->vertex[i].in_degree_edge[j].index == i)
              pos = j;
          if(pos < 0)
            {
              cout << "Error: not symetric on a loop on vertex " << model_graph->vertex[i].index << "\n";
              print_Graph(model_graph);
              system("pause");
              exit( -1);
            }
          for(j = pos + 1;j < model_graph->vertex[i].nb_in_degree_edges;j++)
            {
              model_graph->vertex[i].in_degree_edge[j - 1].index = model_graph->vertex[i].in_degree_edge[j].index;
              model_graph->vertex[i].in_degree_edge[j - 1].u_count = model_graph->vertex[i].in_degree_edge[j].u_count;
            }
          model_graph->vertex[i].nb_in_degree_edges--;
        }
    }
  return (there_is_a_loop);
}

//==================================================================================
bool
SymbolicGaussElimination::Vertex_Elimination(t_model_graph* model_graph, int pos, bool* interpolate, int length_markowitz, bool dynamic)
{
  int i, j, k, min_edge, vertex_to_eliminate = -1, to_add, to_add_index, curr_u_count;
  int i1, i2, j1, j2, i3, i4, i5, size;
  int i_max, j_max;
  bool OK;
  t_vertex* lvertex = model_graph->vertex;
  t_vertex* ilvertex;
  int nb_new = 0, nb_free = 0;
  int a_loop = 0;
  char usi;
#ifdef DIRECT_COMPUTE
  int in_y = 0;
#endif /**DIRECT_COMPUTE**/
  size = model_graph->nb_vertices;
  min_edge = size * size + 1;
  vertex_to_eliminate = -1;
  // minimize the product of in and out degree of the remaining vertices (Markowitz criteria)
  for(i = starting_vertex;i < starting_vertex + length_markowitz;i++)
    if(((lvertex[i].nb_out_degree_edges*lvertex[i].nb_in_degree_edges) < min_edge) && (lvertex[i].nb_in_degree_edges > 0))
      {
        min_edge = lvertex[i].nb_out_degree_edges * lvertex[i].nb_in_degree_edges;
        vertex_to_eliminate = i;
      }
  OK = 0;
  //cout << "save_direct=" << save_direct << "\n";
  if(vertex_to_eliminate >= 0)
    {
#ifdef PRINT_OUT_1
      cout << "--------------------------------------------------------------- \n";
      cout << "Elimination of vertex " << lvertex[vertex_to_eliminate].index << " interpolate=" << *interpolate << "\n";
      cout << "min_edge=" << min_edge << " length_markowitz=" << length_markowitz << "\n";
      print_Graph(model_graph);
#endif /**PRINT_OUT**/
#ifdef DIRECT_COMPUTE
      table_y[vertex_count].index = lvertex[vertex_to_eliminate].index;
#ifdef PRINT_OUT
      cout << "vertex_count=" << vertex_count << " size=" << size << " vertex_to_eliminate=" << vertex_to_eliminate << "\n";
#endif
#endif /**DIRECT_COMPUTE**/
      ilvertex = &(lvertex[vertex_to_eliminate]);
      i_max = ilvertex->nb_in_degree_edges;
      j_max = ilvertex->nb_out_degree_edges;
      for(i = 0;i < i_max;i++)
        {
          i1 = ilvertex->in_degree_edge[i].index;
          i2 = ilvertex->in_degree_edge[i].u_count;
#ifdef DIRECT_COMPUTE
          table_y[vertex_count].u_index[in_y] = i2;
          table_y[vertex_count].y_index[in_y] = lvertex[i1].index;
          in_y++;
          nb_table_u++;
          if(save_direct)
            {
              usi=5;
              //SaveCode << (&usi) << i2;
              SaveCode.write(reinterpret_cast<char *>(&usi), sizeof(usi));
              SaveCode.write(reinterpret_cast<char *>(&i2), sizeof(i2));
#ifdef PRINT_OUT_1
              cout << "push(u[" << i2 << "]) \n";
#endif
            }
          else
            {
              table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u) - 2 * sizeof(int));
              table_u = table_u->pNext;
              table_u->pNext = NULL;
              table_u->type = 5;
              table_u->index = i2;
            }
          //cout << "ok0\n";
#ifdef PRINT_OUT_1
          cout << "push(u[" << i2 << "]) (" << table_u << ")\n";
#endif /**PRINT_OUT**/
#endif /**DIRECT_COMPUTE**/
          for(j = 0;j < j_max;j++)
            {
              to_add = -9999999;
              j1 = ilvertex->out_degree_edge[j].index;
              j2 = ilvertex->out_degree_edge[j].u_count;
              //cout << "ok1\n";
              /* Is there already an edge from vertex i1 to vertex j1? */
              for(k = 0;((k < lvertex[i1].nb_out_degree_edges) && (to_add < 0));k++)
                if(lvertex[i1].out_degree_edge[k].index == j1)
                  {
                    /*yes*/
                    to_add = lvertex[i1].out_degree_edge[k].u_count;
                    to_add_index = k;
                  }
              if(to_add != -9999999)
                {
#ifdef DEBUGR
                  cout << " modification of edge between vertices " << lvertex[i1].index << " and " << lvertex[j1].index << "\n";
#endif /**DEBUGR**/
#ifdef DIRECT_COMPUTE
                  if(save_direct)
                    {
                      nb_table_u++;
                      usi=2;
                      SaveCode.write(reinterpret_cast<char *>(&usi), sizeof(usi));
                      SaveCode.write(reinterpret_cast<char *>(&to_add), sizeof(to_add));
                      SaveCode.write(reinterpret_cast<char *>(&j2), sizeof(j2));
                      SaveCode.write(reinterpret_cast<char *>(&i2), sizeof(i2));
                      //SaveCode << (unsigned short int)(2) << to_add  << j2 << i2;
#ifdef PRINT_OUT_1
                      cout << "u[" << to_add << "]+=u[" << j2 << "]*u[" << i2 << "]; \n";
#endif
                    }
                  else
                    {
                      nb_table_u++;
                      table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u));
                      table_u = table_u->pNext;
                      table_u->pNext = NULL;
                      table_u->type = 2;
                      table_u->index = to_add;
                      table_u->op1 = j2;
                      table_u->op2 = i2;
                    }
#ifdef PRINT_OUT_1
                  cout << "u[" << to_add << "]+=u[" << j2 << "]*u[" << i2 << "]; (" << table_u << ")\n";
#endif /**PRINT_OUT**/
#endif /**DIRECT_COMPUTE**/
                }
              else
                {
                  /*now modifie the in_degree edge from vertex_to_elminate to j1 */
                  if(j1 == i1)
                    {
                      /* its a loop  */
                      /*Store it to operate after*/
                      if(a_loop >= size)
                        {
                          cout << "Error : a_loop (" << a_loop << ") >= " << size << "\n";
                          exit( -1);
                        }
                      s_j2[a_loop] = j2;
                      s_i2[a_loop] = i2;
                      s_i1[a_loop] = i1;
                      a_loop++;
                    }
                  else
                    {
#ifdef DEBUGR
                      cout << " creation of edge between vertices " << lvertex[i1].index << " and " << lvertex[j1].index << "\n";
#endif /**DEBUGR**/
#ifdef SYMPLIFY  /*no reuse of free numbers, to be sure that there is no cyclcial behaviour in numbering*/
                      curr_u_count = get_free_u_list(dynamic);
#else
                      curr_u_count = u_count;
                      u_count++;
#endif
#ifdef DIRECT_COMPUTE
                      nb_table_u++;
                      nb_new++;
                      if(save_direct)
                        {
                          usi=1;
                          SaveCode.write(reinterpret_cast<char *>(&usi), sizeof(usi));
                          SaveCode.write(reinterpret_cast<char *>(&curr_u_count), sizeof(curr_u_count));
                          SaveCode.write(reinterpret_cast<char *>(&j2), sizeof(j2));
                          SaveCode.write(reinterpret_cast<char *>(&i2), sizeof(i2));
#ifdef PRINT_OUT_1
                          cout << "u[" << curr_u_count << "]=u[" << j2 << "]*u[" << i2 << "]; \n";
#endif
                        }
                      else
                        {
                          table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u));
                          table_u = table_u->pNext;
                          table_u->pNext = NULL;
                          table_u->type = 1;
                          table_u->index = curr_u_count;
                          table_u->op1 = j2;
                          table_u->op2 = i2;
                        }
#ifdef PRINT_OUT_1
                      cout << "u[" << curr_u_count << "]=u[" << j2 << "]*u[" << i2 << "]; (" << table_u << ")\n";
#endif /**PRINT_OUT**/
#endif /**DIRECT_COMPUTE**/
                      /*its not a loop so adding a new adge*/
                      /*modify the in_degree edge from vertex_to_elminate to j1*/
#ifdef SORTED
                      k = lvertex[j1].nb_in_degree_edges;
                      while((i1 < lvertex[j1].in_degree_edge[k - 1].index) && (k > 0))
                        {
                          lvertex[j1].in_degree_edge[k].index = lvertex[j1].in_degree_edge[k - 1].index;
                          lvertex[j1].in_degree_edge[k].u_count = lvertex[j1].in_degree_edge[k - 1].u_count;
                          k--;
                        }
                      lvertex[j1].in_degree_edge[k].index = i1;
                      lvertex[j1].in_degree_edge[k].u_count = curr_u_count;
                      (lvertex[j1].nb_in_degree_edges)++;
                      /*Tested*/
                      if(lvertex[j1].max_nb_in_degree_edges<=lvertex[j1].nb_in_degree_edges)
                        {
                          lvertex[j1].max_nb_in_degree_edges=lvertex[j1].nb_in_degree_edges;
                          t_edge *tmp_edge;
                          tmp_edge=lvertex[j1].in_degree_edge;
                          int taille=lvertex[j1].nb_in_degree_edges*sizeof(t_edge);
                          tmp_edge = (t_edge*)realloc(tmp_edge, taille);
                          lvertex[j1].in_degree_edge=tmp_edge;
                        }
                      /*EndTested*/
#else
                      k = lvertex[j1].nb_in_degree_edges;
                      lvertex[j1].in_degree_edge[k].index = i1;
                      lvertex[j1].in_degree_edge[k].u_count = curr_u_count;
                      (lvertex[j1].nb_in_degree_edges)++;
#endif
                      /*now modify the out_degree edge from i1 to vertex_to_elminate */
#ifdef SORTED
                      k = lvertex[i1].nb_out_degree_edges;
                      while((j1 < lvertex[i1].out_degree_edge[k - 1].index) && (k > 0))
                        {
                          lvertex[i1].out_degree_edge[k].index = lvertex[i1].out_degree_edge[k - 1].index;
                          lvertex[i1].out_degree_edge[k].u_count = lvertex[i1].out_degree_edge[k - 1].u_count;
                          k--;
                        }
                      lvertex[i1].out_degree_edge[k].index = j1;
                      lvertex[i1].out_degree_edge[k].u_count = curr_u_count;
                      (lvertex[i1].nb_out_degree_edges)++;
                      /*Tested*/
                      if(lvertex[j1].max_nb_out_degree_edges<lvertex[j1].nb_out_degree_edges)
                        {
                          lvertex[j1].max_nb_out_degree_edges=lvertex[j1].nb_out_degree_edges;
                          lvertex[j1].out_degree_edge=(t_edge*)realloc(lvertex[j1].out_degree_edge,lvertex[j1].nb_out_degree_edges*sizeof(t_edge));
                        }
                      /*EndTested*/
#else
                      k = lvertex[i1].nb_out_degree_edges;
                      lvertex[i1].out_degree_edge[k].index = j1;
                      lvertex[i1].out_degree_edge[k].u_count = curr_u_count;
                      (lvertex[i1].nb_out_degree_edges)++;
#endif
                    }
                }

            }
        }
#ifdef SIMPLIFY
      nb_free_u_list1 = 0;
#endif /**SIMPLIFY**/
      for(i = 0;i < i_max;i++)
        {
          /*now supress the out_degree edge from i1 to vertex_to_elminate */
          i1 = ilvertex->in_degree_edge[i].index;
#ifdef DEBUGR
          cout << " supress the out_degree edge from " << lvertex[i1].index << " to " << ilvertex->index << "\n";
#endif /**DEBUGR**/
          to_add = -1;
          for(k = 0;((k < lvertex[i1].nb_out_degree_edges) && (to_add < 0));k++)
            if(lvertex[i1].out_degree_edge[k].index == vertex_to_eliminate)
              to_add = k;
          if(to_add < 0)
            {
              cout << "error: Model_Graph not correctly filled in out_degree (edge from " << lvertex[i1].index << " to " << lvertex[vertex_to_eliminate].index << ")\n";
              print_Graph(model_graph);
              system("pause");
              exit( -1);
            }
#ifdef SIMPLIFYS
          nb_free++;
          set_free_u_list(lvertex[i1].out_degree_edge[to_add].u_count);
#endif /**SIMPLIFY**/
#ifdef SIMPLIFYT
          if(simplification_allowed)
            {
              free_u_list1[nb_free_u_list1] = lvertex[i1].out_degree_edge[to_add].u_count;
              if(nb_free_u_list1 > max_nb_free_u_list1)
                max_nb_free_u_list1 = nb_free_u_list1;
              nb_free_u_list1++;
            }
#endif /**SIMPLIFY**/
          for(k = to_add + 1;k < lvertex[i1].nb_out_degree_edges;k++)
            {
              lvertex[i1].out_degree_edge[k - 1].index = lvertex[i1].out_degree_edge[k].index;
              lvertex[i1].out_degree_edge[k - 1].u_count = lvertex[i1].out_degree_edge[k].u_count;
            }
          (lvertex[i1].nb_out_degree_edges)--;
        }
      for(j = 0;j < j_max;j++)
        {
          /*now supress the in_degree edge from vertex_to_elminate to j1 */
          j1 = ilvertex->out_degree_edge[j].index;
#ifdef DEBUGR
          cout << " supress the in_degree edge from " << ilvertex->index << " to " << lvertex[j1].index << "\n";
#endif /**DEBUGR**/
          to_add = -1;
          for(k = 0;((k < lvertex[j1].nb_in_degree_edges) && (to_add < 0));k++)
            if(lvertex[j1].in_degree_edge[k].index == vertex_to_eliminate)
              to_add = k;
          if(to_add < 0)
            {
              cout << "error: Model_Graph not correctly filled in in_degree (edge from " << lvertex[vertex_to_eliminate].index << " to " << lvertex[j1].index << ")\n";
              print_Graph(model_graph);
              system("pause");
              exit( -1);
            }
#ifdef SIMPLIFYS
          set_free_u_list(lvertex[j1].in_degree_edge[to_add].u_count);
          nb_free++;
#endif /**SIMPLIFY**/
#ifdef SIMPLIFYT
          if(simplification_allowed)
            {
              free_u_list1[nb_free_u_list1] = lvertex[j1].in_degree_edge[to_add].u_count;
              if(nb_free_u_list1 > max_nb_free_u_list1)
                max_nb_free_u_list1 = nb_free_u_list1;
              nb_free_u_list1++;
            }
#endif /**SIMPLIFY**/
          for(k = to_add + 1;k < lvertex[j1].nb_in_degree_edges;k++)
            {
              lvertex[j1].in_degree_edge[k - 1].index = lvertex[j1].in_degree_edge[k].index;
              lvertex[j1].in_degree_edge[k - 1].u_count = lvertex[j1].in_degree_edge[k].u_count;
            }
          (lvertex[j1].nb_in_degree_edges)--;
        }
      if(a_loop)
        {
          for(i = 0;i < a_loop;i++)
            {
              i1 = s_i1[i];
              if(nstacked)
                {
                  /*re-use the u_count index eliminated during the loop elmination*/
                  k = 0;
                  while((k <= nb_loop_table) && (loop_table_vertex_index[k] != model_graph->vertex[i1].index))
                    k++;
                  if(loop_table_vertex_index[k] != model_graph->vertex[i1].index)
                    {
                      cout << "Error can't find the loop associated to vertex " << model_graph->vertex[i1].index << "\n";
                      k = 0;
                      while((k <= nb_loop_table) && (loop_table_vertex_index[k] != model_graph->vertex[i1].index))
                        {
                          cout << "loop_table_vertex_index[" << k << "]=" << loop_table_vertex_index[k] << " =? model_graph->vertex[" << i1 << "].index=" << model_graph->vertex[i1].index << "\n";
                          k++;
                        }
                      exit( -1);
                    }
                  curr_u_count = loop_table_u_count[k] ;
                }
              else
                {
#ifdef SIMPLIFY
                  curr_u_count = get_free_u_list(dynamic);
#else /**SIMPLIFY**/
                  curr_u_count = u_count;
                  u_count++;
#endif /**SIMPLIFY**/
                }
              nb_table_u++;
              if(save_direct)
                {
                  usi=6;
                  SaveCode.write(reinterpret_cast<char *>(&usi), sizeof(usi));
                  SaveCode.write(reinterpret_cast<char *>(&curr_u_count), sizeof(curr_u_count));
                  SaveCode.write(reinterpret_cast<char *>(&s_i2[i]), sizeof(s_i2[i]));
                  SaveCode.write(reinterpret_cast<char *>(&s_j2[i]), sizeof(s_j2[i]));
#ifdef PRINT_OUT_1
                  cout << "u[" << curr_u_count << "]=1/(1-u[" << s_i2[i] << "]*u[" << s_j2[i] << "]) ;\n";
#endif
                }
              else
                {
                  table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u) );
                  table_u = table_u->pNext;
                  table_u->pNext = NULL;
                  table_u->type = 6;
                  table_u->index = curr_u_count;
                  table_u->op1 = s_i2[i];
                  table_u->op2 = s_j2[i];
                }
              nb_new++;
#ifdef PRINT_OUT_1
              cout << "u[" << curr_u_count << "]=1/(1-u[" << s_i2[i] << "]*u[" << s_j2[i] << "]) (" << table_u << ");\n";
#endif /**PRINT_OUT**/
              i5 = curr_u_count;
              OK = 0;
              for(k = 0;k < lvertex[i1].nb_in_degree_edges;k++)
                {
                  i3 = lvertex[i1].in_degree_edge[k].index;
                  if(i3 != i1)
                    {
                      i4 = lvertex[i1].in_degree_edge[k].u_count;
#ifdef DIRECT_COMPUTE
                      nb_table_u++;
                      if(save_direct)
                        {
                          usi=7;
                          SaveCode.write(reinterpret_cast<char *>(&usi), sizeof(usi));
                          SaveCode.write(reinterpret_cast<char *>(&i4), sizeof(i4));
                          SaveCode.write(reinterpret_cast<char *>(&i5), sizeof(i5));
#ifdef PRINT_OUT_1
                          cout << "u[" << i4 << "]*=u[" << i5 << "]; \n";
#endif
                        }
                      else
                        {
                          table_u->pNext = (t_table_u*)malloc(sizeof(t_table_u) - sizeof(int));
                          table_u = table_u->pNext;
                          table_u->pNext = NULL;
                          table_u->type = 7;
                          table_u->index = i4;
                          table_u->op1 = i5;
                        }
#ifdef PRINT_OUT_1
                      cout << "u[" << i4 << "]*=u[" << i5 << "]; (" << table_u << ")\n";
#endif /**PRINT_OUT**/
#endif /**DIRECT_COMPUTE**/
                    }
                }
#ifdef SIMPLIFY
              if(simplification_allowed)
                {
                  set_free_u_list(i5);
                  nb_free++;
                }
#endif /**SIMPLIFY**/
            }
        }
#ifdef SIMPLIFYS
      if(simplification_allowed)
        for(i = 0;i < nb_free_u_list1;i++)
          {
            set_free_u_list(free_u_list1[i]);
            nb_free++;
          }
#endif /**SIMPLIFY**/
      /*Change all indexes*/
      free(ilvertex->in_degree_edge);
      free(ilvertex->out_degree_edge);
      for(i = vertex_to_eliminate + 1;i < size;i++)
        lvertex[i - 1] = lvertex[i];
      (model_graph->nb_vertices)--;
      OK = 0;
      for(i = 0;i < model_graph->nb_vertices;i++)
        {
          for(j = 0;j < lvertex[i].nb_in_degree_edges;j++)
            {
              OK = 1;
              if(lvertex[i].in_degree_edge[j].index > vertex_to_eliminate)
                (lvertex[i].in_degree_edge[j].index)--;
            }
          for(j = 0;j < lvertex[i].nb_out_degree_edges;j++)
            {
              if(lvertex[i].out_degree_edge[j].index > vertex_to_eliminate)
                (lvertex[i].out_degree_edge[j].index)--;
            }
        }
      if(in_y>max_nb_table_y)
        max_nb_table_y=in_y;
#ifdef DIRECT_COMPUTE
      table_y[vertex_count].nb = in_y;
#ifdef PRINT_OUT_1
      cout << "y[" << table_y[vertex_count].index << "]=";
      for(i = 0;i < table_y[vertex_count].nb;i++)
        {
          cout << "u[" << table_y[vertex_count].u_index[i] << "]*y[" << table_y[vertex_count].y_index[i] << "]";
          if(i + 1 < table_y[vertex_count].nb)
            cout << "+";
          else
            cout << "\n";
        }
#endif /**PRINT_OUT**/
      vertex_count++;
#endif /**DIRECT_COMPUTE**/
    }
#ifdef PRINT_OUT
  cout << "nb_new=" << nb_new << " nb_free=" << nb_free << "\n";
  cout <<"end of vertex elimination\n";
#endif
  return (OK);
}


//==================================================================================
void
SymbolicGaussElimination::Gaussian_Elimination(t_model_graph* model_graph
#ifdef SAVE
                                               , string file_name
#endif
                                               , bool dynamic)
{
  int i, size, j, k, per;
  bool OK, interpolate = false;
  t_table_u *First_prologue_table_u;
  int length_markowitz=0, per_next = 3;
  bool try_to_interpolate=false;
  int cur_pos=0;
  int prologue_nb_table_u=0, first_nb_prologue_save_table_y=0;
  int nb_first_u_blck, nb_second_u_blck=0, nb_third_u_blck=0;
  int nb_first_y_blck, nb_second_y_blck=0, nb_third_y_blck=0;

  size = model_graph->nb_vertices;
#ifdef PRINT_OUT
  cout << "going to open file file_open=" << file_open << " file_name={" << file_name << "}\n";
#endif
  if(file_open)
    SaveCode.open((file_name + ".bin").c_str(), ios::out | ios::in | ios::binary | ios ::ate );
  else
    SaveCode.open((file_name + ".bin").c_str(), ios::out | ios::binary);
  file_open = true;
  if(!SaveCode.is_open())
    {
      cout << "Error : Can't open file \"" << file_name << ".bin\" for writing\n";
      exit( -1);
    }
#ifdef PRINT_OUT
  print_Graph(model_graph);
#endif /**PRINT_OUT**/
  s_i1 = (int*)malloc(size * sizeof(int));
  s_i2 = (int*)malloc(size * sizeof(int));
  s_j2 = (int*)malloc(size * sizeof(int));
#ifdef SIMPLIFY
  simplification_allowed = 1;
  MAX_FREE_U_LIST = size * /*size *//*12*/50;
  free_u_list = (int*)malloc(MAX_FREE_U_LIST* sizeof(int));
#ifdef PRINT_OUT
  cout << "free_u_list=" << free_u_list << " length=" << ceil((double)4 / (double)2*(double)u_count)*sizeof(int) << "\n";
  cout << "free_u_list length0=" << (*((short int*)&(*(((char*)free_u_list) - 4))) & ~(3)) << "\n";
  cout << "free_u_list1=(int*)malloc(" << u_count << ");\n";
#endif
  free_u_list1 = (int*)malloc(u_count * sizeof(int));
#ifdef PRINT_OUT
  cout << "after free alloc u_count=" << u_count << "\n";
#endif
#endif /**SIMPLIFY**/
#ifdef PRINT_OUT
  cout << "periods=" << periods << "\n";
#endif
  if(dynamic)
    u_count = (periods + y_kmin + y_kmax) * u_count / (y_kmin + y_kmax + 2);
  if(nstacked)
    {
      nb_loop_table = model_graph->nb_vertices;
      loop_table_u_count = (int*)malloc(nb_loop_table * sizeof(int));
      loop_table_vertex_index = (int*)malloc(nb_loop_table * sizeof(int));
    }
#ifdef PRINT_OUT
  cout << "after nb_loop_table=" << nb_loop_table  << "\n";
  system("pause");
#endif

  SaveCode.write(reinterpret_cast<char *>(&nb_endo), sizeof(nb_endo));
  SaveCode.write(reinterpret_cast<char *>(&u_count), sizeof(u_count));
  SaveCode.write(reinterpret_cast<char *>(&u_count_init), sizeof(u_count_init));


#ifdef PRINT_OUT
  cout << "u_count=" << u_count << "\n";
  //We start with the elimination of all loops in the graph
  cout << "going to loop elminate\n";
#endif
  save_direct=false;
  if(nstacked)
    save_direct=false;
  else
    {
      i=0;
      SaveCode.write(reinterpret_cast<char *>(&i), sizeof(i));
      SaveCode.write(reinterpret_cast<char *>(&i), sizeof(i));
      SaveCode.write(reinterpret_cast<char *>(&i), sizeof(i));
      SaveCode.write(reinterpret_cast<char *>(&i), sizeof(i));
    }
  if(Loop_Elimination(model_graph))
    {
#ifdef PRINT_OUT
      cout << "->table_u=" << table_u << "\n";
#endif
      if(nstacked)
        {
          //if there are some loops eliminated, we have to update it for the next blocks
          // So we call the interpolation procedure
          //cout << "nb_first_table_u=";
          interpolation(model_graph, table_y,     1, false, 0);
        }
#ifdef PRINT_OUT
      else
        {
          cout << "nb_first_save_table_u=" << nb_first_save_table_u << " nb_table_u=" << nb_table_u << "\n";
        }
#endif
#ifdef PRINT_OUT
      cout << "->table_u=" << table_u << "\n";
#endif
    }
  if(nstacked)
    {
      first_nb_prologue_save_table_y = vertex_count;
      nb_prologue_save_table_y=0;
      First_prologue_table_u = table_u;
      cur_pos=SaveCode.tellp();
      prologue_nb_table_u=nb_table_u;
      i=0;
      SaveCode.write(reinterpret_cast<char *>(&i), sizeof(i));
#ifdef PRINT_OUT
      cout << "--> fixe prologue\n";
#endif
    }
#ifdef PRINT_OUT
#ifdef SIMPLIFY
  cout << "1 free_u_list length0=" << (*((short int*)&(*(((char*)free_u_list) - 4))) & ~(3)) << "\n";
#endif /**SIMPLIFY**/
#endif
#ifdef DEBUGR
  //Check the graph
#ifdef PRINT_OUT
  cout << "going to Check Graph\n";
#endif
  Check_Graph(model_graph);
#endif
  per = 0;
  OK = true;
#ifdef PRINT_OUT
  cout << "nb_endo=" << nb_endo << " y_kmin=" << y_kmin << " y_kmax=" << y_kmax << " size=" << size << "\n";
#endif
  int pos = 0;
  j = 0;
  if(nstacked)
    try_to_interpolate = true;
#ifdef PRINT_OUT
  cout << "try_to_interpolate=" << try_to_interpolate << "\n";
  cout << "y_kmin=" << y_kmin << " y_kmax=" << y_kmax << "\n";
  cout << "(0) nb_table_u=" << nb_table_u << "\n";
#endif
  if(nstacked)
    save_direct=true;
  while(OK)
    {
#ifdef DEBUGR
      Check_Graph(model_graph);
#endif /**DEBUGR**/
      if(!(j % nb_endo))
        length_markowitz = nb_endo;
      else
        length_markowitz--;
      if(nstacked)
        {
#ifdef PRINT_OUT
          cout << "try_to_interpolate=" << try_to_interpolate << "\n";
#endif
          if(try_to_interpolate)
            {
#ifdef PRINT_OUT
              cout << j << " % " << nb_endo << " = " << (j % nb_endo) << "\n";
#endif
              if(!(j % nb_endo))
                {
                  per++;
                  if(per == y_kmin + 1)
                    {
                      nb_prologue_save_table_y=vertex_count-first_nb_prologue_save_table_y;
                      prologue_save_table_y=(t_table_y*)malloc(nb_prologue_save_table_y*sizeof(*prologue_save_table_y));
                      for(i=first_nb_prologue_save_table_y;i<vertex_count;i++)
                        prologue_save_table_y[i-first_nb_prologue_save_table_y]=table_y[i];
                      k=SaveCode.tellp();
                      SaveCode.seekp(cur_pos);
                      i=nb_table_u-prologue_nb_table_u;
                      SaveCode.write(reinterpret_cast<char *>(&i), sizeof(i));
                      SaveCode.seekp(k);
                      if(nstacked)
                        save_direct=false;
                      first_u_blck = table_u;
                      nb_first_u_blck= nb_table_u;
                      nb_first_y_blck= vertex_count;
                      first_y_blck = vertex_count;
#ifdef PRINT_OUT
                      cout << "(1) nb_table_u=" << nb_table_u << "\n";
                      cout << "first_u_blck=" << first_u_blck << "\n";
                      system("pause");
#endif
                    }
                  else if(per == y_kmin + 2)
                    {
                      second_u_blck = table_u;
                      nb_second_u_blck= nb_table_u;
                      second_y_blck = vertex_count;
                      nb_second_y_blck= vertex_count;
#ifdef PRINT_OUT
                      cout << "(2) nb_table_u=" << nb_table_u << "\n";
                      cout << "second_u_blck=" << second_u_blck << "\n";
                      system("pause");
#endif
                    }
                  else if(per >= y_kmin + 3)
                    {
                      third_u_blck = table_u;
                      nb_third_u_blck= nb_table_u;
                      third_y_blck = vertex_count;
                      nb_third_y_blck= vertex_count;
#ifdef PRINT_OUT
                      cout << "(3) nb_table_u=" << nb_table_u << "\n";
                      cout << "third_u_blck=" << third_u_blck << "\n";
                      system("pause");
#endif
                    }
                }
            }
          OK = Vertex_Elimination(model_graph, pos, &interpolate, length_markowitz, dynamic);
        }
      else
        {
#ifdef PRINT_OUT
          cout << "j=" << j << "\n";
#endif
          OK = Vertex_Elimination(model_graph, pos, &interpolate, length_markowitz, dynamic);
        }
#ifdef SIMPLIFY
#endif /**SIMPLIFY**/
      j++;
#ifdef SIMPLIFY
      if(!(j % nb_endo))
        nb_free_u_list = 0;
#endif /**SIMPLIFY**/
      if(nstacked)
        {
#ifdef PRINT_OUT
          cout << "try_to_interpolate=" << try_to_interpolate << "\n";
#endif
          if(try_to_interpolate)
            {
#ifdef PRINT_OUT
              cout << j << " % " << nb_endo << " = " << (j % nb_endo) << "\n";
#endif
              if(!((j) % nb_endo))
                {
#ifdef PRINT_OUT
                  cout << "per (" << per << ") >= y_kmin (" << y_kmin << ")+3\n";
#endif
                  if( per >= y_kmin + 3)
                    {
                      // pctimer_t t1, t2;
#ifdef PRINT_OUT
                      cout << "bef check regularity\n";
                      system("pause");
#endif
                      if(Check_Regularity(first_u_blck, second_u_blck, third_u_blck))
                        {
#ifdef PRINT_OUT
                          cout << "af check regularity OK \n";
                          system("pause");
                          cout << "nb_first_save_table_u=" << nb_first_save_table_u << "\n";
#endif
#ifdef PRINT_OUT
                          cout << "table_u=" << table_u << "  table_u->pNext=" << table_u->pNext << "\n";
                          cout << "(3) nb_table_u=" << nb_table_u << "\n";
#endif
                          // t1 = pctimer();

                          k=SaveCode.tellp();
#ifdef PRINT_OUT
                          cout << "middle_count_loop= " << middle_count_loop << "\n";
#endif
                          i=0;
                          SaveCode.write(reinterpret_cast<char *>(&i), sizeof(i));
                          table_u = interpolation(model_graph, table_y, per_next /*+1*/ + 2, true, per - y_kmin - 3);
#ifdef PRINT_OUT
                          cout << "middle_count_loop= " << middle_count_loop << "\n";
#endif
                          i=SaveCode.tellp();
                          SaveCode.seekp(k);
                          SaveCode.write(reinterpret_cast<char *>(&middle_count_loop), sizeof(middle_count_loop));
                          SaveCode.seekp(i);
                          i=0;
                          SaveCode.write(reinterpret_cast<char *>(&i), sizeof(i));  //last_u
                          OK = false;
                          vertex_count -= nb_endo;
                          // t2 = pctimer();
#ifdef PRINT_OUT
                          // cout << "interpolate=" << 1000*(t2 - t1) << "\n";
                          cout << "table_u=" << table_u << "  table_u->pNext=" << table_u->pNext << "\n";
#endif
#ifdef SAVE
                          last_save_table_u = table_u;
                          nb_last_save_table_y = vertex_count;
#endif /**SAVE**/
                        }
                      else
                        {
                          nb_prologue_save_table_y=nb_second_y_blck-first_nb_prologue_save_table_y;
                          prologue_save_table_y=(t_table_y*)malloc(nb_prologue_save_table_y*sizeof(*prologue_save_table_y));
                          for(i=first_nb_prologue_save_table_y;i<nb_second_y_blck;i++)
                            prologue_save_table_y[i-first_nb_prologue_save_table_y]=table_y[i];
                          k=SaveCode.tellp();
                          SaveCode.seekp(cur_pos);
                          i=nb_second_u_blck-prologue_nb_table_u+1;
                          //cout << "nb_prologue_save_table_u=" << i << "\n";
                          SaveCode.write(reinterpret_cast<char *>(&i), sizeof(i));
                          SaveCode.seekp(k);
                          write_to_file_table_u_b(first_u_blck,second_u_blck, &nb_prologue_save_table_u,true);
                          //cout << "nb_prologue_save_table_u(1)=" << nb_prologue_save_table_u << "\n";
                          nb_first_u_blck=nb_second_u_blck;
                          nb_second_u_blck=nb_third_u_blck+1;
                          nb_first_y_blck=nb_second_y_blck;
                          nb_second_y_blck=nb_third_y_blck;
                          first_u_blck = second_u_blck;
                          second_u_blck = third_u_blck;
                          first_y_blck = second_y_blck;
                          second_y_blck = third_y_blck;
                        }
                    }
                }
#ifdef DEBUGR
              cout << 100*j / size << "%\n";
#endif /**DEBUGR**/
            }
        }
    }
#ifdef PRINT_OUT
  cout << "nstacked=" << nstacked << "\n";
  cout << "nb_first_save_table_y=" << nb_first_save_table_y << "\n";
  cout << "nb_prologue_save_table_y=" << nb_prologue_save_table_y << "\n";
  cout << "nb_middle_save_table_y=" << nb_middle_save_table_y << "\n";
  cout << "nb_last_save_table_y=" << nb_last_save_table_y << "\n";
#endif
  if((nstacked)&&(!nb_last_save_table_y))
    {
      cout << "not synchronized per=" << per << "\n";
      exit(-1);
    }

#ifdef PRINT_OUT
  cout << "end vertex supress\n";
#endif
#ifdef SAVE
  if(!nstacked)
    table_u->pNext = NULL;
  if(nstacked)
    {
      table_u = NULL;
      last_save_table_u = NULL;
      nb_last_save_table_y = 0;
      OK = true;
    }
  else
    {
      SaveCode.write(reinterpret_cast<char *>(&nb_table_u), sizeof(nb_table_u));
      write_to_file_table_u_b(First_table_u->pNext, table_u->pNext, &nb_last_save_table_u, false );
      nb_last_save_table_y = vertex_count;
      last_save_table_y=(t_table_y*)malloc(nb_last_save_table_y*sizeof(*last_save_table_y));
      for(i=0;i<nb_last_save_table_y;i++)
        last_save_table_y[i]=table_y[i];
    }
#ifdef PRINT_OUT
  cout << "goint to write nb_endo=" << nb_endo << "\n";
#endif
#ifdef PRINT_OUT
  cout << "-->write prologue\n";
#endif
#ifdef PRINT_OUT
  cout << "-->write first\n";
#endif
#ifdef PRINT_OUT
  cout << "middle_count_loop=" << middle_count_loop << "\n";
  cout << "-->write middle\n";
#endif
#ifdef PRINT_OUT
  cout << "-->write last\n";
#endif
  nb_last_save_table_u = i;
#ifdef PRINT_OUT
  cout << "nb_prologue_save_table_y=" << nb_prologue_save_table_y << "\n";
#endif
  write_to_file_table_y( prologue_save_table_y, NULL, nb_prologue_save_table_y, 0);
#ifdef PRINT_OUT
  cout << "nb_first_save_table_y=" << nb_first_save_table_y << "\n";
#endif
  write_to_file_table_y( first_save_table_y, NULL, nb_first_save_table_y, 0);
#ifdef PRINT_OUT
  cout << "nb_middle_save_table_y=" << nb_first_save_table_y << "\n";
#endif
  write_to_file_table_y( middle_save_table_y, middle_save_i_table_y, nb_middle_save_table_y, nb_middle_save_table_y);
#ifdef PRINT_OUT
  cout << "//// last_save_table_y ===\n";
  cout << "nb_last_save_table_y=" << nb_last_save_table_y << "\n";
#endif
  write_to_file_table_y( last_save_table_y, NULL, nb_last_save_table_y, 0);
  SaveCode.close();
#endif /**SAVE**/
#ifdef SIMPLIFY
#ifdef PRINT_OUT
  cout << "try_to_interpolate=" << try_to_interpolate << "\n";
#endif
  free(free_u_list);
  free(free_u_list1);
#endif /**SIMPLIFY**/
  free(s_i1);
  free(s_i2);
  free(s_j2);
}

void
SymbolicGaussElimination::file_is_open1()
{
  file_open=true;
}

void
SymbolicGaussElimination::file_is_open()
{
  file_open=true;
  //file_is_open1();
}

void
SymbolicGaussElimination::SGE_compute(Model_Block *ModelBlock, int blck, bool dynamic, string file_name, int endo_nbr)
{
  t_model_graph *model_graph;
  int block_u_count, nb_table_y;
  // pctimer_t t1, t2;
  int i;
  int mean_var_in_equation;

  init_glb();
  model_graph = (t_model_graph*)malloc(sizeof(*model_graph));
  nstacked = dynamic;
#ifdef PRINT_OUT
  periods = ModelBlock->Periods;
  cout << "nstacked=" << nstacked << "\n";
  // t1 = pctimer();
  cout << "periods=" << periods << "\n";
#endif

  u_count = ModelBlock_Graph(ModelBlock, blck, dynamic, model_graph, endo_nbr, &block_u_count, &starting_vertex, &periods, &nb_table_y, &mean_var_in_equation);
  if(dynamic)
    {
      cout << "Mean endogenous variable per equation: " << mean_var_in_equation << ", density indicator=" << double(mean_var_in_equation)/endo_nbr*100 << "%\n";
      cout << "Coding the model...";
    }
  Time = periods;
#ifdef PRINT_OUT
  // t2 = pctimer();
  // cout << "/* ModelBlock_Graph : " << 1000*(t2 - t1) << " milliseconds u_count : " << u_count << "*/\n";
#endif
  int size = model_graph->nb_vertices;
  nb_endo = ModelBlock->Block_List[blck].Size;
  y_kmin = ModelBlock->Block_List[blck].Max_Lag;
  y_kmax = ModelBlock->Block_List[blck].Max_Lead;
  periods = ModelBlock->Periods;
  if(periods<=3)
    nstacked=false;
  //cout << "periods=" << periods << " y_kmin=" << y_kmin << " y_kmax=" << y_kmax << "\n";
  u_count_init = u_count;
#ifdef PRINT_OUT
  cout << "size : " << size << "\n";
#endif
  table_y = (t_table_y*)malloc((size + 1) * sizeof(t_table_y));
  for(i = 0;i <= size;i++)
    {
      table_y[i].u_index = (int*)malloc(nb_table_y * sizeof(int));
      table_y[i].y_index = (int*)malloc(nb_table_y * sizeof(int));
      table_y[i].index = i;
      table_y[i].nb = 0;
    }
  table_u = (t_table_u*)malloc(sizeof(*table_u)-3*sizeof(int));
  table_u->pNext = NULL;
  First_table_u = table_u;
#ifdef PRINT_OUT
  cout << "oł en est-on ? nb_table_y=" << nb_table_y << "\n";
  system("pause");
  cout << "u_count_init=" << u_count_init << "\n";
  cout << "table_u=" << table_u << "\n";
  t1 = pctimer();
#endif
  Gaussian_Elimination(model_graph, file_name, dynamic);
  free_model_graph(model_graph);
#ifdef PRINT_OUT
  cout << "max_nb_table_y=" << max_nb_table_y << "\n";
  cout << "max_nb_in_degree_edges=" << max_nb_in_degree_edges << "\n";
  cout << "max_nb_out_degree_edges=" << max_nb_out_degree_edges << "\n";
  t2 = pctimer();
  cout << "/* Gaussian Elimination : " << 1000*(t2 - t1) << " milliseconds u_count : " << u_count << "*/\n";
#endif
}