toeplitz.c

/*
** file: toeplitz.c, version 1.0beta, May 2012  
**
** Project:  Midapack library, ANR MIDAS'09
** Purpose:  Provide Toeplitz algebra tools suitable for Cosmic Microwave Background (CMB)
**           data analysis.
**
** Authors:  Pierre Cargemel, Frederic Dauvergne, Maude Le Jeune, Antoine Rogier, 
**           Radek Stompor (APC, Paris)
**
** 
** Copyright (c) 2010-2012 APC CNRS Université Paris Diderot
** 
** 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 http://www.gnu.org/licenses/gpl.html
**
** For more information about ANR MIDAS'09 project see
** http://www.apc.univ-paris7.fr/APC_CS/Recherche/Adamis/MIDAS09/index.html
** 
** ACKNOWLEDGMENT: This work has been supported in part by the French National 
** Research Agency (ANR) through COSINUS program (project MIDAS no. ANR-09-COSI-009).
**
**
** Revision 1.0  2012/05/01?  fd@APC
** Official release 1.0, including only the Toeplitz algebra module.
**
*/

#include "toeplitz.h"


//=========================================================================
//Global parameters


int NFFT = 1;


int FFTW_FLAG  = FFTW_ESTIMATE;

//Global parameter just to know the rank for printing when VERBOSE mode is on
int PRINT_RANK = -1;


//=========================================================================


int print_error_message(int error_number, char const *file, int line)
{
  char *str_mess;
  str_mess = (char *) malloc(100 * sizeof(char));
  if(error_number == 1)    
    sprintf (str_mess, "Error on line %d of %s. Toeplitz band width > vector size\n", line, file);
  if(error_number == 2)    
    sprintf (str_mess, "Error on line %d of %s. Bad allocation.\n", line, file);
  if(error_number == 3)    
    sprintf (str_mess, "Error on line %d of %s. Error at fftw multithread initialization.\n", line, file);
  if(error_number == 7)
    sprintf (str_mess, "Error on line %d of %s.\n", line, file);
  fprintf(stderr, "%s", str_mess);
  printf("%s", str_mess);
  return error_number;
  
}


//=========================================================================


int optimal_blocksize(int n, int lambda, int bs_flag) {
  int bs;  //computed optimal block size
  int min_bs;  //minimum block size used for the computation
  int min_pow2;  //minimum power of two index used for the block size computation

  min_bs = 3*lambda;
  min_pow2 = (int) ceil( log(min_bs)/log(2) );
  bs = pow(2, min_pow2);


  if (bs_flag==1) {  //flag to use a different formula to find an optimal block size
  while(bs < 2*(lambda-1)*log(bs+1) && bs<n) {
    bs = bs*2;  }}

  if (bs > n)       //This is to avoid block size much bigger than the matrix. Append mostly
    bs = 3*lambda;  //when the matrix is small compared to his bandwith

  if(VERBOSE)
    printf("Computed optimal blocksize is %d (with lambda = %d)\n", bs, lambda);

  return bs;
}


//=========================================================================


int tpltz_init(int n, int lambda, int *nfft, int *blocksize, fftw_complex **T_fft, double *T, fftw_complex **V_fft, double **V_rfft, fftw_plan *plan_f, fftw_plan *plan_b)
{
  //initialize block size
  *blocksize = optimal_blocksize(n, lambda, 0);
  //initialize nfft
  *nfft = NFFT;
  //initialize fftw plan allocation flag
  int fftw_flag = FFTW_FLAG;

  //initialize fftw for omp threads
  fftw_init_omp_threads();
  //initialize fftw array and plan for T (and make it circulant first)
  circ_init_fftw(T, (*blocksize), lambda, T_fft);
  //initialize fftw array and plan for V  
  rhs_init_fftw(nfft, (*blocksize), V_fft, V_rfft, plan_f, plan_b, fftw_flag);

  if (VERBOSE)
    printf("Initialization finished successfully\n");

  return 0;
}


//=========================================================================


int fftw_init_omp_threads()
{
  int status;

  //initialize fftw omp threads
#ifdef fftw_MULTITHREADING
  status = fftw_init_threads();
  if (status==0)
    return print_error_message (3, __FILE__, __LINE__);

#pragma omp parallel
  n_thread = omp_get_num_threads();
  if (VERBOSE)
    printf("Using %d threads\n", n_thread);

  fftw_plan_with_nthreads(n_thread);
  if (VERBOSE)
    printf("Using multithreaded FFTs with %d threads\n", n_thread);
#endif

  return 0;
}


//=========================================================================


int rhs_init_fftw(int *nfft, int fft_size, fftw_complex **V_fft, double **V_rfft, fftw_plan *plan_f, fftw_plan *plan_b, int fftw_flag)
{
  //allocate fftw arrays and plans for V
  *V_fft  = (fftw_complex*) fftw_malloc((*nfft)*(fft_size/2+1) * sizeof(fftw_complex) );
  *V_rfft = (double*) fftw_malloc((*nfft)*fft_size * sizeof(double) );
  if (*V_fft==0 || *V_rfft==0)
    return print_error_message (2, __FILE__, __LINE__);

  *plan_f = fftw_plan_many_dft_r2c(1, &fft_size, (*nfft), *V_rfft, &fft_size, 1, fft_size, *V_fft, NULL, 1, fft_size/2+1, fftw_flag );
  *plan_b = fftw_plan_many_dft_c2r(1, &fft_size, (*nfft), *V_fft, NULL, 1, fft_size/2+1, *V_rfft, &fft_size, 1, fft_size, fftw_flag );


  return 0;
}


//=========================================================================


int circ_init_fftw(double *T, int fft_size, int lambda, fftw_complex **T_fft)
{
  //routine variable
  int i;
  int circ_fftw_flag = FFTW_ESTIMATE;
  //allocation for T_fft
  *T_fft = (fftw_complex*) fftw_malloc( (fft_size/2+1) * sizeof(fftw_complex) );
  if (*T_fft==0)
    return print_error_message (2, __FILE__, __LINE__);
  double *T_circ = (double*) (*T_fft);

  //inplace fft
  fftw_plan plan_f_T;
  plan_f_T   = fftw_plan_dft_r2c_1d( fft_size, T_circ, *T_fft, circ_fftw_flag );

  //make T circulant
#pragma omp parallel for 
  for(i=0; i<fft_size+2;i++) 
    T_circ[i] = 0.0;

  T_circ[0] = T[0];
  for(i=1;i<lambda;i++) {
    T_circ[i] = T[i]; 
    T_circ[fft_size-i] = T[i];    } 

  fftw_execute(plan_f_T);
  fftw_destroy_plan(plan_f_T);

  return 0;
}


//=========================================================================


int tpltz_cleanup(fftw_complex **T_fft, fftw_complex **V_fft, double **V_rfft,fftw_plan *plan_f, fftw_plan *plan_b){
  fftw_destroy_plan(*plan_f);
  fftw_destroy_plan(*plan_b);
  fftw_free(*T_fft);
  fftw_free(*V_fft);
  fftw_free(*V_rfft);
#ifdef fftw_MULTITHREADING
  fftw_cleanup_threads();
#endif
  fftw_cleanup();
}


//=========================================================================


