toeplitz_gappy_seq.dev.c

Go to the documentation of this file.

 
 
//=========================================================================
// Alternave version of the sequential routine for the gaps - in dev
// Need to change the name to gap_reduce_matrix for more explicit purpose
// fd@apc
int gstmm(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, int *id0gap, int *lgap,
          int ngap) {
 
    // 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__);
 
    int     nnew    = 0;
    int     lcolnew = 0;
    int     lcol;
    double *Vcol;
    int     icol, igap;
    ;
    int id0col;
    int id0out = 0;
    int lnew;
 
    int nfullcol;
    nfullcol = (l - (n - id0 % n) - (id0 + l) % n)
             / n; // check how many full columns input data have
 
    int ifirstgap, ilastgap;
 
    if (id0 % n != 0) {
        // first column if not full
        id0out = 0;
        id0col = id0 % n;
        Vcol   = (*V) + id0col;
        lcol   = n - id0col;
 
        // find the first and last gap on the column
        for (k = 0; k < ngap; k++) {
            if (lgap[k] != 0 && id0col < id0gap[k] + lgap[k]) break;
        }
        ifirstgap = k;
 
        for (k = ngap - 1; k >= 0; k--) {
            if (lgap[k] != 0 && id0gap[k] <= id0col) break;
        }
        ilastgap = k;
 
 
        for (i = 0; i < ngap; i++) printf("id0gap[%d]=%d\n", i, id0gap[i]);
        for (i = 0; i < ngap; i++) printf("lgap[%d]=%d\n", i, lgap[i]);
 
        printf("---\n");
 
        printf("lcol=%d\n", lcol);
        for (i = 0; i < lcol; i++) printf("Vcol[%d]=%f\n", i, Vcol[i]);
 
        // reduce the gap to lambda zeros on the column
        gap_reduce(&Vcol, id0col, lcol, lambda, id0gap, lgap, ngap, &lcolnew,
                   id0out);
 
        id0out += lcolnew;
 
        printf("---\n");
 
        for (i = 0; i < lcol; i++) printf("Vcolout[%d]=%f\n", i, Vcol[i]);
        printf("===\n");
 
        printf("---\n");
 
        for (i = 0; i < l; i++) printf("(*V)[%d]=%f\n", i, (*V)[i]);
    }
 
    // generic loop on the full column
    for (icol = 0; icol < nfullcol; icol++) {
 
        id0out = lcolnew;
        id0col = 0;
        Vcol   = (*V) + id0col;
        lcol   = n - id0col;
 
        gap_reduce(&Vcol, id0col, lcol, lambda, id0gap, lgap, ngap, &lcolnew,
                   id0out);
 
        id0out += lcolnew;
 
        if (icol == 0) // take the nnew directly if they is a full column
            nnew = lcolnew;
    }
 
 
    // last column if not full and more than one column
    if (idf % n != n && m > 1) {
 
        id0out = 0;
        id0col = 0;
        Vcol   = (*V) + id0col;
        lcol   = idf % n;
 
        // find the first and last gap on the column
        ifirstgap = 0; // we already know it, no need to compute
 
        for (k = ngap - 1; k >= 0; k--) {
            if (lgap[k] != 0 && id0gap[k] <= id0col) break;
        }
        ilastgap = k;
 
        // reduce the gap to lambda zeros on the column
        gap_reduce(&Vcol, id0col, lcol, lambda, id0gap, lgap, ngap, &lcolnew,
                   id0out);
    }
 
 
    // cleaning the extra terms
    if (nnew == 0) { // compute the nnew if there is no full column
        for (igap = 0; igap < ngap; igap++)
            nnew += id0gap[igap] - min(id0gap[igap], lambda);
    }
 
    for (i = nnew * m; i < n * m; i++) (*V)[i] = 0;
 
    printf("===\n");
    for (i = 0; i < l; i++) printf("(*V)[%d]=%f\n", i, (*V)[i]);
 
 
    // now do the product
    //   int id0new = ..;
 
    //  stmm( V, nnew, m, id0new, lnew, T_fft, lambda, V_fft, V_rfft, plan_f,
    //  plan_b, blocksize, nfft)
 
 
    // gap extend - same thing as above and reverse gap_reduce routine
    //  gap_extend()
 
    // maybe create gap_blockreduce to clarify this routine.
 
    return 0;
}
 
