static INT volatile *ctrbtab;
static pthread_mutex_t *mutex_cblk; /** Mutex pour proteger l'acquisition d'un bloc-colonne **/
static pthread_mutex_t *mutex_ctrb; /** Mutex pour proteger la modification du compteur de contribution d'un bloc-colonne **/
static pthread_mutex_t *mutex_block; /** Mutex pour proteger l'ajout dans un blok **/


/**** This function allows the static pivoting of small diagonals entries ****/
int API_CALL(LDLt_piv)(INT n, INT stride, FLOAT *coefftab, FLOAT normA, FLOAT epsilon)
{
  /******************************************************************************************/
  /* This function computes the LDLt    factorization of a dense  matrix stored in a single */
  /* vector by column                                                                       */
  /* On entry :                                                                             */
  /* n : dimension of the matrix                                                            */
  /* stride : stride to pass from a column to another in the vector                         */
  /* normA : the norm of matrix (used for static numerical pivoting)                        */
  /* epsilon : the tolerance that will be used in the accelerator (also used for pivoting)  */
  /* On return :                                                                            */
  /*  The matrix contains in its lower triangular part the L factor and the D factor is     */
  /* stored in the diagonal                                                                 */
  /******************************************************************************************/

  INT i, j;
  FLOAT *colptr, *facolptr;
  FLOAT piv;
  FLOAT f;
  FLOAT thresh;
  int UN = 1;
  int len;
  

 /** Threshold for pivoting **/
  thresh = sqrt(epsilon)*normA;

  for(i=0;i<n-1;i++)
    {
      /** Find the index of the next column  **/
      colptr = coefftab + i*stride;
      
      /** Static pivoting **/
      if(fabs(colptr[i]) < thresh)
	{
	  fprintf(stderr, "Pivot \n");
	  colptr[i] = thresh;
	}

      piv = 1.0/colptr[i];

      /** Update lower triangular part L **/
      for(j=i+1;j<n;j++)
	{
	  f = -colptr[j]*piv;
	  facolptr = coefftab + stride * j;
	  len = n-j;
	  SOPALIN_AXPY(len, f, colptr+j, UN, facolptr+j, UN);
	  colptr[j] = -f;
	}
    }

  return 0;
}  

 
void* API_CALL(factorisation_LDLt_thread)(INT me)
{
  /*************************************************************/
  /* This function performs in place the L.Lt factorization of */
  /* a symmetric sparse matrix structured in dense blocks of   */
  /* coefficients                                              */
  /*************************************************************/
  SymbolMatrix *symbmtx;
  SolverMatrix *solvmtx;
  INT p, q, k, m;
  FLOAT *W, *F; /** buffers de travail **/
  INT UN=1;
  FLOAT alpha, beta;
  INT stride;
  INT cdim, rdim, hdim, mdim;
  FLOAT *cc; 
  FLOAT *bc; 
  FLOAT *fc; 
  FLOAT *wc; 
  INT facestride;
  INT bloknum, facecblknum;
  INT decalcol;
  Clock timer;  /*** Un timer pour mesurer le temps de factorisation **/

  solvmtx = datacode;
  symbmtx = &(solvmtx->symbmtx);


 

  
  printf("Debut de la factorisation pour thread %d \n", me);
  clockInit(&timer);
  clockStart(&timer);
  


  /*** Allocation de buffers de travail ***/
  W = (FLOAT *)malloc(sizeof(FLOAT)*solvmtx->coefmax);
  assert(W != NULL);
  F = (FLOAT *)malloc(sizeof(FLOAT)*solvmtx->coefmax);
  assert(F != NULL);
  
  for(k=0;k<symbmtx->cblknbr;k++)
    {
      if(pthread_mutex_trylock(mutex_cblk+k) != 0)
	{
	  /** Le bloc colonne est deja pris en charge par une autre thread : on passe au suivant **/
	  /*fprintf(stderr, "Pas possible prendre ctrb %d \n" ,k);*/
	  continue;
	}
      
      /*** Attente active tant que toutes les contributions n'ont pas ete recu ***/
      while(ctrbtab[k] != 0);


      cdim   = symbmtx->cblktab[k].lcolnum - symbmtx->cblktab[k].fcolnum + 1; /** Largeur du bloc colonne **/
      stride = solvmtx->cblktab[k].stride;
   


      /**********************************/
      /* Factorisation du bloc diagonal */
      /**********************************/
      /** Calcul du pointeur dans coeftab vers le  début des coefficients du bloc diagonal **/ 
      cc     = solvmtx->coeftab + solvmtx->bloktab[ symbmtx->cblktab[k].bloknum ].coefind;
      API_CALL(LDLt_piv)(cdim, stride, cc, 0.0, 0.0);
      

      if(k == symbmtx->cblknbr)
	continue;  /** Pas de bloc extradiagonal le calcul est terminer **/

      /************************************************************/
      /* "Diviser" les blocs extra diagonaux par le bloc diagonal */
      /************************************************************/
      {
	/** On effectue l'operation : M = M.(L^-1)t  avec M = l'ensemble compacte des blocs extra-diagonaux du bloc colonne **/
	char *side = "R";
	char *uplo = "L";
	char *trans = "T";
	char *diag = "U";

	alpha = 1.0;
	
	/** Calcul du pointeur dans coeftab vers le  début des coefficients du premier bloc extra diagonal **/ 
	bloknum = symbmtx->cblktab[k].bloknum+1; /** Indice du premier bloc extra-diagonal **/
	bc = solvmtx->coeftab + solvmtx->bloktab[bloknum].coefind;


	rdim = stride - cdim; /** Hauteur de la surface compactée des bloc extra-diagonaux **/
	SOPALIN_TRSM(side,  uplo, trans, diag, rdim, cdim, alpha, cc, stride, bc, stride);
	
	/** Calcul de F dans un buffer temporaire: on divise toutes les colonnes par leur terme diagonal  **/
	for(m=0;m<cdim;m++)
	  {
	    /** Pointeur vers la m-ieme colonne dans le bloc colonne **/
	    bc = solvmtx->coeftab + solvmtx->bloktab[bloknum].coefind + m*stride;

	    /** Pointeur vers la m-ieme colonne dans le buffer F **/
	    fc = F + m*stride;
	    alpha = 1.0/ cc[m*stride + m]; /** alpha est le m-ieme terme diagonal du bloc diagonal **/
	    SOPALIN_COPY(rdim, bc, UN, fc, UN); /** On copie les bloc extra diagonaux dans F **/
	    SOPALIN_SCAL(rdim, alpha, fc, UN); /** On divise la m-ieme colonne par le terme diagonal **/
	  }
      
      }
    
      /** Boucle sur les blocs extra-diagonaux **/
      rdim = stride - cdim;  /** rdim est la hauteur des bloc extra-diagonaux >= p dans le bloc colonne k **/
                             /** On l'initialise à la heuteur totale des bloc extra diagonaux **/

      for(p=symbmtx->cblktab[k].bloknum+1;p<symbmtx->cblktab[k+1].bloknum;p++)
	{

	  /*** Stocker dans W le produit des blocs {A(q,k), q>=p} avec Ft(p,k) ***/
	  char *opA = "N"; /** Ne pas transposer A(**) **/
	  char *opF = "T"; /** Transpose F **/

	  /** Calcul du pointeur dans coeftab vers le début des coefficients du bloc extra-diagonal p **/ 
	  bc = solvmtx->coeftab + solvmtx->bloktab[p].coefind;

	  /** Calcul du pointeur dans F vers le debut des coefficients du bloc extra-diagonal p correspondant **/
	  bloknum = symbmtx->cblktab[k].bloknum + 1; /** Indice du premier bloc extra-diagonal**/
	  fc = F + solvmtx->bloktab[p].coefind -  solvmtx->bloktab[ bloknum ].coefind;

	  /** Calcul de la hauteur du bloc extra diagonal p **/
	  hdim = (symbmtx->bloktab[p].lrownum - symbmtx->bloktab[p].frownum+1);
	  alpha = 1.0;
	  beta = 0.0;
	  
	  SOPALIN_GEMM(opA, opF, rdim, hdim, cdim, alpha, bc, stride, fc, stride, beta, W, stride);
	  
	  /****************************************************************************************/
	  /** Pour tout les blocs extra-diagonaux {A(j,k), j>= i} faire A(i,j) = A(i,j) - W(i, ) **/
	  /****************************************************************************************/
	  facecblknum = symbmtx->bloktab[p].cblknum; /** Indice du bloc colonne en face du bloc extra-diagonal p **/

	  /** Stride du cblk en face du bloc extra-diagonal p **/
	  facestride = solvmtx->cblktab[facecblknum].stride;

	  /** Nombre de colonne "à gauche" de la zone modifiée dans le bloc colonne facecblknum **/
	  decalcol = symbmtx->bloktab[p].frownum - symbmtx->cblktab[facecblknum].fcolnum;

	  /*** Largueur de la zone modifiee == hauteur du bloc p dans le bloc colonne k **/
	  mdim = symbmtx->bloktab[p].lrownum - symbmtx->bloktab[p].frownum+1;

	  bloknum = symbmtx->cblktab[facecblknum].bloknum; /** Indice du bloc diagonal du bloc colonne facebloknum **/

	  for(q=p;q<symbmtx->cblktab[k+1].bloknum;q++) /** Pour tous les bloc extradiagonaux A(q,k) avec q >= p **/
	    {

	      /** On trouve le bloc A(bloknum,facebloknum) en face du bloc extradiagonal A(q,k) **/
	      while(! (symbmtx->bloktab[bloknum].frownum <= symbmtx->bloktab[q].frownum
		       && symbmtx->bloktab[bloknum].lrownum >= symbmtx->bloktab[q].lrownum) )
		{
		  bloknum++;
		  assert(bloknum < symbmtx->cblktab[facecblknum+1].bloknum);
		}

	      /** Calcul du pointeur de debut de la zone modifiée dans le bloc A(bloknum, facebloknum) **/
	      bc = solvmtx->coeftab + solvmtx->bloktab[bloknum].coefind + decalcol*facestride + symbmtx->bloktab[q].frownum - symbmtx->bloktab[bloknum].frownum;
	      
	      /** Calcul du pointeur de debut du bloc correspondant à A(q,k) dans W **/
	      wc = W + solvmtx->bloktab[q].coefind - solvmtx->bloktab[p].coefind;
	      
	      /*** Update de la contribution ***/

	      pthread_mutex_lock(mutex_block+bloknum);
	      hdim = symbmtx->bloktab[q].lrownum - symbmtx->bloktab[q].frownum+1; /** Hauter du bloc colonne q **/
	      alpha = -1.0;
	      for(m=0; m < mdim; m++)
		{
		  SOPALIN_AXPY(hdim, alpha, wc, UN, bc, UN);
		  wc += stride;
		  bc += facestride;
		}
	      pthread_mutex_unlock(mutex_block+bloknum); 
	    }
	  
	  /***********************************************************************************/
	  /*** On decremente le compteur de contribution du bloc-colonne en face du bloc p ***/
	  /***********************************************************************************/
	  pthread_mutex_lock(mutex_ctrb+facecblknum);
	  ctrbtab[facecblknum]--;
	  assert(ctrbtab[facecblknum]>=0);
	  pthread_mutex_unlock(mutex_ctrb+facecblknum);



	  /** Mettre à jour la hauteur des bloc extra diagonaux qui restent **/
	  rdim -= symbmtx->bloktab[p].lrownum - symbmtx->bloktab[p].frownum+1;

	}


      /** On recopie F dans les bloc extradiagonaux **/
      rdim = stride - cdim;  /* hauteur des blocs extradiagonaux **/
      bloknum = symbmtx->cblktab[k].bloknum+1; /** Indice du premier bloc extra-diagonal **/
      for(m=0;m<cdim;m++)
	{
	  /** Pointeur vers la m-ieme colonne dans le bloc colonne **/
	  bc = solvmtx->coeftab + solvmtx->bloktab[bloknum].coefind + m*stride;
	  
	  /** Pointeur vers la m-ieme colonne dans le buffer F **/
	  fc = F + m*stride;
	  SOPALIN_COPY(rdim, fc, UN, bc, UN); /** On copie la colonne de F **/
	}

    }

  fprintf(stdout, "Thread %ld a fini la factorisation \n", (long)me);

  /** Libération des buffers **/
  free(W);
  free(F);




  clockStop(&timer);
  fprintf(stdout, "Temps de la factorisation %g \n", clockVal(&timer));

  /*** Resolution des systèmes triangulaires ****/
  API_CALL(up_down_smp)((void *)me);
  API_CALL(gmres_smp)((void *)me);


  return NULL;
}