int copy_block(int ninrow, int nincol, double *Vin, int noutrow, int noutcol, double *Vout, int inrow, int incol, int nblockrow, int nblockcol, int outrow, int outcol, double norm, int set_zero_flag)
{
  int i, j, p, offsetIn, offsetOut;

  //do some size checks first
  if( (nblockcol > nincol) || (nblockrow > ninrow) || (nblockcol > noutcol) || (nblockrow > noutrow)) {
    printf("Error in routine copy_block. Bad size setup.\n");
    return print_error_message(7, __FILE__, __LINE__);  
  }

  if(set_zero_flag) {
#pragma omp parallel for private(i)
    for(i=0;i<noutrow*noutcol;i++)  //could use maybe memset but how about threading
      Vout[i] = 0.0;  
  }

  offsetIn = ninrow*incol+inrow;
  offsetOut = noutrow*outcol+outrow;

#pragma omp parallel for private(i,j,p)
  for(i=0;i<nblockcol*nblockrow;i++) {    //copy the block
    j = i/nblockrow;
    p = i%nblockrow;
    Vout[offsetOut+j*noutrow+p] = Vin[offsetIn+j*ninrow+p]*norm;
  }

  return 0;
}


//=========================================================================


int scmm_direct(int fft_size, fftw_complex *C_fft, int ncol, double *V_rfft, double **CV, fftw_complex *V_fft, fftw_plan plan_f_V, fftw_plan plan_b_CV)
{
  //routine variables
  int sizeT = fft_size/2+1;
  int i, idx;

  //perform forward FFT
  fftw_execute(plan_f_V); //input in V_rfft; output in V_fft

#pragma omp parallel for private(idx)
  for(i=0;i<ncol*sizeT;i++) {
    idx = i%sizeT;
    V_fft[i][0] = C_fft[idx][0]*V_fft[i][0]-C_fft[idx][1]*V_fft[i][1];
    V_fft[i][1] = C_fft[idx][0]*V_fft[i][1]+C_fft[idx][1]*V_fft[i][0];    }

  //perform  backward FFts
  fftw_execute(plan_b_CV); //input in V_fft; output in V_rfft 

  return 0;
}


//=========================================================================


int scmm_basic(double **V, int blocksize, int m, fftw_complex *C_fft, int lambda, double **CV, fftw_complex *V_fft, double *V_rfft, int nfft, fftw_plan plan_f_V, fftw_plan plan_b_CV)
{
  //routine variables
  int i,k; //loop index
  int nloop = (int) ceil((1.0*m)/nfft);  //number of subblocks

  // Loop over set of columns
  int ncol = min(nfft, m);  //a number of columns to be copied from the data to working matrix
                            //equal the number of simultaneous FFTs

#pragma omp parallel for private(i)
  for( i=0;i<blocksize*ncol;i++)
    V_rfft[i] = 0.0;  //could use maybe memset but how about threading

  for(k=0;k<nloop;k++) {   //this is the main loop over the set of columns
    if (k==nloop-1)   //last loop ncol may be smaller than nfft
      ncol = m-(nloop-1)*nfft;  

  //init fftw matrices. 
  //extracts a block of ncol full-length columns from the data matrix and embeds in a bigger
  //matrix padding each column with lambda zeros. Note that all columns will be zero padded 
  //thanks to the "memset" call above

  copy_block(blocksize, m, (*V), blocksize, ncol, V_rfft, 0, k*nfft, blocksize, ncol, 0, 0, 1.0, 0);
  //note: all nfft vectors are transformed below ALWAYS in a single go (if ncol < nfft) the extra
  //useless work is done. 

  scmm_direct(blocksize, C_fft, ncol, V_rfft, CV, V_fft, plan_f_V, plan_b_CV);
  //note: the parameter CV is not really used

  //extract the relevant part from the result
  copy_block(blocksize, ncol, V_rfft, blocksize, m, (*CV), 0, 0, blocksize, ncol, 0, k*nfft, 1.0/((double) blocksize), 0);

  }  //end of loop over the column-sets


  return 0;
}


//=========================================================================


int stmm_core(double **V, int n, int m, fftw_complex *T_fft, int blocksize, int lambda, fftw_complex *V_fft, double *V_rfft, int nfft, fftw_plan plan_f, fftw_plan plan_b, int flag_offset)
{

  //routine variable 
  int status;
  int i,j,k,p;  //loop index 
  int currentsize;
  int blocksize_eff = blocksize-2*lambda;  //just a good part after removing the overlaps
  int nbloc = ceil( (1.0*n)/blocksize_eff);  //a number of subblock of slide/overlap algorithm

  if (flag_offset==1)
    nbloc = ceil((1.0*(n-2*lambda))/blocksize_eff);


  double *V_bloc, *TV_bloc;
  V_bloc  = (double *) calloc(blocksize*m, sizeof(double));
  TV_bloc = (double *) calloc(blocksize*m, sizeof(double));      
  if((V_bloc==0)||(TV_bloc==0))
    return print_error_message(2, __FILE__, __LINE__);

  int offset=0;
  if (flag_offset==1)
    offset=lambda;

  //note: lambda-1 instead of lambda should be enough to work.
  //The same for mpi routines, we generally use one extra term in the communications

  int iV = 0;  //"-lambda+offset";  //first index in V 
  int iTV = offset;  //first index in TV 

  //"k=0";
  //first subblock separately as it requires some padding. prepare the block of the data vector
  //with the overlaps on both sides
  currentsize = min( blocksize-lambda+offset, n-iV); 
  //note: if flag_offset=0, pad first lambda elements with zeros (for the first subblock only)
  // and if flag_offset=1 there is no padding with zeros.
  copy_block( n, m, *V, blocksize, m, V_bloc, 0, 0, currentsize, m, lambda-offset, 0, 1.0, 0);

  //do block computation 
  status = scmm_basic(&V_bloc, blocksize, m, T_fft, lambda, &TV_bloc, V_fft, V_rfft, nfft, plan_f, plan_b);

  if (status!=0) {
    printf("Error in stmm_core.");
    return print_error_message(7, __FILE__, __LINE__);  }


  //now copy first the new chunk of the data matrix **before** overwriting the input due to overlaps !
  iV = blocksize_eff-lambda+offset;

  if(nbloc > 1) {
    currentsize  = min( blocksize, n-iV);  //not to overshoot          

    int flag_reset = (currentsize!=blocksize);  //with flag_reset=1, always "memset" the block.
    copy_block( n, m, *V, blocksize, m, V_bloc, iV, 0, currentsize, m, 0, 0, 1.0, flag_reset);
  }

  //and now store the ouput back in V
  currentsize  = min( blocksize_eff, n-iTV);       // to trim the extra rows
  copy_block( blocksize, m, TV_bloc, n, m, *V, lambda, 0, currentsize, m, iTV, 0, 1.0, 0);


  iTV += blocksize_eff;
  //now continue with all the other subblocks    
  for(k=1;k<nbloc;k++) {

    //do bloc computation 
    status = scmm_basic(&V_bloc, blocksize, m, T_fft, lambda, &TV_bloc, V_fft, V_rfft, nfft, plan_f, plan_b);
    if (status!=0) break;


    iV += blocksize_eff;
    //copy first the next subblock to process 
    if(k != nbloc-1) {
      currentsize = min(blocksize, n-iV);  //not to overshoot          

      int flag_resetk = (currentsize!=blocksize);  //with flag_reset=1, always "memset" the block.
      copy_block( n, m, *V, blocksize, m, V_bloc, iV, 0, currentsize, m, 0, 0, 1.0, flag_resetk);
    }

    //and then store the output in V 
    currentsize  = min( blocksize_eff, n-iTV);  //not to overshoot               
    copy_block( blocksize, m, TV_bloc, n, m, *V, lambda, 0, currentsize, m, iTV, 0, 1.0, 0);
    iTV += blocksize_eff;

  }//end bloc computation 


  free(V_bloc);
  free(TV_bloc);

  return status;
}



//=========================================================================


