dependencies

hi,

can anyone point me in the right direction to take reduce the
dependencies of this algorithm:

for (i=0; i < m; i++)
{

n1 = n2;
n2 = n2 + n2;
e = -6.283185307179586/n2;
a = 0.0;

for (j=0; j < n1; j++)
{
c = cos(a);
s = sin(a);
a = a + e;

for (k=j; k < n; k=k+n2)
{

{
t1 = c*x[k+n1] - s*y[k+n1];
t2 = s*x[k+n1] + c*y[k+n1];
x[k+n1] = x[k] - t1;
y[k+n1] = y[k] - t2;
x[k] = x[k] + t1;
y[k] = y[k] + t2;
}

}
}
}

On 2010-02-09, uche <uraniumore235@hotmail.com> wrote:
> hi,
>
> can anyone point me in the right direction to take reduce the
> dependencies of this algorithm:

I replaced all of your tabs with double-spaces. I apologise if you
liked them, but in my defense some newserver or archiver would have
eventually stripped them out anyway.

>
> for (i=0; i < m; i++)
> {
>
> n1 = n2;
> n2 = n2 + n2;
> e = -6.283185307179586/n2;

Well, n2 is completely unnecessary here. All it does is confuse
things. And you should have a PI constant. So here you would
have simply

e = -PI / n1;

and at the end of the loop,

n1 *= 2;

> a = 0.0;
>

Where do all these variables come from?

> for (j=0; j < n1; j++)
> {
> c = cos(a);
> s = sin(a);
> a = a + e;

How about
c = 1;
s = 0;
or actually, drop c and s because they're both meaningless constants.

>
> for (k=j; k < n; k=k+n2)

for(k = j; k < n; k += 2*n1)

Something seems very wrong about all of your loops, aside from
code style. Given that I have no clue what you're doing, I'll
leave it at that.

> {
>
> {
> t1 = c*x[k+n1] - s*y[k+n1];
> t2 = s*x[k+n1] + c*y[k+n1];

t1 = x[k + n1];
t2 = y[k + n1];

> x[k+n1] = x[k] - t1;
> y[k+n1] = y[k] - t2;
> x[k] = x[k] + t1;
> y[k] = y[k] + t2;

It seems to me that you need an awful lot of array bounds
checking to make this code safe.

> }
>
> }
> }
> }

HTH. HAND.

Andrew

uche wrote:

> hi,
>
> can anyone point me in the right direction to take reduce the
> dependencies of this algorithm:
>

Huh? The only dependencies I see are on math.h / lib (ala sin() & cos()).

Maybe you mean optimize?

--
*From the 1966 TV series:*
Robin: You can't get away from Batman that easy!
Batman: Easily.
Robin: Easily.
Batman: Good grammar is essential, Robin.

Andrew Poelstra <apoelstra@localhost.localdomain> writes:

> On 2010-02-09, uche <uraniumore235@hotmail.com> wrote:
>>
>> can anyone point me in the right direction to take reduce the
>> dependencies of this algorithm:
>
> I replaced all of your tabs with double-spaces. I apologise if you
> liked them, but in my defense some newserver or archiver would have
> eventually stripped them out anyway.
>
>>
>> for (i=0; i < m; i++)
>> {
>>
>> n1 = n2;
>> n2 = n2 + n2;
>> e = -6.283185307179586/n2;
>
> Well, n2 is completely unnecessary here. All it does is confuse
> things. And you should have a PI constant. So here you would
> have simply
>
> e = -PI / n1;
>
> and at the end of the loop,
>
> n1 *= 2;
>
>> a = 0.0;
>>
> Where do all these variables come from?
>
>> for (j=0; j < n1; j++)
>> {
>> c = cos(a);
>> s = sin(a);
>> a = a + e;
>
> How about
> c = 1;
> s = 0;
> or actually, drop c and s because they're both meaningless
> constants.

Eh? 'a' is incremented in this loop is it not? It ranges from 0 to
just shy of -PI unless I've missed something you've seen.

There was not enough context for me to unravel what this code was
doing. To the OP: some explanation may get you more helpful advice.

<snip>
--
Ben.

On 2010-02-10, Ben Bacarisse <ben.usenet@bsb.me.uk> wrote:
> Andrew Poelstra <apoelstra@localhost.localdomain> writes:
>> On 2010-02-09, uche <uraniumore235@hotmail.com> wrote:
>>>
>>> <snip>
>>>
>>> for (j=0; j < n1; j++)
>>> {
>>> c = cos(a);
>>> s = sin(a);
>>> a = a + e;
>>
>> How about
>> c = 1;
>> s = 0;
>> or actually, drop c and s because they're both meaningless
>> constants.
>
> Eh? 'a' is incremented in this loop is it not? It ranges from 0 to
> just shy of -PI unless I've missed something you've seen.
>

Oh, you're right. I noticed that the letter 'a' did not occur
after the line
a = a + e
and then deleted it for the wrong reason.

> There was not enough context for me to unravel what this code was
> doing. To the OP: some explanation may get you more helpful advice.
>
><snip>

Aye.

On Feb 9, 10:00*pm, Andrew Poelstra <apoels...@localhost.localdomain>
wrote:
> On 2010-02-10, Ben Bacarisse <ben.use...@bsb.me.uk> wrote:
>
>
>
> > Andrew Poelstra <apoels...@localhost.localdomain> writes:
> >> On 2010-02-09, uche <uraniumore...@hotmail.com> wrote:
>
> >>> <snip>
>
> >>> * * * for (j=0; j < n1; j++)
> >>> * * * {
> >>> * * * * c = cos(a);
> >>> * * * * s = sin(a);
> >>> * * * * a = a + e;
>
> >> How about
> >> * c = 1;
> >> * s = 0;
> >> or actually, drop c and s because they're both meaningless
> >> * constants.
>
> > Eh? *'a' is incremented in this loop is it not? *It ranges from 0 to
> > just shy of -PI unless I've missed something you've seen.
>
> Oh, you're right. I noticed that the letter 'a' did not occur
> after the line
> * a = a + e
> and then deleted it for the wrong reason.
>
> > There was not enough context for me to unravel what this code was
> > doing. *To the OP: some explanation may get you more helpful advice.
>
> ><snip>
>
> Aye.

Hi,

I got this code for this source: http://cnx.org/content/m12016/latest/

The code is towards the bottom ..

Any would like to strip it down so that I can do some meaningful
parallelism using openMP.

On Feb 10, 5:18*am, uche <uraniumore...@hotmail.com> wrote:
> On Feb 9, 10:00*pm, Andrew Poelstra <apoels...@localhost.localdomain>
> wrote:
>
>
>
> > On 2010-02-10, Ben Bacarisse <ben.use...@bsb.me.uk> wrote:
>
> > > Andrew Poelstra <apoels...@localhost.localdomain> writes:
> > >> On 2010-02-09, uche <uraniumore...@hotmail.com> wrote:
>
> > >>> <snip>
>
> > >>> * * * for (j=0; j < n1; j++)
> > >>> * * * {
> > >>> * * * * c = cos(a);
> > >>> * * * * s = sin(a);
> > >>> * * * * a = a + e;
>
> > >> How about
> > >> * c = 1;
> > >> * s = 0;
> > >> or actually, drop c and s because they're both meaningless
> > >> * constants.
>
> > > Eh? *'a' is incremented in this loop is it not? *It ranges from 0to
> > > just shy of -PI unless I've missed something you've seen.
>
> > Oh, you're right. I noticed that the letter 'a' did not occur
> > after the line
> > * a = a + e
> > and then deleted it for the wrong reason.
>
> > > There was not enough context for me to unravel what this code was
> > > doing. *To the OP: some explanation may get you more helpful advice..
>
> > ><snip>
>
> > Aye.
>
> Hi,
>
> I got this code for this source:http://cnx.org/content/m12016/latest/
>
> The code is towards the bottom ..
>
> Any would like to strip it down so that I can do some meaningful
> parallelism using openMP.

Although I am not too familiar with OpenMP's programming model, I
believe domain decomposition is what you're after. FFT can't be
decomposed functionally anyway. So for a distributed memory parallel
processing model, you'd distribute a multi-dimensional array across
the rows (for instance), transform all the dimensions of the data that
are local, then to perform the transpose by exchanging messages
between worker processes/threads. If the transpose aren't required to
be performed in-place, it's significantly simpler.

GL.

On Feb 10, 2:28*am, Michael Foukarakis <electricde...@gmail.com>
wrote:
> On Feb 10, 5:18*am, uche <uraniumore...@hotmail.com> wrote:
>
>
>
>
>
> > On Feb 9, 10:00*pm, Andrew Poelstra <apoels...@localhost.localdomain>
> > wrote:
>
> > > On 2010-02-10, Ben Bacarisse <ben.use...@bsb.me.uk> wrote:
>
> > > > Andrew Poelstra <apoels...@localhost.localdomain> writes:
> > > >> On 2010-02-09, uche <uraniumore...@hotmail.com> wrote:
>
> > > >>> <snip>
>
> > > >>> * * * for (j=0; j < n1; j++)
> > > >>> * * * {
> > > >>> * * * * c = cos(a);
> > > >>> * * * * s = sin(a);
> > > >>> * * * * a = a + e;
>
> > > >> How about
> > > >> * c = 1;
> > > >> * s = 0;
> > > >> or actually, drop c and s because they're both meaningless
> > > >> * constants.
>
> > > > Eh? *'a' is incremented in this loop is it not? *It ranges from0 to
> > > > just shy of -PI unless I've missed something you've seen.
>
> > > Oh, you're right. I noticed that the letter 'a' did not occur
> > > after the line
> > > * a = a + e
> > > and then deleted it for the wrong reason.
>
> > > > There was not enough context for me to unravel what this code was
> > > > doing. *To the OP: some explanation may get you more helpful advice.
>
> > > ><snip>
>
> > > Aye.
>
> > Hi,
>
> > I got this code for this source:http://cnx.org/content/m12016/latest/
>
> > The code is towards the bottom ..
>
> > Any would like to strip it down so that I can do some meaningful
> > parallelism using openMP.
>
> Although I am not too familiar with OpenMP's programming model, I
> believe domain decomposition is what you're after. FFT can't be
> decomposed functionally anyway. So for a distributed memory parallel
> processing model, you'd distribute a multi-dimensional array across
> the rows (for instance), transform all the dimensions of the data that
> are local, then to perform the transpose by exchanging messages
> between worker processes/threads. If the transpose aren't required to
> be performed in-place, it's significantly simpler.
>
> GL.- Hide quoted text -
>
> - Show quoted text -

Great help. Is there some paper or diagram that I can take a look at
in order to reprogram the code from above ?

On Feb 10, 3:51*pm, uche <uraniumore...@hotmail.com> wrote:
> On Feb 10, 2:28*am, Michael Foukarakis <electricde...@gmail.com>
> wrote:
>
>
>
> > On Feb 10, 5:18*am, uche <uraniumore...@hotmail.com> wrote:
>
> > > On Feb 9, 10:00*pm, Andrew Poelstra <apoels...@localhost.localdomain>
> > > wrote:
>
> > > > On 2010-02-10, Ben Bacarisse <ben.use...@bsb.me.uk> wrote:
>
> > > > > Andrew Poelstra <apoels...@localhost.localdomain> writes:
> > > > >> On 2010-02-09, uche <uraniumore...@hotmail.com> wrote:
>
> > > > >>> <snip>
>
> > > > >>> * * * for (j=0; j < n1; j++)
> > > > >>> * * * {
> > > > >>> * * * * c = cos(a);
> > > > >>> * * * * s = sin(a);
> > > > >>> * * * * a = a + e;
>
> > > > >> How about
> > > > >> * c = 1;
> > > > >> * s = 0;
> > > > >> or actually, drop c and s because they're both meaningless
> > > > >> * constants.
>
> > > > > Eh? *'a' is incremented in this loop is it not? *It ranges from 0 to
> > > > > just shy of -PI unless I've missed something you've seen.
>
> > > > Oh, you're right. I noticed that the letter 'a' did not occur
> > > > after the line
> > > > * a = a + e
> > > > and then deleted it for the wrong reason.
>
> > > > > There was not enough context for me to unravel what this code was
> > > > > doing. *To the OP: some explanation may get you more helpful advice.
>
> > > > ><snip>
>
> > > > Aye.
>
> > > Hi,
>
> > > I got this code for this source:http://cnx.org/content/m12016/latest/
>
> > > The code is towards the bottom ..
>
> > > Any would like to strip it down so that I can do some meaningful
> > > parallelism using openMP.
>
> > Although I am not too familiar with OpenMP's programming model, I
> > believe domain decomposition is what you're after. FFT can't be
> > decomposed functionally anyway. So for a distributed memory parallel
> > processing model, you'd distribute a multi-dimensional array across
> > the rows (for instance), transform all the dimensions of the data that
> > are local, then to perform the transpose by exchanging messages
> > between worker processes/threads. If the transpose aren't required to
> > be performed in-place, it's significantly simpler.
>
> > GL.- Hide quoted text -
>
> > - Show quoted text -
>
> Great help. Is there some paper or diagram that I can take a look at
> in order to reprogram the code from above ?

I'm not aware of any - maybe you could try google. I did an MPI
implementation a while ago but a quick glance through my archive
doesn't seem to find it . The FFTW library (http://www.fftw.org/)
has a parallel implementation, maybe you could check its source code
out to see if it fits your purposes.

In article <48064520-911a-4474-a650-60dce94b947d@z17g2000yqh.googlegroups.com>, uche <uraniumore235@hotmail.com> writes:

> I got this code for this source: http://cnx.org/content/m12016/latest/
>
> The code is towards the bottom ..
>
> Any would like to strip it down so that I can do some meaningful
> parallelism using openMP.

I think gcc might do some of that for you:

http://gcc.gnu.org/wiki/Graphite


But you may be able to transform the code manually as well (I'm saying
this without really looking at the code):

http://en.wikipedia.org/wiki/Polytope_model

Cheers,
lacos