How to make it faster?

[istillshine]

When I used gprof to see which function consumed most running time, I
identified the following one. sz was less than 5000 on average, but
foo had been called about 1,000,000 times. I have tried using
"register sum = 0.0" and saw some improvement. My question is how to
improve foo further to make it faster.

double foo(double *a, double *b, int sz)
{

double sum = 0.0;
int i;

for (i = 0; i < sz; i++)
sum += a[i]*b[i];

return sum;
}

[Hallvard B Furuseth]

is*********@gma il.com writes:

>>Or compute exp(sum) instead of sum, if that doesn't overflow:
(...)

I came up with a solution.
exp_sum *= (z1/z2)*HUGE_VAL;
return sz*log(HUGE_VAL ) + log(exp_sum)

How large this HUGE_VAL should be?

Depends on your data set and the range of result values. Looks to me
like if z1/z2 are normally about the same size, HUGE_VAL should be 1...
I have little experience with numerical programming, but I suppose a
good value would be one which normally yields exp_sum == ca 1, i.e.
HUGE_VAL = exp(result/sz). But if you need this, maybe that means
the data set is too variable and you had better do log() in the loop.

--
Hallvard

[Hallvard B Furuseth]

is*********@gma il.com writes:

exp_sum *= (z1/z2)*HUGE_VAL;
return sz*log(HUGE_VAL ) + log(exp_sum)

Eh. I'm getting tired, but shouldn't that be exp_sum *= (z1/z2) /
HUGE_VAL or return log(exp_sum) - sz*log(HUGE_VAL )?

And a similar error carried over to my answer. (The idea was to
try to collect a normal value of the answer in HUGE_VAL so
so log(exp_sum) return approximately 0 for a normal data set.)

--
Hallvard

[istillshine]

On Apr 9, 1:19 pm, Hallvard B Furuseth <h.b.furus...@u sit.uio.no>
wrote:

istillsh...@gma il.com writes:

Can also get rid of the multiplications with exp0:
Instead of exp1, exp2, x0, ... z2, keep variables
with values (current exp1)/exp0 ... (current z2)/exp0.

I don't understand this

exp1 = exp(x->data[0] - 1.0);
exp2 = exp(x->data[1] - 1.0);
for (i = 0; i < sz; i++) {
x0 = g0[i];
y0 = f0[i];
(no other changes in the loop)
}

Variables <exp,x,y,z><0,1 ,2now have values (original value / exp0).
sum and df->data do not change (except due to rounding), since
they get values (foo/exp0) / (bar/exp0) instead of foo/bar.

Got it. So smart.

[istillshine]

On Apr 9, 2:29 pm, istillsh...@gma il.com wrote:

>
Combined with the previous 3 secs' saving, the running time has
reduced to 11 secs from 71 secs. Thank you. I really learned
something new and useful.

After setting -O2 flag in compiling, it only took 7 secs. The running
time is satisfactory. Your ideas are indispensable.

[istillshine]

On Apr 9, 2:29 pm, istillsh...@gma il.com wrote:

Combined with the previous 3 secs' saving, the running time has
reduced to 11 secs from 71 secs. Thank you. I really learned
something new and useful.

After setting -O2 flag in compiling, it only took 7 secs.

[istillshine]

On Apr 9, 2:50 pm, Paul Hsieh <websn...@gmail .comwrote:

#define ITER 5

Change this to:

#define ITER 64

Just trust me, and try it out.

I tried using "#define ITER 64". Overflow occurred in this setting.
I got "1.#J". I don't know what "1.#J" stands for? When I changed
the type of sump and sumd to long double, the problem disappeared.
The program saved 1 sec because of larger ITER.

I also tried another way:

#define ITER 64

register double sum = 0.0;

if (i % ITER == 0) {
sum += log(prodn/prodd); /* introducing log, but not so frequent */
prodn = 1.0;
prodd = 1.0;
}

In the way, the running time remained the same: 11 secs before setting
-O2, and 7 secs after setting -O2. However, the return value was a
little different.

[Paul Hsieh]

On Apr 9, 3:51 pm, istillsh...@gma il.com wrote:

On Apr 9, 2:50 pm, Paul Hsieh <websn...@gmail .comwrote:

#define ITER 5

Change this to:

#define ITER 64

Just trust me, and try it out.

I tried using "#define ITER 64". Overflow occurred in this setting.
I got "1.#J". I don't know what "1.#J" stands for?

Its a kind of NaN which is in fact an overflow condition of some
kind. In that case, you should try 4 or 8. But the key is to use a
power of 2. This will allow your compiler to use an "as-if" rule to
reduce the cost of executing the "%" operator. Otherwise you can use
a target counter:

if (i >= tgtI) {
tgtI += ITER;
prodn /= prodd;
prodd = 1.0;
}

That would probably work just as well.

[...] When I changed
the type of sump and sumd to long double, the problem disappeared.
The program saved 1 sec because of larger ITER.

I also tried another way:

#define ITER 64

register double sum = 0.0;

if (i % ITER == 0) {
sum += log(prodn/prodd); /* introducing log, but not so frequent */
prodn = 1.0;
prodd = 1.0;

}

In the way, the running time remained the same: 11 secs before setting
-O2, and 7 secs after setting -O2. However, the return value was a
little different.

Well I would suggest trying other values of ITER first.
Unfortunately, long double sounds like its use a larger likely more
expensive data type which may likely have impact on the performance of
rest of your code. So I would recommend against doing that.

--
Paul Hsieh
http://www.pobox.com/~qed/
http://bstring.sf.net/

[christian.bau]

On Apr 9, 7:50*pm, Paul Hsieh <websn...@gmail .comwrote:

Glad to hear it. *Transcendental calls are the worst in terms of
performance. *After that comes divides and modulos. *After that is
branch mis-predictions. *If you can optimize *all that, you will at
least be within 2x-4x of optimal (barring algorithm and IO
improvements) pretty much all time.

Usually the total cost is computational cost plus cost of data access.
Since the computational cost has by now been reduced by a factor ten,
it would now be a good idea to look at things like cache usage, but
that would require a bigger picture.

It was mentioned that g [i] is always one of 0.0, 0.25, 0.5 and 0.75.
A trivial change (changing these arrays from double to float) would
save 25% of memory bandwidth. A simple experiment would be to measure
the time of one program run, then do the same number of calls to foo,
but always with the same set of data. If that is a lot faster, then
cache improvements can help.

[user923005]

On Apr 11, 3:56*pm, "christian. bau" <christian....@ cbau.wanadoo.co .uk>
wrote:

On Apr 9, 7:50*pm, Paul Hsieh <websn...@gmail .comwrote:

Glad to hear it. *Transcendental calls are the worst in terms of
performance. *After that comes divides and modulos. *After that is
branch mis-predictions. *If you can optimize *all that, you will at
least be within 2x-4x of optimal (barring algorithm and IO
improvements) pretty much all time.

Usually the total cost is computational cost plus cost of data access.
Since the computational cost has by now been reduced by a factor ten,
it would now be a good idea to look at things like cache usage, but
that would require a bigger picture.

It was mentioned that g [i] is always one of 0.0, 0.25, 0.5 and 0.75.
A trivial change (changing these arrays from double to float) would
save 25% of memory bandwidth. A simple experiment would be to measure
the time of one program run, then do the same number of calls to foo,
but always with the same set of data. If that is a lot faster, then
cache improvements can help.

If the problem were scaled so that g[i] contained characters of
0,1,2,3 {which represent multiples of 0.25} it would reduce the
requirement for 4 byte float to 1 byte for char (assuming 32 bit float
and 8 bit char).