int stmm_reshape(double **V, int n, int m, int id0, int l, fftw_complex *T_fft, int lambda, fftw_complex *V_fft, double *V_rfft, fftw_plan plan_f, fftw_plan plan_b, int blocksize, int nfft)
{

  //routine variable 
  int i,j,k,p;  //loop index 
  int idf       = id0+l-1;
  int cfirst    = id0/n;  //first column index 
  int clast     = idf/n;  //last column index 
  int clast_1   = (idf+1)/n;
  int m_eff     = clast - cfirst + 1 ;  //number of columns 
  int rfirst    = id0%n;
  int rlast     = idf%n;

  int fft_size;
  int flag_nfft=1;  //flag to avoid nfft reshaping when nfft=1. By default, it is active,
                    //thus if nfft==1
  int flag_offset;

  if (l<lambda) // test to avoid communications errors
    return print_error_message (1, __FILE__, __LINE__);


  //allocation for the vector shape V1 
  double *V1;
  int v1_size;

  v1_size = l + (m_eff-1)*lambda;  //we have  l = (m_eff-2)*n+ (n-rfirst) + (rlast+1) (just for information)
  V1 = (double *) calloc(v1_size, sizeof(double));


  //reshape the matrix V into a vector shape V1. It is add lambda zeros
  //at the end of each column.
#pragma omp parallel for private(i,j,p)   
  for(i=0;i<l;i++) {
    j = (i+rfirst)/n;
    p = (i+rfirst)%n;
    V1[p-rfirst+j*(n+lambda)] = (*V)[i]; }
  //note: not really need copy if m=1. Will be improve as memory use.


  //shortcut possible for nfft==1 with flag_nfft==1. You need also to put flag_offset=0 in the core routine.
  //to make it work correctly.
  if (flag_nfft==1 && nfft==1) {
    flag_offset=0;
    int status = stmm_core(&V1, v1_size, nfft, T_fft, blocksize, lambda, V_fft, V_rfft, nfft, plan_f, plan_b, flag_offset);
    if (status!=0)
      return status;

  }
  else {  //perform sliding windows algorithm on a matrix shape optimize for nfft 

  fft_size = ceil(1.0*v1_size/nfft)+2*lambda;  //maybe use a better name  "fft_nnew"
  //  int blocksize_eff = blocksize-2*lambda;
  //  int nbloc = ceil(1.0*fft_size/blocksize_eff);  //just for information
  //  int fft_size = nbloc*blocksize_eff; //if we want to have exactly the size correspondance
                                       //(not really usefull).

  //reallocation of V with the computed matrix size above 
  (*V) = realloc((*V), fft_size*nfft * sizeof(double));

  //convert the vector shape into a matrix shape optimize for nfft
  vect2nfftblock(V1, v1_size, (*V), fft_size, nfft, lambda);

  //perform sequential multiplication by small blocks 
  flag_offset=1;
  int status = stmm_core(V, fft_size, nfft, T_fft, blocksize, lambda, V_fft, V_rfft, nfft, plan_f, plan_b, flag_offset);
  if (status!=0)
    return status;

  nfftblock2vect((*V), fft_size, nfft, lambda, V1, v1_size);

  }//end else of  if (flag_nfft==1)

  //reallocation of (*V) to his input size
  *V = realloc(*V, l * sizeof(double));

  //reshape the vector shape V1 into the input matrix V. It is remove the lambda zeros
  //at the end of each column used for the computation.
#pragma omp parallel for private(j,p) 
  for(i=0;i<l;i++) {
    j=(i+rfirst)/n;
    p=(i+rfirst)%n;
    (*V)[i] = V1[p-rfirst+j*(n+lambda)]; }


  free(V1);

//ps: we could consider theses transformations as pointers functions to avoid many copies
//of datas. This will be improved in futurs developments as the memory use.


  return 0;
}


//=========================================================================


int vect2nfftblock(double *V1, int v1_size, double *V2, int fft_size, int nfft, int lambda)
{
  int i,j,p;
  //reshape the vector V1 into a matrix shape (*V) for nfft optimization
  // double *V2; V2 = (double *) calloc(fft_size*nfft, sizeof(double));
  //maybe use V2 instead to avoid reallocation. Anyway, find a way to better optimize memory use.

  //copy the lambda elements on the edges of each column of the reshape matrix
  for (j=0; j<nfft; j++){
    for (i=0; i<lambda; i++) {  //top overlap
      if (j*(fft_size-2*lambda)+i-lambda>=0 && j*(fft_size-2*lambda)+i-lambda<v1_size)
        (V2)[j*fft_size+i] = V1[j*(fft_size-2*lambda)+i-lambda];
      else
        (V2)[j*fft_size+i] = 0.0; }
    for (i=0; i<lambda; i++) {  //bottom overlap
      if ((j+1)*(fft_size-2*lambda)+i<v1_size)
        (V2)[fft_size-lambda+i+(j*fft_size)] = V1[(j+1)*(fft_size-2*lambda)+i];
      else
        (V2)[fft_size-lambda+i+(j*fft_size)] = 0.0;  }
  }

  //copy the middle part of the optimize nfft matrix shape (*V)
#pragma omp parallel for private(j,p) 
  for(i=0;i<v1_size;i++) {
    j = i/(fft_size-2*lambda);
    p = i%(fft_size-2*lambda);
    (V2)[lambda+p+(j*fft_size)] = V1[i]; }  //core copy 

  for(i=v1_size;i<nfft*(fft_size-2*lambda);i++) {
    j = i/(fft_size-2*lambda);
    p = i%(fft_size-2*lambda);
    (V2)[lambda+p+(j*fft_size)] = 0.0; }  //padd zeros to extra terms (nfft*fft_size != v1_size)


  return 0;
}


//=========================================================================


int nfftblock2vect(double *V2, int fft_size, int nfft, int lambda, double *V1, int v1_size)
{
  int i,j,p;

  //reshape the matrix shape (*V) used for nfft optimization into the vector shape V1.
  //It is extract the middle part without the lambda edges elements of each column.
#pragma omp parallel for private(j,p) 
  for(i=0;i<v1_size;i++) {
    j = i/(fft_size-2*lambda);
    p = i%(fft_size-2*lambda);
    V1[i] = V2[lambda+p+(fft_size*j)];   }

  //almost equivalent to
  //copy_block(fft_size, nfft, (*V), (fft_size-2*lambda), nfft, V1, lambda, 0, 0, 0, 1.0, 0);
  //not good if (fft_size-2*lambda)*nfft <> v1_size
  //         i.e. if ceil(1.0*v1_size/nfft) <> v1_size/nfft
  //We need to extend V1 memory to use this

  return 0;
}


//=========================================================================