//=========================================================================
 
 
int gap_reduce(double **V, int id0, int l, int lambda, int *id0gap, int *lgap,
               int ngap, int *newl, int id0out) {
    // routine variables
    int i, k, ii;
    int lg;
    int newid0gap_ip1;
    int newlgap_ip1;
 
 
    double *Vout;
    Vout = (*V) - (id0 - id0out);
 
 
    if (id0gap[0] > id0) {
        lg = id0gap[0] - id0;
        printf("lg=%d\n", lg);
        memmove(&(Vout)[0], &(*V)[0], lg * sizeof(double));
        // to do if the first element isn't a gap
    }
 
    int offset_first = id0gap[0] - id0 + min(lambda, lgap[0]);
    offset_first     = max(0, offset_first);
 
 
    //  for (k=0; k<offset_first; k++)
    //    (Vout)[k] = 0;
 
    printf("offset_first=%d\n", offset_first);
 
    for (i = 0; i < l; i++) printf("(Vout)[%d]=%f\n", i, (Vout)[i]);
    printf("---\n");
 
    printf("l=%d\n", l);
 
 
    for (i = 0; i < (ngap - 1); i++) {
        lg = id0gap[i + 1] - id0 - (id0gap[i] - id0 + lgap[i]);
        printf("lg=%d\n", lg);
        memmove(&(Vout)[offset_first], &(*V)[id0gap[i] - id0 + lgap[i]],
                lg * sizeof(double));
 
        for (ii = 0; ii < l; ii++) printf("(Vout)[%d]=%f\n", ii, (Vout)[ii]);
 
        newid0gap_ip1 = offset_first + lg;
        newlgap_ip1   = min(lambda, lgap[i + 1]);
        for (k = newid0gap_ip1; k < newid0gap_ip1 + newlgap_ip1; k++)
            (Vout)[k] = 0;
 
        printf("lg=%d ; lambda=%d ; lgap[i+1]=%d\n", lg, lambda, lgap[i + 1]);
        offset_first += lg + min(lambda, lgap[i + 1]);
    } // end gaps loop
 
    printf("---\n");
    printf("offset_first=%d\n", offset_first);
    for (i = 0; i < l; i++) printf("(Vout)[%d]=%f\n", i, (Vout)[i]);
 
    i = (ngap - 1);
 
    if (id0gap[i] - id0 + lgap[i] < l) {
        printf("toc\n");
        lg = l - (id0gap[i] - id0 + lgap[i]);
        memmove(&(Vout)[offset_first], &(*V)[id0gap[i] - id0 + lgap[i]],
                lg * sizeof(double));
        offset_first += lg;
        *newl = offset_first;
    } else {
        *newl = offset_first - min(lambda, lgap[ngap - 1])
              + min(lambda, (id0 + l - id0gap[ngap - 1]));
    }
 
    printf("---\n");
    printf("l=%d, *newl=%d, offset_first=%d\n", l, *newl, offset_first);
 
    for (i = offset_first; i < l; i++) (Vout)[i] = 0;
 
    printf("*newl=%d\n", *newl);
 
 
    return 0;
}
//=========================================================================
 
 
int gap_masking(double **V, int l, int *id0gap, int *lgap, int ngap) {
 
    int i, k;
    int lg;
 
    int offset_first = id0gap[0];
    for (i = 0; i < (ngap - 1); i++) {
        lg = id0gap[i + 1] - (id0gap[i] + lgap[i]);
        memmove(&(*V)[offset_first], &(*V)[id0gap[i] + lgap[i]],
                lg * sizeof(double));
        offset_first += lg;
    }
 
    i  = (ngap - 1);
    lg = l - (id0gap[i] + lgap[i]);
 
    memmove(&(*V)[offset_first], &(*V)[id0gap[i] + lgap[i]],
            lg * sizeof(double));
    offset_first += lg;
 
 
    for (i = offset_first; i < l; i++) (*V)[i] = 0;
 
 
    return 0;
}
 
 
//=========================================================================
 
 
int gap_filling(double **V, int l, int *id0gap, int *lgap, int ngap) {
 
    int i, k;
    int lg;
 
    int lgaptot = 0;
    for (i = 0; i < ngap; i++) lgaptot += lgap[i];
 
    int id0_Vmv;
    int offset_last;
 
    id0_Vmv     = id0gap[ngap - 1] + lgap[ngap - 1];
    offset_last = id0_Vmv - lgaptot;
    lg          = l - id0_Vmv;
 
    if (lg > 0)
        memmove(&(*V)[id0_Vmv], &(*V)[offset_last], lg * sizeof(double));
 
    for (i = ngap - 2; i >= 0; i--) {
        id0_Vmv = id0gap[i] + lgap[i];
        lg      = id0gap[i + 1] - id0_Vmv;
        offset_last -= lg;
 
        memmove(&(*V)[id0_Vmv], &(*V)[offset_last], lg * sizeof(double));
    }
 
    for (i = 0; i < ngap; i++)
        for (k = 0; k < lgap[i]; k++) (*V)[id0gap[i] + k] = 0;
 
    return 0;
}
 
 
// a naive version - in dev
int gap_masking_naive(double **V, int id0, int local_V_size, int m, int nrow,
                      int *id0gap, int *lgap, int ngap) {
    int i, j, k;
    int icol;
 
    int offsetV = id0;
    k           = 0;
 
    for (i = 0; i < local_V_size; i++) {
        icol = (i + id0) % nrow;
 
        if (icol < id0gap[k]) {
            (*V)[offsetV] = (*V)[i];
            offsetV       = offsetV + 1;
        } else if (icol >= id0gap[k] + lgap[k]) {
            k = k + 1;
            i = i - 1;
            //    (*V)[offsetV]=(*V)[i];
            //    offsetV=offsetV+1;
        } else { // do not copy, just skip
        }
 
 
        return 0;
    }