toeplitz_seq.c

Go to the documentation of this file.

#include "toeplitz.h"

//r1.1 - Frederic Dauvergne (APC)
//This is the sequential version of the mpi routines for the API.
//The stbmm and gstbmm are the same as the mpi_stbmm and mpi_gstbmm but without
//any communication. stmm is a simplifed call of the sequential product including
//initialization and cleaning. 


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

int stmm(double **V, int n, int m, double *T, int lambda, Flag flag_stgy)
{

//fftw variables
  fftw_complex *V_fft, *T_fft;
  double *V_rfft;
  fftw_plan plan_f, plan_b;

//product parameters
  int nfft, blocksize;

  FILE *file;
  file = stdout;


  tpltz_init(n, lambda, &nfft, &blocksize, &T_fft, T, &V_fft, &V_rfft, &plan_f, &plan_b, flag_stgy);

  //Toeplitz computation
  if(VERBOSE)
    fprintf(file, "Before stmm_main call : nfft = %d, blocksize = %d\n", nfft, blocksize);
  stmm_main(V, n, m, 0, n*m, T, T_fft, lambda, V_fft, V_rfft, plan_f, plan_b, blocksize, nfft, flag_stgy);

  tpltz_cleanup(&T_fft, &V_fft, &V_rfft, &plan_f, &plan_b);


  return 0;
}


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

int stbmm(double **V, int nrow, int m_cw, int m_rw, Block *tpltzblocks, int nb_blocks, int64_t idp, int local_V_size, Flag flag_stgy)
{


  int nb_blocks_local=nb_blocks;
  int nb_blocks_all=nb_blocks;
 // int idp=0;
 // int local_V_size=nrow;
 //MPI_Comm comm=NULL;

  mpi_stbmm(V, nrow, m_cw, m_rw, tpltzblocks, nb_blocks, nb_blocks, idp, local_V_size, flag_stgy, NULL);



  return 0;
}


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

// This matrix V contains defined gaps which represents the useless data for the comutation. The gaps indexes are defined in the global time space as the generized toeplitz matrix, 
// meaning the row dimension. Each of his diagonal blocks is a symmetric, band-diagonal Toeplitz matrix, which can be different for each block.
int gstbmm(double **V, int nrow, int m_cw, int m_rw, Block *tpltzblocks, int nb_blocks, int64_t idp, int local_V_size, int64_t *id0gap, int *lgap, int ngap, Flag flag_stgy)
{
  int nb_blocks_local=nb_blocks;
  int nb_blocks_all=nb_blocks;
 // int idp=0;
 // int local_V_size=nrow;
 //MPI_Comm comm=NULL;

  mpi_gstbmm(V, nrow, m_cw, m_rw, tpltzblocks, nb_blocks, nb_blocks, idp, local_V_size, id0gap, lgap, ngap, flag_stgy, NULL);



  return 0;
}


int gstbmm0(double **V, int nrow, int m, int m_rowwise, Block *tpltzblocks, int nb_blocks_local, int nb_blocks_all, int id0p, int local_V_size, int64_t *id0gap, int *lgap, int ngap, Flag flag_stgy)
{

  int rank=0;
  int i,j,k;   //some indexes

  int flag_skip_build_gappy_blocks = flag_stgy.flag_skip_build_gappy_blocks;

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

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


//allocation for the gappy structure of the diagonal block Toeplitz matrix
  int nb_blocks_gappy;
  int nb_blockgappy_max;
  int Tgappysize_max;

  Block *tpltzblocks_gappy;

//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
//to remake  stbmm(V, nrow, m, m_rowwise, tpltzblocks, nb_blocks_local, nb_blocks_all, id0p, local_V_size, flag_stgy);

  }
  else { //build gappy blocks strategy

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

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

//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
  tpltzblocks_gappy = (Block *) calloc(nb_blockgappy_max,sizeof(Block));


//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 = flag_stgy.flag_param_distmin_fixed;
  build_gappy_blocks(nrow, m, tpltzblocks, nb_blocks_local, nb_blocks_all, id0gap, lgap, ngap, tpltzblocks_gappy, &nb_blocks_gappy, flag_param_distmin_fixed);


if (VERBOSE) {
    fprintf(file, "[%d] nb_blocks_gappy=%d\n", rank, nb_blocks_gappy);
    for(i=0;i<nb_blocks_gappy;i++)
      fprintf(file, "[%d] idvgappy[%d]=%d ; ngappy[%d]=%d\n", rank, i, tpltzblocks_gappy[i].idv, i, tpltzblocks_gappy[i].n );
}
//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
//to remake  stbmm(V, nrow, m, m_rowwise, tpltzblocks_gappy, nb_blocks_local, nb_blocks_all, id0p, local_V_size, flag_stgy);

  } //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, m_rowwise, id0gap, lgap, ngap);


  return 0;
}