int stmm(double **V, int n, int m, int id0, int l, fftw_complex *T_fft, int lambda, fftw_complex *V_fft, double *V_rfft, fftw_plan plan_f, fftw_plan plan_b, int blocksize, int nfft)
{

  //routine variable 
  int i,j,k,p;  //loop index 
  int idf       = id0+l-1;
  int cfirst    = id0/n;  //first column index 
  int clast     = idf/n;  //last column index 
  int clast_1   = (idf+1)/n;
  int m_eff     = clast - cfirst + 1 ;  //number of columns 
  int rfirst    = id0%n;
  int rlast     = idf%n;

  if (l<lambda) // test to avoid communications errors
    return print_error_message (1, __FILE__, __LINE__);

//the middle
int nloop_middle;  //change it to number of full column to improve memory
  if (m_eff>2)
    nloop_middle = (m_eff-2)/nfft;  
  else
    nloop_middle = 0;

  int m_middle = nfft*nloop_middle;
  int vmiddle_size = n*m_middle;
  double *Vmiddle;


  if (nloop_middle==0) {
    stmm_reshape(V, n, m, id0, l, T_fft, lambda, V_fft, V_rfft, plan_f, plan_b, blocksize, nfft);

  }
  else { //(nloop_middle>0)          
  
    int offset_middle = n-id0;
    Vmiddle = (*V)+offset_middle;

    int flag_offset=0;
    stmm_core(&Vmiddle, n, m_middle, T_fft, blocksize, lambda, V_fft, V_rfft, nfft, plan_f, plan_b, flag_offset);


//edge  (first+last columns + extra column from the euclidian division)
  int v1_size = min(l,(n-rfirst));  
  int v2_size = l-(v1_size+vmiddle_size);
  if (v2_size<0)
    v2_size=0;

  int vedge_size = v1_size+v2_size;
  double *Vedge;
  Vedge = (double *) calloc(vedge_size, sizeof(double));


#pragma omp parallel for private(i)
  for(i=0;i<v1_size;i++)
    Vedge[i] = (*V)[i];


  if (v2_size>0) {
  int offset_edge = v1_size+vmiddle_size;
#pragma omp parallel for private(i)
  for(i=0;i<v2_size;i++)
    Vedge[v1_size+i] = (*V)[offset_edge+i];
  }

  stmm_reshape(&Vedge, n, m, id0, vedge_size, T_fft, lambda, V_fft, V_rfft, plan_f, plan_b, blocksize, nfft);

//Copy result to V:
#pragma omp parallel for private(i)
  for(i=0;i<v1_size;i++)
    (*V)[i]=Vedge[i];

  if (v2_size>0) {
  int offset_edge2 = v1_size+vmiddle_size;
#pragma omp parallel for private(i)
  for(i=0;i<v2_size;i++)
    (*V)[offset_edge2+i]=Vedge[v1_size+i];
  }

  if (vmiddle_size>0) {
  int offset_middle2 = n-rfirst;
#pragma omp parallel for private(i)
  for(i=0;i<vmiddle_size;i++)
    (*V)[offset_middle2+i]=Vmiddle[i];
  }

  free(Vedge);

  } //End (nloop_middle>0)

//ps: this loop can be improved by considering the number of full column instead.
//We can also improve the memory use and the number of copy.

  return 0;
}


//=========================================================================
 

int mpi_stmm(double **V, int n, int m, int id0, int l, double *T, int lambda, MPI_Comm comm) 
{

  //mpi variables 
  int rank;   //rank process 
  int size;   //number of processes 
  MPI_Status status; 
  MPI_Comm_rank(comm, &rank);     
  MPI_Comm_size(comm, &size);     

 
  //routine variables 
  int i,j,k; // some index 
  int idf = id0+l;  // first index of scattered V for rank "rank + 1"; 
  int cfirst = id0/n;   // first column index 
  int clast  = idf/n; // last column index 
  int clast_r = (idf-1)/n; 
  int m_eff  = clast_r - cfirst + 1 ; 
  double *V1, *Lambda;

  // Mpi communication conditions 
  // Mpi comm is needed when columns are truncated 
  int right = rank + 1; 
  int left  = rank - 1; 
  int v1_size = l + 2*lambda; // size including comm 
  if (rank==0 || cfirst*n==id0){ // no left comm 
    v1_size -= lambda; 
    left = MPI_PROC_NULL;}  
  if (rank==(size-1) || clast*n==idf){ // no right comm 
    v1_size -= lambda; 
    right = MPI_PROC_NULL;} 

  // init data to send 
  Lambda=(double *) malloc(2*lambda * sizeof(double)); 
  if (Lambda==0) 
    return print_error_message (2, __FILE__, __LINE__); 
 
  for(i=0;i<lambda;i++)    { 
    Lambda[i]=(*V)[i]; 
    Lambda[i+lambda]=(*V)[i+l-lambda];    } 
  if (VERBOSE) 
    printf("[rank %d] Left comm with %d | Right comm with %d\n", rank, left, right); 

  //send and receive data 
  MPI_Sendrecv_replace(Lambda, lambda, MPI_DOUBLE, left, MPI_USER_TAG, right, MPI_USER_TAG, comm, &status);  //1st comm 
  MPI_Sendrecv_replace(&(*(Lambda+lambda)), lambda, MPI_DOUBLE, right, MPI_USER_TAG, left, MPI_USER_TAG, comm, &status);  //2nd comm 
  
 
  if (l<lambda)  //After sendrecv to avoid problems of communication for others processors
    return print_error_message (1, __FILE__, __LINE__); 

  //copy received data  
  if(left==MPI_PROC_NULL && right==MPI_PROC_NULL)  // 0--0 : nothing to do
    V1 = *V;
  else if(left==MPI_PROC_NULL) { // 0--1 : realloc
    *V = realloc(*V, v1_size * sizeof(double));
    if(*V == NULL) 
      return print_error_message (2, __FILE__, __LINE__); 
    V1 = *V;  }
  else  // 1--1 or 1--0 : new allocation
    V1  = (double *) malloc(v1_size * sizeof(double)); 
  
  if (left!=MPI_PROC_NULL){ 
    for(i=0;i<lambda;i++) 
      V1[i] = Lambda[i+lambda]; 
    id0 -= lambda;} 
  if (right!=MPI_PROC_NULL){ 
    for(i=0;i<lambda;i++) 
      V1[i+v1_size-lambda] = Lambda[i];      } 
  
  // Copy input matrix V 
  int offset = 0; 
  if (left!=MPI_PROC_NULL){ 
    offset = lambda;
#pragma omp parallel for 
    for(i=offset;i<l+offset;i++) 
      V1[i] = (*V)[i-offset]; }
  
  fftw_complex *V_fft, *T_fft; 
  double *V_rfft; 
  fftw_plan plan_f, plan_b; 


  //Compute matrix product 
  int nfft,  blocksize;
  tpltz_init(v1_size, lambda , &nfft, &blocksize, &T_fft, T, &V_fft, &V_rfft, &plan_f, &plan_b);

  if(VERBOSE)
  printf("[rank %d] Before middle-level call : blocksize=%d, nfft=%d\n", rank, blocksize, nfft);

  stmm(&V1, n, m, id0, v1_size, T_fft, lambda, V_fft, V_rfft, plan_f, plan_b, blocksize, nfft);

  tpltz_cleanup(&T_fft, &V_fft, &V_rfft,&plan_f, &plan_b ); 
  
  // Copy output matrix TV
  offset = 0;
  if (left!=MPI_PROC_NULL)
    offset = lambda;
  if(left==MPI_PROC_NULL && right==MPI_PROC_NULL) // 0--0
    *V=V1;
  else if(left==MPI_PROC_NULL) { // 0--1
    V1 = realloc(V1, l * sizeof(double));
    if(V1 == NULL)
      return print_error_message (2, __FILE__, __LINE__);
    *V = V1;     }
  else { // 1--0 or 1--1
#pragma omp parallel for
    for(i=offset;i<l+offset;i++)
      (*V)[i-offset] = V1[i];    }

  if(left!=MPI_PROC_NULL)
    free(V1);
  
  return 0; 
} 


//=========================================================================


