smdvquad.c

#include "smdvquad.h"

void smdv_quad(const mxArray *Q, const double **x, double* y) {
    const double *pQ, *px;
    mwIndex *irQ, *jcQ;
    mwSize col, total = 0;
    mwIndex star_row_idx, stop_row_idx, current_row_index;
    mwSize n;

    /* Get the starting positions of all four data arrays. */
    pQ = mxGetPr(Q);
    irQ = mxGetIr(Q);
    jcQ = mxGetJc(Q);
    n = mxGetN(Q);

    /* Retrieve the values of x */
    px = *x;

    /* Do le magic  */
    *y = 0;
    for (col = 0; col < n; col++) { // Iterate over all columns

        star_row_idx = jcQ[col];
        stop_row_idx = jcQ[col + 1];

        if (star_row_idx == stop_row_idx)
            continue;
        else {
            for (current_row_index = star_row_idx; current_row_index < stop_row_idx; current_row_index++) {
                (*y) += pQ[total++] * px[irQ[current_row_index]] * px[col];
            }
        }
    }
}