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 = MPI_COMM_NULL;
 
    mpi_stbmm(V, nrow, m_cw, m_rw, tpltzblocks, nb_blocks, nb_blocks, idp,
              local_V_size, flag_stgy, comm);
 
 
    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 = MPI_COMM_NULL;
 
    mpi_gstbmm(V, nrow, m_cw, m_rw, tpltzblocks, nb_blocks, nb_blocks, idp,
               local_V_size, id0gap, lgap, ngap, flag_stgy, comm);
 
 
    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]=%ld ; 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;
}