int mpi_stbmm(double **V, int *n, int m, int nrow, double *T, int nb_blocks_local, int nb_blocks_all, int *lambda, int *idv, int idp, int local_V_size, MPI_Comm comm)
{

  //MPI parameters 
  int rank;   //process rank 
  int size;   //process number 

  MPI_Status status;
  MPI_Comm_rank(comm, &rank);
  MPI_Comm_size(comm, &size);

  PRINT_RANK=rank;

  FILE *file;
  file = stdout;

  //neighbour if the process for mpi communication conditions if a block is shared
  int right = rank + 1;
  int left  = rank - 1;

  int i,j,k;  //some indexes 


  //Define the indices for each process
  int idv0, idvn;  //indice of the first and the last block of V for each processes

  int *nnew;
  nnew = (int*) calloc( nb_blocks_local, sizeof(int));
  int idpnew;
  int local_V_size_new;  


  int status_params = get_overlapping_blocks_params( nb_blocks_local, idv, n, local_V_size, nrow, idp, &idpnew, &local_V_size_new, nnew, &idv0, &idvn);

  if (VERBOSE==1)
    printf("status_params=%d\n", status_params);

  if( status_params == 0) {
  free( nnew);
    return(0); // no work to be done
  }


  if (VERBOSE==1) {  //print on screen news parameters definition if VERBOSE
    fprintf(file, "news parameters definition:\n");
    fprintf(file, "[%d] idp=%d ; idpnew=%d", idp, idpnew);
    fprintf(file, "[%d] local_V_size=%d ; local_V_size_new=%d\n", rank, local_V_size, local_V_size_new);
    for(i=0;i<nb_blocks_local;i++)
      fprintf(file, "[%d] n[%d]=%d ; nnew[%d]=%d\n", rank, i, n[i], i, nnew[i]);
  }


  int vShft=idpnew-idp;   //new first element of relevance in V

  //Define the column indices:
  //index of the first and the last column of V for the current process
  int idvm0 = idpnew/nrow;
  int idvmn = (idpnew+local_V_size_new-1)/nrow;
  //number of columns of V for the current process
  int ncol_rank = idvmn-idvm0+1;
  //number of blocks for the current process with possibly repetitions 
  int nb_blocs_rank;

  if(ncol_rank == 1) // Empty process not allowed 
    nb_blocs_rank = idvn - idv0 + 1;
  else
    nb_blocs_rank = (ncol_rank-2)*nb_blocks_local + (nb_blocks_local-idv0) + (idvn+1);  //in this case nb_blocks_local = nblocs_all

  if (VERBOSE==1) {  //print on screen 
  fprintf(file, "[%d] nb_blocs_rank=%d, nb_blocks_local=%d\n", rank, nb_blocs_rank, nb_blocks_local);
  }

  //Define the indices for the first and the last element in each blocks
  int idvp0 = idpnew%nrow-idv[idv0];  //index of the first element of the process in the first block
  int idvpn;  //reverse index of the last element of the process in the last block
              //It's the number of remaining elements needed to fully complete the last block
  idvpn = idv[idvn]+nnew[idvn]-1 - (idpnew+local_V_size_new-1)%nrow ;


  //Define the offsets for the first and last blocks of the process for V1
  int v1_size = local_V_size_new;  //length of the new vector V1
  int offset0=0;
  int offsetn=0;

  if(idvp0 != 0) {
    offset0 = min( idvp0, lambda[idv0]);
    v1_size += offset0; }

  if(idvpn != 0) {
    offsetn = min(idvpn, lambda[idvn]);
    v1_size += offsetn; }
//ps: no need v1_size anymore, I keep it just for information.


//initialization for send and receive data
  //vector which receives V datas
  double *LambdaOut=(double *) calloc(lambda[idv0]+lambda[idvn], sizeof(double));
//  double *LambdaOut=(double *) calloc(toSendLeft+toSendRight, sizeof(double));

  if (LambdaOut==0)
    return print_error_message (2, __FILE__, __LINE__);

  int toSendLeft = 0;

  if( offset0 != 0) {
          toSendLeft = min( idv[idv0]+nnew[idv0]-idpnew%nrow, lambda[idv0]);
      for(i=0;i<toSendLeft;i++)
         LambdaOut[i]=(*V)[i];
  }

  int toSendRight = 0;

  if( offsetn != 0) {
    toSendRight = min( (idpnew+local_V_size)%nrow-idv[idvn], lambda[idvn]);
    for(i=0;i<toSendRight;i++)
      LambdaOut[i+lambda[idv0]]=(*V)[i+local_V_size-toSendRight];
 
  }

  if(rank==0 || offset0==0)
    left = MPI_PROC_NULL;
  if(rank==size-1 || offsetn==0)
    right = MPI_PROC_NULL;

  double *LambdaIn=(double *) calloc(lambda[idv0]+lambda[idvn], sizeof(double));
//  double *LambdaIn=(double *) calloc(offset0+offsetn, sizeof(double));
  if (LambdaIn==0)
         return print_error_message (2, __FILE__, __LINE__);

//send and receive data
  MPI_Sendrecv( LambdaOut, toSendLeft, MPI_DOUBLE, left,  MPI_USER_TAG, (LambdaIn+lambda[idv0]), offsetn, MPI_DOUBLE, right, MPI_USER_TAG, comm, &status);
  MPI_Sendrecv( (LambdaOut+lambda[idv0]), toSendRight, MPI_DOUBLE, right,  MPI_USER_TAG, LambdaIn, offset0, MPI_DOUBLE, left, MPI_USER_TAG, comm, &status);

  free(LambdaOut);


//size of the first and the last bloc for the current process
  int v0rank_size, vnrank_size;
  if (nb_blocs_rank == 1) {  //only one block 
    v0rank_size = ((idpnew+local_V_size_new-1)%nrow +1) - idpnew%nrow + offset0 + offsetn;
    vnrank_size = 0; //just for convenience - no really need it
  }
  else { //more than one block
    v0rank_size = idv[idv0] + nnew[idv0] - idpnew%nrow + offset0;
    vnrank_size = ((idpnew+local_V_size_new-1)%nrow +1) - idv[idvn] + offsetn;
  }


//---------------------------------------
//initialization for the blocks loop

  int idv1=0;     //index of *V1

  int mid;  //local number of column for the current block
  //index of the first element of the process inside the first block
  int offset_id0;
  offset_id0 = idvp0;

//fftw variables
  fftw_complex *V_fft, *T_fft;
  double *V_rfft;
  fftw_plan plan_f, plan_b;
//init local block vector
  double *V1block;
  int lambdaShft;

//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//loop on the blocks inside the process
  int nfft, blocksize;
  int iblock;  //index for the loop on the blocks
//  int loopindex;
  int id; //indice of the current block

  int vblock_size;
  int id0block;

  for(iblock=idv0;iblock<idv0+nb_blocs_rank;iblock++) {
    id = iblock%nb_blocks_local;  //index of current block


  if(nnew[id]>0) { //the block is ok

  mid=1;  //always one!

  for( lambdaShft=k=0; k<id; k++)
    lambdaShft += lambda[k];


//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//first case : First block of the process
  if(iblock==idv0) {
  if(VERBOSE)
    fprintf(file, "[%d] First block...\n", rank);

  vblock_size=v0rank_size;
  id0block=(offset_id0-offset0);

  V1block = (double *) calloc(vblock_size, sizeof(double));


#pragma omp parallel for                        
  for (j=0;j<offset0;j++)
    V1block[j] = LambdaIn[j];

//if (nb_blocs_rank == 1) currentsize_middlepart=vblock_size-offset0-offsetn = local_V_size_new
//else currentsize_middlepart=vblock_size-offset0
  int currentsize_middlepart=min(vblock_size-offset0, local_V_size_new);
#pragma omp parallel for   
  for (j=0;j<currentsize_middlepart;j++)
    V1block[offset0+j] = (*V)[j+vShft];

if (nb_blocs_rank == 1) {
#pragma omp parallel for   
  for (j=0;j<offsetn;j++)
    V1block[vblock_size-offsetn+j] = LambdaIn[j+lambda[idv0]];
}

  //init Toeplitz arrays
  tpltz_init(vblock_size, lambda[id], &nfft, &blocksize, &T_fft, (T+lambdaShft), &V_fft, &V_rfft, &plan_f, &plan_b);
  //Toeplitz computation


  if(VERBOSE)
    fprintf(file, "[%d] Before middle-level call : nfft = %d, blocksize = %d\n", rank, nfft, blocksize);
  stmm(&V1block, nnew[id], mid, id0block, vblock_size, T_fft, lambda[id], V_fft, V_rfft, plan_f, plan_b, blocksize, nfft);
  tpltz_cleanup(&T_fft, &V_fft, &V_rfft, &plan_f, &plan_b);


  int currentsize=min(vblock_size-offset0, local_V_size_new);
#pragma omp parallel for
  for (j=0;j<currentsize;j++)
    (*V)[vShft+j] = V1block[offset0+j];


  free(V1block);

  }//end (First case)


//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//Generic case : Generic block of the process
  else if(iblock!=idv0 && iblock!=idv0+nb_blocs_rank-1) {
  if(VERBOSE)
    fprintf(file, "[%d] generic block...\n");

  vblock_size=nnew[id];
  id0block=0;

  V1block = (double *) calloc(vblock_size, sizeof(double));

  idv1 = idv[id]-idp%nrow - vShft + offset0 +nrow*( (iblock/nb_blocks_local) );

#pragma omp parallel for                        
  for (j=0;j<vblock_size;j++)
    V1block[j] = (*V)[j+idv1-offset0+vShft];

  //init Toeplitz arrays
  tpltz_init(nnew[id], lambda[id], &nfft, &blocksize, &T_fft, (T+lambdaShft), &V_fft, &V_rfft, &plan_f, &plan_b);
  //Toeplitz computation
  if(VERBOSE)
    fprintf(file, "[%d] Before middle-level call : nfft = %d, blocksize = %d\n", rank, nfft, blocksize);
  stmm(&V1block, nnew[id], mid, id0block, vblock_size, T_fft, lambda[id], V_fft, V_rfft, plan_f, plan_b, blocksize, nfft);
  tpltz_cleanup(&T_fft, &V_fft, &V_rfft, &plan_f, &plan_b);

#pragma omp parallel for                        
  for (j=0;j<vblock_size;j++)
    (*V)[j+idv[id] - idp%nrow + nrow*( (iblock/nb_blocks_local) )] = V1block[j];


  free(V1block);

  }  //end (Generic case)

//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
// Last case : Last block of the process
  else if(iblock==idv0+nb_blocs_rank-1 && iblock!= idv0) {
  if(VERBOSE)  
    fprintf(file, "[%d] last block...\n");

  vblock_size=vnrank_size;
  id0block=0;

  V1block = (double *) calloc(vblock_size, sizeof(double));

  idv1 = idv[id] - idp%nrow - vShft + offset0  + nrow*( (iblock/nb_blocks_local) );

#pragma omp parallel for   
  for (j=0;j<vblock_size-offsetn;j++)
    V1block[j] = (*V)[j+idv1-offset0+vShft];
#pragma omp parallel for   
  for (j=0;j<offsetn;j++)
    V1block[vblock_size-offsetn+j] = LambdaIn[j+lambda[idv0]];


  //init Toeplitz arrays
  tpltz_init(vblock_size, lambda[id], &nfft, &blocksize, &T_fft, (T+lambdaShft), &V_fft, &V_rfft, &plan_f, &plan_b);
  //Toeplitz computation
  if(VERBOSE)
    printf("[%d] Before middle-level call : nfft = %d, blocksize = %d\n", rank, nfft, blocksize);
  stmm(&V1block, nnew[id], mid, id0block, vblock_size, T_fft, lambda[id], V_fft, V_rfft, plan_f, plan_b, blocksize, nfft);
  tpltz_cleanup(&T_fft, &V_fft, &V_rfft, &plan_f, &plan_b);


#pragma omp parallel for                                                 
  for (j=0;j<vnrank_size-offsetn;j++)
    (*V)[idv[id]-idp%nrow+j+nrow*( (iblock/nb_blocks_local) )] = V1block[j];


  free(V1block);

  }//end of last block
  else { break; }//error  //we can put the generic case here instead of between first and last cases
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
  }//end of if(nnew[id]>0)
  }//end of loop over the blocks


  return 0;
}


