at page 8 of the paper : [Reducing floating point error in dot product using the superblock family of algorithms- Anthony M. Castaldo , R.Clint Whaley, and Anthony T. Chronopouos][1]
there is the following code:
typedef double scalar;
typedef scalar *Vec;
scalar dotProd(Vec X, Vec Y, int nb)
{ int n = X.length;
int nblks = n / nb;
int nsblks = sqrt(nblks);
int blksInSblk = nblks / nsblks;
scalar dot = 0.0, sdot, cdot;
for (s=0; s < nsblks; s++)
{ sdot = 0.0;
for (b=0; b < blksInSblk; b++)
{
cdot = X[0] * Y[0];
for (i=1; i < nb; i++)
cdot += X[i] * Y[i];
sdot += cdot;
X += nb; Y += nb;
}
dot += sdot;
}
return dot;
}
I do not understand this line of the pseudocode:
X += nb; Y+= nb;
What does it mean if X
and Y
are two vectors and nb
is an integer?
Addendum: if not using pointers, how to "increment" by nb
elements?
int n = x.size();
int nblks = n / nb;
int nsblks = sqrt(nblks);
int blksInSblk = nblks / nsblks;
double dot = 0.0;
double sdot = 0.0;
double cdot = 0.0;
for(int s=0;s<nsblks;s++) {
std::cout << "iteration s= " << s << std::endl;
for(int b=0;b<blksInSblk;b++) {
std::cout << "iteration b= " << b << std::endl;
cdot += x(0) * y(0);
std::cout << "cdot += x(0) * y(0) = " << cdot << std::endl;
int pointingTo = 0;
for(int i=pointingTo;i<nb;i++) {
cdot += x(i) * y(i);
}
sdot += cdot;
// Increment the pointer to x by nb elements:
// x += nb;
pointingTo += nb;
// Increment the pointer to y by nb elements;
// y += nb
//}
}
dot += sdot;
}
The code you show can be interpreted as C code if you assume
typedef double scalar;
typedef scalar *Vec;
and add a line like int s, b, i;
somewhere before your first loop.
In that case, it is pretty clear that X += nb
is incrementing the pointer to X
by the size of nb
scalars. This is consistent with what the algorighm seems to be doing: summing up the dot product in chunks of size nb
instead of all at once.