toeplitz_block.c

Go to the documentation of this file.

#include "toeplitz.h"


//r1.1 - Frederic Dauvergne (APC)
//This is the routines related to the Toeplitz blocks diagonal routine.
//There is a sequential equivalent routine in the file toeplitz_seq.c

//todo:
//- remove the nooptimize communication

//=========================================================================
#ifdef W_MPI


int mpi_stbmm(double **V, int64_t nrow, int m, int m_rowwise, Block *tpltzblocks, int nb_blocks_local, int nb_blocks_all, int64_t idp, int local_V_size, Flag flag_stgy, MPI_Comm comm)
{
#else  //for sequential use only
int mpi_stbmm(double **V, int64_t nrow, int m, int m_rowwise, Block *tpltzblocks, int nb_blocks_local, int nb_blocks_all, int64_t idp, int local_V_size, Flag flag_stgy)
{
#endif


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

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

#else
  rank=0;
  size=1;
#endif

  PRINT_RANK=rank;

  FILE *file;
  file = stdout;

  int i,j,k;  //some indexes


  //identification of the mpi neighbours process to communicate when there is a shared block
  int right = rank+1;
  int left  = rank-1;


  //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));
  int64_t idpnew;
  int local_V_size_new;  
  int n_rowwise=local_V_size;

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


  if(PRINT_RANK==0 && VERBOSE>2)
    printf("status_params=%d\n", status_params);

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

  if (tpltzblocks[idv0].lambda==0 || tpltzblocks[idvn].lambda==0)
    return print_error_message (2, __FILE__, __LINE__);


  if(PRINT_RANK==0 && VERBOSE>2) {  //print on screen news parameters definition if VERBOSE
    fprintf(file, "new parameters caracteristics:\n");
    fprintf(file, "[%d] idp=%ld ; idpnew=%ld\n", rank, 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, (tpltzblocks[i].n), i, nnew[i]);
    for(i=0;i<nb_blocks_local;i++)
      fprintf(file, "[%d] tpltzblocks[%d].idv=%ld\n", rank, i, tpltzblocks[i].idv);
}


  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_blocks_rank;

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

  if(PRINT_RANK==0 && VERBOSE>2)
  fprintf(file, "[%d] nb_blocks_rank=%d, nb_blocks_local=%d\n", rank, nb_blocks_rank, nb_blocks_local);

  //Define the indices for the first and the last element in each blocks
  int idvp0 = idpnew%nrow-tpltzblocks[idv0].idv;  //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 = tpltzblocks[idvn].idv+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 offset0, offsetn;
  int distcorrmin_idv0 = (tpltzblocks[idv0].lambda)-1;
  int distcorrmin_idvn = (tpltzblocks[idvn].lambda)-1;

  //if(idvp0 != 0) 
    offset0 = min( idvp0, distcorrmin_idv0); 
  //if(idvpn != 0) 
    offsetn = min(idvpn, distcorrmin_idvn);  


  int toSendLeft=0;
  int toSendRight=0;

#ifdef W_MPI
  if(offset0!=0) {
    toSendLeft = min( tpltzblocks[idv0].idv+nnew[idv0]-idpnew%nrow, distcorrmin_idv0); 
  }
  if( offsetn != 0) {
    toSendRight = min( (idpnew+local_V_size_new)%nrow-tpltzblocks[idvn].idv, distcorrmin_idvn); 
  }

 int flag_optimlambda=1; //to allocate only the memory place needed

 int lambdaOut_offset;
 int lambdaIn_offset;
 double *LambdaOut;
 int lambdaOut_size, lambdaIn_size;

 if (flag_optimlambda==1) {
  LambdaOut=(double *) calloc((toSendLeft+toSendRight)*m_rowwise, sizeof(double));
  lambdaOut_offset = toSendLeft*m_rowwise;
  lambdaIn_offset = offset0*m_rowwise;
  lambdaOut_size = (toSendLeft+toSendRight)*m_rowwise ; 
  lambdaIn_size = (offset0+offsetn)*m_rowwise;
 }
 else {
  LambdaOut=(double *) calloc((tpltzblocks[idv0].lambda+tpltzblocks[idvn].lambda)*m_rowwise, sizeof(double));
  lambdaOut_offset = tpltzblocks[idv0].lambda*m_rowwise;
  lambdaIn_offset = tpltzblocks[idv0].lambda*m_rowwise;
  lambdaOut_size = (tpltzblocks[idv0].lambda+tpltzblocks[idvn].lambda)*m_rowwise;
  lambdaIn_size = (tpltzblocks[idv0].lambda+tpltzblocks[idvn].lambda)*m_rowwise;
 }


  if(offset0!=0) {
    for (j=0;j<m_rowwise;j++)
    for (i=0;i<toSendLeft;i++)
      LambdaOut[i+j*toSendLeft]=(*V)[i+j*n_rowwise]; //good because toSendLeft=0 if it 
  }                                                   //doesnt start on a the first block. 
  if( offsetn != 0) {
    for (j=0;j<m_rowwise;j++)
    for (i=0;i<toSendRight;i++)
      LambdaOut[i+j*toSendRight+lambdaOut_offset]=(*V)[i+j*n_rowwise+local_V_size-toSendRight];
  }                                        //good too using same argument than for offset0!=0
                                           //if local_V_size!=local_V_size_new+vShft mean there is extra
                                           //terms a the end and so offsetn=0
                                           //idpnew+local_V_size_new = idp+local_V_size and vShft=idpnew-idp
                                           //so local_V_size=vShft+local_V_size_new
  if(rank==0 || offset0==0)
    left = MPI_PROC_NULL;
  if(rank==size-1 || offsetn==0)
    right = MPI_PROC_NULL;

  double *LambdaIn=(double *) calloc(lambdaIn_size, sizeof(double));


 int flag_blockingcomm=0;  //to use blocking comm
 MPI_Request requestLeft_r, requestLeft_s;
 MPI_Request requestRight_r, requestRight_s;

  if (flag_blockingcomm==1) {
//send and receive data
  MPI_Sendrecv( LambdaOut, toSendLeft*m_rowwise, MPI_DOUBLE, left,  MPI_USER_TAG, (LambdaIn+lambdaIn_offset), offsetn*m_rowwise, MPI_DOUBLE, right, MPI_USER_TAG, comm, &status);
  MPI_Sendrecv( (LambdaOut+lambdaOut_offset), toSendRight*m_rowwise, MPI_DOUBLE, right,  MPI_USER_TAG, LambdaIn, offset0*m_rowwise, MPI_DOUBLE, left, MPI_USER_TAG, comm, &status);

  }
  else {
//to the Left
  MPI_Irecv((LambdaIn+lambdaIn_offset), offsetn*m_rowwise, MPI_DOUBLE, right, MPI_USER_TAG, comm, &requestLeft_r);
  MPI_Isend(LambdaOut, toSendLeft*m_rowwise, MPI_DOUBLE, left, MPI_USER_TAG, comm, &requestLeft_s);

//to the Right
  MPI_Irecv(LambdaIn, offset0*m_rowwise, MPI_DOUBLE, left, MPI_USER_TAG, comm, &requestRight_r);
  MPI_Isend((LambdaOut+lambdaOut_offset), toSendRight*m_rowwise, MPI_DOUBLE, right, MPI_USER_TAG, comm, &requestRight_s);

 }

#endif


//size of the first and the last block for the current process
  int v0rank_size, vnrank_size;
  if (nb_blocks_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 = tpltzblocks[idv0].idv + nnew[idv0] - idpnew%nrow + offset0;
    vnrank_size = ((idpnew+local_V_size_new-1)%nrow +1) - tpltzblocks[idvn].idv + offsetn;
  }


#ifdef W_MPI

if (flag_blockingcomm!=1) {
  //MPI_Wait for lambda comm
  MPI_Wait(&requestLeft_r, &status);
  MPI_Wait(&requestLeft_s, &status);
  MPI_Wait(&requestRight_r, &status);
  MPI_Wait(&requestRight_s, &status);

}


  free(LambdaOut);

#endif


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

  int idv1=0;     //old index of *V1
  int idv2=0;     //index 


  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;

  int jj;


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


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

#ifdef W_MPI
//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//first case : First block of the process
  if(iblock==idv0) {
  if(PRINT_RANK==0 && VERBOSE>2)
    fprintf(file, "[%d] First block...\n", rank);

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

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

  for (j=0;j<m_rowwise;j++) { 
#pragma omp parallel for //num_threads(NB_OMPTHREADS_STBMM)
  for (i=0;i<offset0;i++)
    V1block[i+j*vblock_size] = LambdaIn[i+j*offset0];
  }
//note: check if copyblock could be used instead.


//if (nb_blocks_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);

  for (j=0;j<m_rowwise;j++) {
#pragma omp parallel for //num_threads(NB_OMPTHREADS_STBMM)
  for (i=0;i<currentsize_middlepart;i++)
    V1block[offset0+i+j*vblock_size] = (*V)[i+vShft+j*n_rowwise];
  }

if (nb_blocks_rank == 1) {
  for (j=0;j<m_rowwise;j++) {
#pragma omp parallel for //num_threads(NB_OMPTHREADS_STBMM)
  for (i=0;i<offsetn;i++) {
    V1block[vblock_size-offsetn+i+j*vblock_size] = LambdaIn[i+lambdaIn_offset+j*offsetn];
  }}
}


  //init Toeplitz arrays
  tpltz_init(vblock_size, tpltzblocks[id].lambda, &nfft, &blocksize, &T_fft, (tpltzblocks[id].T_block), &V_fft, &V_rfft, &plan_f, &plan_b, flag_stgy);

  //Toeplitz computation
  if(PRINT_RANK==0 && VERBOSE>2)
    fprintf(file, "[%d] Before stmm_main call : nfft = %d, blocksize = %d\n", rank, nfft, blocksize);
  stmm_main(&V1block, vblock_size, m_rowwise, 0, m_rowwise*vblock_size, (tpltzblocks[id].T_block), T_fft, tpltzblocks[id].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);


  int currentsize=min(vblock_size-offset0, local_V_size_new);
  for (j=0;j<m_rowwise;j++) {
#pragma omp parallel for //num_threads(NB_OMPTHREADS_STBMM)
  for (i=0;i<currentsize;i++)
    (*V)[vShft+i+j*n_rowwise] = V1block[offset0+i+j*vblock_size];
  }

  free(V1block);

  }//end (First case)


//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//Generic case : Generic block of the process
  else if(iblock!=idv0 && iblock!=idv0+nb_blocks_rank-1) {
#endif

  if(PRINT_RANK==0 && VERBOSE>2)
    fprintf(file, "[%d] generic block...\n");

  vblock_size=nnew[id];
  id0block=0;

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

  idv1 = (tpltzblocks[id].idv)-idp%nrow - vShft + offset0 +nrow*( (iblock/nb_blocks_local) );  //no need
//  idv2 = idv[id]-idp%nrow + nrow*( (iblock/nb_blocks_local) );
  idv2 = (tpltzblocks[id].idv)-(idpnew)%nrow+vShft + nrow*( (iblock/nb_blocks_local) );

  for (j=0;j<m_rowwise;j++) {
#pragma omp parallel for //num_threads(NB_OMPTHREADS_STBMM)
  for (i=0;i<vblock_size;i++)
    V1block[i+j*vblock_size] = (*V)[i+idv2+j*n_rowwise];
//    V1block[i] = (*V)[i+idv1-offset0+vShft];
  }

  //init Toeplitz arrays
  tpltz_init(nnew[id], tpltzblocks[id].lambda, &nfft, &blocksize, &T_fft, (tpltzblocks[id].T_block), &V_fft, &V_rfft, &plan_f, &plan_b, flag_stgy);

  //Toeplitz computation
  if(PRINT_RANK==0 && VERBOSE>2)
    fprintf(file, "[%d] Before stmm_main call : nfft = %d, blocksize = %d\n", rank, nfft, blocksize);
  stmm_main(&V1block, vblock_size, m_rowwise, 0, m_rowwise*vblock_size, (tpltzblocks[id].T_block), T_fft, tpltzblocks[id].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);


  for (j=0;j<m_rowwise;j++) {
#pragma omp parallel for //num_threads(NB_OMPTHREADS_STBMM)
  for (i=0;i<vblock_size;i++) {
    (*V)[i+idv2+j*n_rowwise] = V1block[i+j*vblock_size];
  }}


  free(V1block);

#ifdef W_MPI
  }  //end (Generic case)

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

  vblock_size=vnrank_size;
  id0block=0;

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

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


  for (j=0;j<m_rowwise;j++) {
#pragma omp parallel for //num_threads(NB_OMPTHREADS_STBMM)
  for (i=0;i<vblock_size-offsetn;i++)
    V1block[i+j*vblock_size] = (*V)[i+idv2+j*n_rowwise];
//    V1block[i] = (*V)[i+idv1-offset0+vShft];
  }

  for (j=0;j<m_rowwise;j++) {
#pragma omp parallel for //num_threads(NB_OMPTHREADS_STBMM)
  for (i=0;i<offsetn;i++)
    V1block[vblock_size-offsetn+i+j*vblock_size] = LambdaIn[i+lambdaIn_offset+j*offsetn];
  }


  //init Toeplitz arrays
  tpltz_init(vblock_size, tpltzblocks[id].lambda, &nfft, &blocksize, &T_fft, (tpltzblocks[id].T_block), &V_fft, &V_rfft, &plan_f, &plan_b, flag_stgy);

  //Toeplitz computation
  if(PRINT_RANK==0 && VERBOSE>2)
    fprintf(file, "[%d] Before stmm_main call : nfft = %d, blocksize = %d\n", rank, nfft, blocksize);

  stmm_main(&V1block, vblock_size, m_rowwise, 0, vblock_size*m_rowwise, (tpltzblocks[id].T_block), T_fft, tpltzblocks[id].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);

  for (j=0;j<m_rowwise;j++) {
#pragma omp parallel for //num_threads(NB_OMPTHREADS_STBMM)                      
  for (i=0;i<vnrank_size-offsetn;i++) {
    (*V)[idv2+i+j*n_rowwise] = V1block[i+j*vblock_size];
  }}


  free(V1block);

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


  free(LambdaIn);


  return 0;
}

//#endif

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


int get_overlapping_blocks_params(int nbloc, Block *tpltzblocks, int local_V_size, int64_t nrow, int64_t idp, int64_t *idpnew, int *local_V_size_new, int *nnew, int *ifirstBlock, int *ilastBlock)
{

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


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


  if( nfullcol_data > 0) {

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

    }
    else {  //two columns partially observed   

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

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

  if(PRINT_RANK==0 && VERBOSE>2)
    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 < tpltzblocks[ib].idv+nnew[ib]) break;
    }

    if (ib<nbloc && tpltzblocks[ib].idv <= idptmp%nrow) {
      *ifirstBlock = ib;
      *idpnew = idptmp;
    }
    else if (ib<nbloc && tpltzblocks[ib].idv > idptmp%nrow) {
      *ifirstBlock = ib;
  //   int64_t extrabegining = tpltzblocks[ib].idv-idp%nrow;  //note I put int64 just to be sure. Never used
//      *idpnew = idp+extrabegining;//tpltzblocks[ib].idv;
      int idvfirstcolumn = idptmp/nrow;
      *idpnew = tpltzblocks[ib].idv+idvfirstcolumn*nrow;
    }
    else { //ib=nb_blocs
      idptmp += nrow-idptmp%nrow; //(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 && tpltzblocks[ib].idv <= idptmp%nrow) break;
    }


    if (ib>=0 && idptmp%nrow < tpltzblocks[ib].idv+nnew[ib]) {
      *ilastBlock = ib;
      *local_V_size_new = local_V_size-(*idpnew)+idp;
    }
    else if (ib>=0 && tpltzblocks[ib].idv+nnew[ib] <= idptmp%nrow) {
      *ilastBlock = ib;
      //int64_t extraend = (local_V_size-1+idp)%nrow+1-(tpltzblocks[ib].idv+nnew[ib]);  //note I put int64 just to be sure
      //*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;  //compute twice ... ? remove this one

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

    }
    else {
      idptmp = idptmp - (idptmp%nrow)-1;//(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);
}