//====================================================================


// loop over the blocks to find which are useful for the local data

int get_overlapping_blocks_params(int nbloc, int *idv, int *n, int local_V_size, int nrow, int idp, int *idpnew, int *local_V_size_new, int *nnew, int *ifirstBlock, int *ilastBlock)
{

  int ib, nblockOK=0, firstrow, lastrow, nfullcol_data;
  int idptmp;

  nfullcol_data = (local_V_size-(nrow-idp%nrow)-(idp+local_V_size)%nrow)/nrow;  //check how many full columns input data have

//  int rank = PRINT_RANK;
//  FILE *file;
//  file = stdout;



  if( nfullcol_data > 0) {

  for( ib=0; ib<nbloc; ib++) {
    if( idv[ib] < nrow) {
      nnew[ib] = min( n[ib], nrow-idv[ib]);  //block used for the product
      nblockOK++;
    }
  }   

  }
  else {  //no full column observed 

    firstrow = idp%nrow;
    lastrow = (idp+local_V_size-1)%nrow;

    if( firstrow < lastrow) {  //just one column partially observed   

    for( ib=0; ib<nbloc; ib++) {
    if( (idv[ib]+n[ib] > firstrow) && (idv[ib] < lastrow+1)) {
      nnew[ib] = min( n[ib], nrow-idv[ib]);  //block used for the product
      nblockOK++;
    }
    }

    }
    else {  //two columns partially observed   

      for( ib=0; ib<nbloc; ib++) {
        if( (idv[ib]+n[ib] > firstrow) && (idv[ib] < nrow)) {  //intersects first partial column
          nnew[ib] = min( n[ib], nrow-idv[ib]);  //block used for the product
          nblockOK++;
        }

        if( (idv[ib] < lastrow+1) && (idv[ib]+n[ib] > 0)) {  //intersects second partial column
          nnew[ib] = min( n[ib], nrow-idv[ib]);  //block used for the product
          nblockOK++;  //may overcount but we do not care
        }  //could use else insteed!
      }
     }
  }

  if (VERBOSE)
    printf("nblockOK=%d\n", nblockOK);


  if( nblockOK == 0) return(0);  //no blocks overlapping with the data 

  //find the first and last relevant blocks for the begining and end of the local data  V

 //first block
  idptmp = idp;

  for( *ifirstBlock = -1; *ifirstBlock == -1;     ) {
    for(ib=0;ib<nbloc;ib++) {
      if(nnew[ib] != 0 && idptmp%nrow < idv[ib]+nnew[ib]) break;
    }

    if (ib<nbloc && idv[ib] <= idptmp%nrow) {
      *ifirstBlock = ib;
      *idpnew = idptmp;
    }
    else if (ib<nbloc && idv[ib] > idptmp%nrow) {
      *ifirstBlock = ib;
      int extrabegining = idv[ib]-idp%nrow;
//      *idpnew = idp+extrabegining;//idv[ib];
      int idvfirstcolumn = idptmp/nrow;
      *idpnew = idv[ib]+idvfirstcolumn*nrow;
    }
    else { //ib=nb_blocs
      idptmp += (int) (nrow-idptmp%nrow);
//          idtmp = (int) ceil((1.0*idpnew)/(1.0*nrow))*nrow; // go to the first element of the next column
    }}


 //last block
  idptmp = idp+local_V_size-1;

  for( *ilastBlock = -1; *ilastBlock == -1; ) {
    for(ib=nbloc-1;ib>=0;ib--) {
      if(nnew[ib] != 0 && idv[ib] <= idptmp%nrow) break;
    }


    if (ib>=0 && idptmp%nrow < idv[ib]+nnew[ib]) {
      *ilastBlock = ib;
      *local_V_size_new = local_V_size-(*idpnew)+idp;
    }
    else if (ib>=0 && idv[ib]+nnew[ib] <= idptmp%nrow) {
      *ilastBlock = ib;
      int extraend = (local_V_size-1+idp)%nrow+1-(idv[ib]+nnew[ib]);
      //*local_V_size_new = (local_V_size+idp)%nrow-(idv[*ilastBlock]+nnew[*ilastBlock]);
     //idv[*ilastBlock]+nnew[*ilastBlock]-(*idpnew);
      *local_V_size_new = local_V_size-(*idpnew)+idp-extraend;

      int idvlastcolumn = idptmp/nrow;
      *local_V_size_new = idv[ib]+nnew[ib]+idvlastcolumn*nrow - (*idpnew);

    }
    else {
      idptmp = (int) idptmp - (idptmp%nrow)-1;
//        idtmp = (int) floor( (1.0*idpnew)/(1.0*nrow))*nrow-1; // go to the last element of the previous column
    }}


    return(1);
}


