Skip to content

Commit

Permalink
Initialize vector index outside the loop, calculate it once inside th…
Browse files Browse the repository at this point in the history
…e loop, and access it as needed.
  • Loading branch information
aadler committed Aug 21, 2023
1 parent 3bbb127 commit 8e80096
Showing 1 changed file with 16 additions and 7 deletions.
23 changes: 16 additions & 7 deletions src/compHorner.c
Original file line number Diff line number Diff line change
Expand Up @@ -220,9 +220,9 @@ extern SEXP eftHorner_c(SEXP x, SEXP a) {
// The key here is to remember R stores matrices as contiguous vectors in
// COLUMN-MAJOR order. So to get to cell [i, j]---assuming C convention that
// everything starts at 0---you start by iterating over COLUMNS first and then
// actually want (numrows * rowIndex + colIndex). See the R code for how it
// works in R, but it's vectorized, so operations are done for all columns in
// one row at a time.
// actually want vectorindex = (numrows * rowIndex + colIndex). See the R code
// for how it works in R, but it's vectorized, so operations are done for all
// columns in one row at a time.
//
// Also, looking at how twoProd and twoSum work, only need to store A$x
// equivalent in variable as everything else is a result and not an input too.
Expand All @@ -232,16 +232,19 @@ extern SEXP eftHorner_c(SEXP x, SEXP a) {
// new results as the old results and then decrement the row counter.
//
// (AA: 2023-08-20)
int vecidx;
for (int j = nm1; j-- > 0; ) {
for (int i = 0; i < m; ++i) {
// Proper position of value given flattened column-major matrix.
vecidx = nm1 * i + j;
// First element of twoProd
Ax = ps0[i] * px[i];
// Second element of twoProd
ppiM[nm1 * i + j] = twoPrody(ps0[i], px[i]);
ppiM[vecidx] = twoPrody(ps0[i], px[i]);
// First element of twoSum
s1[i] = Ax + pa[j];
// Second element of twoSum
psigM[nm1 * i + j] = twoSumy(Ax, pa[j]);
psigM[vecidx] = twoSumy(Ax, pa[j]);
}
// Storing the new as old.
for (int i = 0; i < m; ++i) {
Expand Down Expand Up @@ -289,16 +292,22 @@ extern SEXP hornerSum_c(SEXP x, SEXP p, SEXP np, SEXP q) {

// See above for why this addressing schema is needed.
// Replace the 0's with the sum of the LAST rows of p and q (pi and Sigma).
int vecidx;
for (int i = 0; i < m; ++i) {
pr0[i] = pp[prows * (i + 1) - 1] + pq[prows * (i + 1) - 1];
// Proper position of value given flattened column-major matrix.
vecidx = prows * (i + 1) - 1;
pr0[i] = pp[vecidx] + pq[vecidx];
}

// See above for why this addressing schema is needed.
// Now build "upwards".
if (prows > 1) {
int vecidx;
for (int j = prows - 1; j-- > 0; ) {
for (int i = 0; i < m; ++i) {
r1[i] = pr0[i] * px[i] + (pp[prows * i + j] + pq[prows * i + j]);
// Proper position of value given flattened column-major matrix.
vecidx = prows * i + j;
r1[i] = pr0[i] * px[i] + (pp[vecidx] + pq[vecidx]);
}
for (int i = 0; i < m; ++i) {
pr0[i] = r1[i];
Expand Down

0 comments on commit 8e80096

Please sign in to comment.