toeplitz_seq.c

/*
** file: toeplitz_seq.c, version 1.0, 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"
//#include "toeplitz_seq.h"
extern int NFFT;
extern int FFTW_FLAG;


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


int 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 parameters 
  int rank=-1;   //sequential use 

  FILE *file;
  file = stdout;

  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;


//---------------------------------------
//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];

//-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
//Generic case : Generic block of the process
  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 of if(nnew[id]>0)
  }//end of loop over the blocks


  return 0;
}


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


int 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)
{

  int rank=-1; //sequential use  

  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
  stbmm(V, n, m, nrow, T, nb_blocks_local, nb_blocks_all, lambda, idv, id0p, local_V_size);

  }
  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

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

  } //end flag_skip_build_gappy_blocks==1


//put zeros on V for the gaps location again. Indeed, some gaps are just replaced by zeros numbers 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;
}