void* API_CALL(solveur_thread)(SolverMatrix *solvmtx)
{
  INT i, k;
  int ret,status;
  Backup b;
  void *arguments;
  SymbolMatrix *symbmtx;
  
  arguments = (void *)malloc(sizeof(void *)*2);
  symbmtx = &(datacode->symbmtx);
  /******* ALLOCATIONS DES STRUCTURES POUR PARALLELISER LA FACTO *******/
  
  mutex_cblk = (pthread_mutex_t *)malloc(sizeof(pthread_mutex_t)*symbmtx->cblknbr);
  mutex_ctrb = (pthread_mutex_t *)malloc(sizeof(pthread_mutex_t)*symbmtx->cblknbr);
  mutex_block = (pthread_mutex_t *)malloc(sizeof(pthread_mutex_t)*symbmtx->bloknbr);
  ctrbtab = (INT *)malloc(sizeof(INT) * symbmtx->cblknbr);
  bzero((INT *)ctrbtab, sizeof(INT)*symbmtx->cblknbr);
  for(k=0;k<symbmtx->cblknbr-1;k++)
    {
      INT facecblknum;
      for(i=symbmtx->cblktab[k].bloknum+1;i<symbmtx->cblktab[k+1].bloknum;i++)
	ctrbtab[symbmtx->bloktab[i].cblknum]++;
    }

  /******* Initialisation des mutex ******/
  for(i=0; i < symbmtx->cblknbr; i++)
    {
      pthread_mutex_init(&mutex_cblk[i],NULL);
      pthread_mutex_init(&mutex_ctrb[i],NULL);
    }

  for(i=0; i < symbmtx->bloknbr; i++)
    pthread_mutex_init(&mutex_block[i],NULL);


  /****** CREATION DES THREADS ********/
  for (i=0;i<SOLV_THRDNBR;i++)
    {
      ret = pthread_create(&calltab[i],NULL,(void *(*)(void *) )API_CALL(factorisation_LDLt_thread), (void *)i);
      if (ret) {printf("error in thread create\n"); EXIT(MOD_SOPALIN,THREAD_ERR);}
    }


  /****** ATTENTE FIN ET DESTRUCTION DES THREADS **********/
  for (i=0;i<SOLV_THRDNBR;i++)
    {
      ret = pthread_join(calltab[i],(void**)&status);
      fprintf(stdout, "La thread %ld a fini \n", (long)i);
      if (ret) {printf("error in thread join\n"); EXIT(MOD_SOPALIN,THREAD_ERR);}      
    }


  /******* Destruction des mutex ******/
  for(i=0; i < symbmtx->cblknbr; i++)
    {
      pthread_mutex_destroy(&mutex_cblk[i]);
      pthread_mutex_destroy(&mutex_ctrb[i]);
    }

  for(i=0; i < symbmtx->bloknbr; i++)
    pthread_mutex_destroy(&mutex_block[i]);

  free(mutex_cblk);
  free(mutex_ctrb);
  free(mutex_block);
  free(ctrbtab);
}