//====================================================================


int mpi_gstbmm(double **V, int *n, int m, int nrow, double *T, int nb_blocks_local, int nb_blocks_all, int *lambda, int *idv, int id0p, int local_V_size, int *id0gap, int *lgap, int ngap, MPI_Comm comm)
{

  //MPI parameters
  int rank;  //process rank
  int size;  //process number

  MPI_Status status;
  MPI_Comm_rank(comm, &rank);
  MPI_Comm_size(comm, &size);

  int i,j,k;   //some indexes

  int flag_skip_build_gappy_blocks=0;

  FILE *file;
  file = stdout;

//put zeros at the gaps location
  reset_gaps( V, id0p, local_V_size, m, nrow, id0gap, lgap, ngap);

//allocation for the gappy structure of the diagonal block Toeplitz matrix
  int nb_blocks_gappy;
  double *Tgappy;
  int *idvgappy;
  int *ngappy;
  int *lambdagappy;

  int nb_blockgappy_max; 
  int Tgappysize_max;

//some computation usefull to determine the max size possible for the gappy variables
  int Tsize=0;
  int lambdamax=0;

  if (VERBOSE) 
    fprintf(file, "[%d] flag_skip_build_gappy_blocks=%d\n", rank, flag_skip_build_gappy_blocks);

  if (flag_skip_build_gappy_blocks==1) {  //no build gappy blocks strategy, just put zeros at gaps location

  //compute the product using only the input Toeplitz blocks structure with zeros at the gaps location
  mpi_stbmm(V, n, m, nrow, T, nb_blocks_local, nb_blocks_all, lambda, idv, id0p, local_V_size, MPI_COMM_WORLD);

  }
  else { //build gappy blocks strategy

  for(Tsize=i=0;i<nb_blocks_local;i++)
    Tsize += lambda[i];

  for(i=0;i<nb_blocks_local;i++) {
    if (lambda[i]>lambdamax)
      lambdamax = lambda[i];
  }

//compute max size possible for the gappy variables
  nb_blockgappy_max = nb_blocks_local+ngap;
  Tgappysize_max = Tsize + lambdamax*ngap;

//allocation of the gappy variables with max size possible
  Tgappy = (double *) calloc(Tgappysize_max,sizeof(double));
  ngappy = (int *) calloc(nb_blockgappy_max, sizeof(int));
  lambdagappy = (int *) calloc(nb_blockgappy_max, sizeof(int));
  idvgappy = (int *) calloc(nb_blockgappy_max, sizeof(int));


//build gappy Toeplitz block structure considering significant gaps locations, meaning we skip
//the gaps lower than the minimum correlation distance. You can also use the flag_param_distmin_fixed
//parameter which allows you to skip the gap lower than these value. Indeed, sometimes it's
//better to just put somes zeros than to consider two separates blocks.
//ps: This criteria could be dependant of the local lambda in futur impovements.
  int flag_param_distmin_fixed=0;
  build_gappy_blocks(n, m, nrow, T, nb_blocks_local, nb_blocks_all, lambda, idv, id0gap, lgap, ngap, &nb_blocks_gappy, Tgappy, idvgappy, ngappy, lambdagappy, flag_param_distmin_fixed);

  if (VERBOSE) {
    fprintf(file, "[%d] nb_blocks_gappy=%d\n", rank, nb_blocks_gappy);
      fprintf(file, "[%d] idvgappy[%d]=%d \n", rank, i, idvgappy[i]);
    for(i=0;i<nb_blocks_gappy;i++)
      fprintf(file, "[%d] ngappy[%d]=%d \n", rank, i, ngappy[i]);
  }
//ps: we could reallocate the gappy variables to their real size. Not sure it's worth it.

//compute the product using the freshly created gappy Toeplitz blocks structure
  mpi_stbmm(V, ngappy, m, nrow, Tgappy, nb_blocks_gappy, nb_blocks_gappy, lambdagappy, idvgappy, id0p, local_V_size, MPI_COMM_WORLD);

  } //end flag_skip_build_gappy_blocks==1


//put zeros on V for the gaps location again. Indeed, some gaps are just replaced by zeros
//in input, it's created some fakes results we need to clear after the computation.
  reset_gaps( V, id0p, local_V_size, m, nrow, id0gap, lgap, ngap);


  return 0;
}


//====================================================================

//put zeros on V for the gaps location
int reset_gaps(double **V, int id0,int local_V_size, int m, int nrow, int *id0gap, int *lgap, int ngap)
{
  int i,j,k;

  for (j=0 ; j<m; j++) {
  for (k=0 ; k<ngap; k++)
   for (i=0 ; i<lgap[k]; i++)

   if (id0gap[k]+i+j*nrow>=id0 && id0gap[k]+i+j*nrow-id0 <id0+local_V_size)
     (*V)[id0gap[k]+i+j*nrow-id0] = 0.;

  }


  return 0;
}

//====================================================================

int build_gappy_blocks(int *n, int m, int nrow, double *T, int nb_blocks_local, int nb_blocks_all, int *lambda, int *idv, int *id0gap, int *lgap, int ngap, int *nb_blocks_gappy_final, double *Tgappy, int *idvgappy, int *ngappy, int *lambdagappy, int flag_param_distmin_fixed)
{

  int i,j,k;
  int id,ib;
  int idtmp;
  int igapfirstblock, igaplastblock;

  int param_distmin=0;
  if (flag_param_distmin_fixed!=0)
    param_distmin = flag_param_distmin_fixed;

  int lambdaShft;

  int igaplastblock_prev=-1;
  int lambdaShftgappy=0;
  int offset_id = 0;

  int flag_igapfirstinside, flag_igaplastinside;
  int nbloc = nb_blocks_local;
  int nblocks_gappy=0;

  int idvtmp_firstblock;


  int nb_blockgappy_max = nb_blocks_local+ngap;
  int Tgappysize_max; 

  int ngappy_tmp;
  int lgap_tmp;

  int flag_gapok=0;

  int distcorr_min;

  int Tsize=0;
  int lambdamax=0;
  for(i=0;i<nb_blocks_local;i++) {
    Tsize += lambda[i];

    if (lambda[i]>lambdamax)  
      lambdamax = lambda[i];
  }

  Tgappysize_max = Tsize + lambdamax*ngap; 

  int Tgappysize=0;
  int k_prev=-1;

  for (k=0;k<ngap;k++) {

  //find block for the gap begining 
  for( igapfirstblock = -1; igapfirstblock == -1; ) {
    idtmp = id0gap[k];

    for(ib=0;ib<nbloc;ib++) {
      if(n[ib] != 0 && idtmp%nrow < idv[ib]+n[ib]) break;  }

    if (ib<nbloc && idv[ib] <= idtmp) {
      igapfirstblock = ib;  //the block contained the id0gap
      flag_igapfirstinside = 1;
    }
    else if (ib<nbloc && idv[ib] > idtmp) {
      igapfirstblock = ib;  //first block after the id0gap
      flag_igapfirstinside = 0;
    }
    else { //ib=nbloc 
      igapfirstblock = -2;  //no block defined
      flag_igapfirstinside = 0; 
    }}

  //find block for the end of the gap - reverse way
  for( igaplastblock = -1; igaplastblock == -1; ) {
    idtmp = id0gap[k]+lgap[k]-1;

    for(ib=nbloc-1;ib>=0;ib--) {
      if(n[ib] != 0 && idv[ib] <= idtmp) break;    }

    if (ib>=0 && idtmp < idv[ib]+n[ib]) {
      igaplastblock = ib;
      flag_igaplastinside = 1;
    }
    else if (ib>=0 && idv[ib]+n[ib] <= idtmp) {
      igaplastblock = ib;
      flag_igaplastinside = 0;
    }
    else { //ib=-1
      igaplastblock = -2;  //no block defined.
      flag_igaplastinside = 0; 
    }}


    if (igapfirstblock==igaplastblock)
      distcorr_min = lambda[igapfirstblock];
    else
      distcorr_min = 0;


//igapfirstblock != -2 && igaplastblock != -2 not really need but it's a shortcut
  if (lgap[k]> max(distcorr_min, param_distmin) && igapfirstblock != -2 && igaplastblock != -2) {

  idvtmp_firstblock = max( idv[igapfirstblock] , id0gap[k_prev]+lgap[k_prev]);

//test if the gap is ok for block reduce/split 
  if (igapfirstblock!=igaplastblock) {

    flag_gapok = 1;  //reduce the gap in each block. no pb if we add max() inside the ifs.
  }
  else if (id0gap[k]-idvtmp_firstblock>=lambda[igapfirstblock] && idv[igaplastblock]+n[igaplastblock]-(id0gap[k]+lgap[k])>=lambda[igaplastblock]) {

    flag_gapok = 1;
  }
  else if (igapfirstblock==igaplastblock){

  int ngappyleft_tmp = id0gap[k]-idvtmp_firstblock;
  int leftadd = max(0,lambda[igapfirstblock] - ngappyleft_tmp); 
  int ngappyright_tmp = idv[igaplastblock]+n[igaplastblock]-(id0gap[k]+lgap[k]);
  int rightadd = max(0,lambda[igapfirstblock] - ngappyright_tmp);
  int restgap = lgap[k] - (leftadd+rightadd);  

//  flag_gapok = (restgap>=0); 
  flag_gapok = (restgap >= max(0, param_distmin));

  }
  else {
  flag_gapok = 0;
  }


 //create gappy blocks if criteria is fullfill
  if (flag_gapok==1) {

  //copy the begining blocks
    for (id=igaplastblock_prev+1;id<igapfirstblock;id++) {

      for(lambdaShft=j=0; j<id; j++)  lambdaShft += lambda[j];

      for (i=0;i<lambda[id];i++)
        Tgappy[i+lambdaShftgappy] =  T[i+lambdaShft];

      lambdagappy[nblocks_gappy] = lambda[id];
      ngappy[nblocks_gappy] = n[id];
      idvgappy[nblocks_gappy] = idv[id];
      nblocks_gappy = nblocks_gappy + 1;
      lambdaShftgappy = lambdaShftgappy + lambda[id];

    }

  //clear last blockgappy if same block again - outside the "if" for border cases with n[]==0
    if (igaplastblock_prev==igapfirstblock && k!= 0) {
      lambdaShftgappy=lambdaShftgappy-lambda[igapfirstblock];
      nblocks_gappy = nblocks_gappy - 1;
//      idvtmp_firstblock = id0gap[k-1]+lgap[k-1]; //always exist because igaplastblock_prev!=-1
                                      //so not first turn - it's replace "idv[igapfirstblock]"
    }

  //reduce first block if defined
    if (flag_igapfirstinside==1 && (id0gap[k]-idvtmp_firstblock)>0) {  //check if inside and not on the border - meaning n[] not zero

      for(lambdaShft=j=0; j<igapfirstblock; j++)  lambdaShft += lambda[j];

      for (i=0;i<lambda[id];i++)   //copy the lambda's block
        Tgappy[i+lambdaShftgappy] =  T[i+lambdaShft];

      lambdagappy[nblocks_gappy] = lambda[id];
      ngappy[nblocks_gappy] = id0gap[k]-idvtmp_firstblock;
      ngappy[nblocks_gappy] = max( ngappy[nblocks_gappy], lambda[igapfirstblock]);

      idvgappy[nblocks_gappy] = idvtmp_firstblock;
      nblocks_gappy = nblocks_gappy + 1;
      lambdaShftgappy = lambdaShftgappy + lambda[igapfirstblock];

    }

  //reduce last block if defined
    if (flag_igaplastinside==1  && (idv[igaplastblock]+n[igaplastblock]-(id0gap[k]+lgap[k]))>0 ) {  //check if inside and not on the border - meaning n[] not zero

      for(lambdaShft=j=0; j<igaplastblock; j++)  lambdaShft += lambda[j];

      for (i=0;i<lambda[id];i++)   //duplicate the lambda's block
        Tgappy[i+lambdaShftgappy] =  T[i+lambdaShft];

      lambdagappy[nblocks_gappy] = lambda[id];
      ngappy[nblocks_gappy] = idv[igaplastblock]+n[igaplastblock]-(id0gap[k]+lgap[k]);
      int rightadd0 = max(0,lambda[igapfirstblock] - ngappy[nblocks_gappy]);

      ngappy[nblocks_gappy] = max( ngappy[nblocks_gappy] , lambda[igaplastblock]);

 //     idvgappy[nblocks_gappy] = id0gap[k]+lgap[k];
      idvgappy[nblocks_gappy] = id0gap[k]+lgap[k]-rightadd0;

      nblocks_gappy = nblocks_gappy + 1;
      lambdaShftgappy = lambdaShftgappy + lambda[igaplastblock];

    }

  igaplastblock_prev = igaplastblock;
  k_prev = k;
 
  }//end if (flag_gapok)
  }//end if (lgap[k]>param_distmin)
  }//end gap loop


 //now continu to copy the rest of the block left
  for (id=igaplastblock_prev+1;id<nb_blocks_local;id++) {

    for(lambdaShft=j=0; j<id; j++)  lambdaShft += lambda[j];

    for (i=0;i<lambda[id];i++)
      Tgappy[i+lambdaShftgappy] =  T[i+lambdaShft];

    lambdagappy[nblocks_gappy] = lambda[id];
    ngappy[nblocks_gappy] = n[id];
    idvgappy[nblocks_gappy] = idv[id];
    nblocks_gappy = nblocks_gappy + 1;
    lambdaShftgappy = lambdaShftgappy + lambda[id];

  }

  *nb_blocks_gappy_final = nblocks_gappy;  //just for output


  return 0;
}