I'm trying to use FFTW to compute fast summations, and I've run into
an issue:

int numFreq3 = numFreq*numFreq*numFreq;

FFTW_complex* dummy_sq_fft =
(FFTW_complex*)FFTW_malloc( sizeof(FFTW_complex)*numFreq3 );

FFTW_complex* dummy_sq =
(FFTW_complex*)FFTW_malloc( sizeof(FFTW_complex)*numFreq3 );

FFTW_complex* orig =
(FFTW_complex*)FFTW_malloc( sizeof(FFTW_complex)*numFreq3 );

FFTW_plan dummyPlan = FFTW_plan_dft_3d( numFreq, numFreq, numFreq,
orig, dummy_sq_fft, FFTW_FORWARD, FFTW_MEASURE );

FFTW_plan dummyInvPlan = FFTW_plan_dft_3d( numFreq, numFreq, numFreq,
dummy_sq_fft, dummy_sq, FFTW_BACKWARD, FFTW_MEASURE );

for(int i= 0; i < numFreq3; i++) {

orig[ i ][ 0 ] = sparseProfile02[ 0 ][ i ][ 0 ];

//img. part == 0

orig[ i ][ 1 ] = sparseProfile02[ 0 ][ i ] [ 1 ];

}

FFTW_execute(dummyPlan);

FFTW_execute(dummyInvPlan);

int count = 0;

for(int i=0; i<numFreq3; i++) {

double factor;
if(sparseProfile02[ 0 ] [ i ] [ 0 ])
   factor = dummy_sq[ i ][ 0 ]/sparseProfile02[ 0 ][ i ][ 0 ];

else factor = dummy_sq[ i ][ 0 ];

if(factor < 0) {

  count++;

  }

}

std::cout<<"Count "<<count<<"\n";

FFTW_free(dummy_sq_fft);
FFTW_free(dummy_sq);

FFTW_destroy_plan(dummyPlan);
FFTW_destroy_plan(dummyInvPlan);
(Here sparseProfile02[0] is of type FFTW_complex*, and contains only
positive real data.)

Since we have dummy_sq = IFFT(FFT(sparseProfile02[ 0 ])), we must have
dummy_sq = n^3*sparseProfile02. But this is true only some of the
time; in fact, values on the dummy_sq grid are negative whenever
corresponding values on the sparseProfile02 grid are zero (but not
vice-versa). Does anyone know why this is happening?

thanks,

rk

On 15/09/2010 21:36, brrk wrote:
> I'm trying to use FFTW to compute fast summations, and;;

This is grossly off topic: C++, and algorithmic.
CLC is for spliting hair about C minutiae. C++, and
almost anything involving using double to represent
an actual quantity (rather than a curios such as NaN)
is a no-no. I am NOT sarcastic.

That said, it looks at fist glance that you could have
a numerical statbility problem, similar to the well
studied one occuring in integer multiplication using
FFT with float or double.

For a start, think about what "equality" or "negative"
means for the result of a computation involving double.


   Francois Grieu

Francois Grieu <fgrieu@gmail.com> writes:
> On 15/09/2010 21:36, brrk wrote:
>> I'm trying to use FFTW to compute fast summations, and;;
>
> This is grossly off topic: C++, and algorithmic.
> CLC is for spliting hair about C minutiae. C++, and
> almost anything involving using double to represent
> an actual quantity (rather than a curios such as NaN)
> is a no-no. I am NOT sarcastic.

You may not be sarcastic, but you're largely mistaken.

The only C++ content in the post was the use of "std::cout << ...".
A bit of Googling indicates that FFTW is a C library.

And we've had extensive topical discussions here regarding
floating-point data and what it represents.

Having said that, there probably are better places to discuss this
particular topic.  Quesions about a third-party library will usually get
better reponses in a forum that dicusses the particular library than in
one that discusses the language it's written in.  (fftw.org doesn't
appear to host discussion forums, but they do encourage e-mail feedback;
you might consider asking them.)

(And the OP's code would be much easier to read with some indentation.)

> That said, it looks at fist glance that you could have
> a numerical statbility problem, similar to the well
> studied one occuring in integer multiplication using
> FFT with float or double.
>
> For a start, think about what "equality" or "negative"
> means for the result of a computation involving double.


-- 
Keith Thompson (The_Other_Keith) kst-u@mib.org  <http://www.ghoti.net/~kst>
Nokia
"We must do something.  This is something.  Therefore, we must do this."
    -- Antony Jay and Jonathan Lynn, "Yes Minister"

brrk <brrk.3001@gmail.com> wrote:
> I'm trying to use FFTW to compute fast summations, and I've run into
> an issue:

<Re-indented to improve readbility>

> int numFreq3 = numFreq*numFreq*numFreq;
> FFTW_complex* dummy_sq_fft =
>          (FFTW_complex*)FFTW_malloc( sizeof(FFTW_complex)*numFreq3 );
> FFTW_complex* dummy_sq =
>          (FFTW_complex*)FFTW_malloc( sizeof(FFTW_complex)*numFreq3 );
> FFTW_complex* orig =
>          (FFTW_complex*)FFTW_malloc( sizeof(FFTW_complex)*numFreq3 );

> FFTW_plan dummyPlan =
>     FFTW_plan_dft_3d( numFreq, numFreq, numFreq,
>                       orig, dummy_sq_fft, FFTW_FORWARD, FFTW_MEASURE );
> FFTW_plan dummyInvPlan =
>     FFTW_plan_dft_3d( numFreq, numFreq, numFreq,
>                       dummy_sq_fft, dummy_sq, FFTW_BACKWARD, FFTW_MEASURE );

> for (int i= 0; i < numFreq3; i++) {
>     orig[ i ][ 0 ] = sparseProfile02[ 0 ][ i ][ 0 ];
>     //img. part == 0
>     orig[ i ][ 1 ] = sparseProfile02[ 0 ][ i ] [ 1 ];
> }

> FFTW_execute(dummyPlan);
> FFTW_execute(dummyInvPlan);

> int count = 0;
> for (int i=0; i<numFreq3; i++) {
>     double factor;
>     if (sparseProfile02[ 0 ] [ i ] [ 0 ])
>         factor = dummy_sq[ i ][ 0 ]/sparseProfile02[ 0 ][ i ][ 0 ];
>     else
>         factor = dummy_sq[ i ][ 0 ];
>
>     if (factor < 0) {
>         count++;
>     }
> }

> std::cout<<"Count "<<count<<"\n";

> FFTW_free(dummy_sq_fft);
> FFTW_free(dummy_sq);

> FFTW_destroy_plan(dummyPlan);
> FFTW_destroy_plan(dummyInvPlan);

> (Here sparseProfile02[0] is of type FFTW_complex*, and contains only
> positive real data.)

> Since we have dummy_sq = IFFT(FFT(sparseProfile02[ 0 ])), we must have
> dummy_sq = n^3*sparseProfile02. But this is true only some of the
> time; in fact, values on the dummy_sq grid are negative whenever
> corresponding values on the sparseProfile02 grid are zero (but not
> vice-versa). Does anyone know why this is happening?

This is C++, I hope you realize (and this is comp.lang.c, not
comp.lang.c++). And it's about the result of some library, that
might be written in C but maybe not), not about the C language.

But what do you expect? That you get back the exact same values
you started with after some rather complex computations? That's
not how it works - computations with floating point values in-
troduce rounding errors and what you see here might be just a
result of that. How far off are the elements that were 0 at the
start compared to the deviations for elements that were non-
zero? If for example an element that was 1 before the forward
and backward transformation becomes maybe 0.9999999987 (after
adjustment for the n^3 factor) and an element that was 0 before
becomes -1.27e-8 than that's absolutely normal. Only if e.g.
elements that started off with 1 remain more or less 1 but
elements that were 0 to begin with become e.g. -0.7 then that
really would be something to worry about. In that case (and
after you made sure that you got everything else right) you
probably should consider talking to the guys that developed
that library, their email address is fftw@fftw.org. If you
do so better send them a complete, compilable program plus
the data you were using.
                              Regards, Jens
-- 
  \   Jens Thoms Toerring  ___      jt@toerring.de
   \__________________________      http://toerring.de

In article <a44973f0-bfba-4cb3-8ea7-
7e49fee5120c@k30g2000vbn.googlegroups.com>, brrk.3001@gmail.com says...
[snip]
You want news:sci.math.num-analysis

On Sep 15, 3:36=A0pm, j...@toerring.de (Jens Thoms Toerring) wrote:
> brrk <brrk.3...@gmail.com> wrote:
> > I'm trying to use FFTW to compute fast summations, and I've run into
> > an issue:
>
> <Re-indented to improve readbility>
>
>
>
>
>
> > int numFreq3 =3D numFreq*numFreq*numFreq;
> > FFTW_complex* dummy_sq_fft =3D
> > =A0 =A0 =A0 =A0 =A0(FFTW_complex*)FFTW_malloc( sizeof(FFTW_complex)*num=
Freq3 );
> > FFTW_complex* dummy_sq =3D
> > =A0 =A0 =A0 =A0 =A0(FFTW_complex*)FFTW_malloc( sizeof(FFTW_complex)*num=
Freq3 );
> > FFTW_complex* orig =3D
> > =A0 =A0 =A0 =A0 =A0(FFTW_complex*)FFTW_malloc( sizeof(FFTW_complex)*num=
Freq3 );
> > FFTW_plan dummyPlan =3D
> > =A0 =A0 FFTW_plan_dft_3d( numFreq, numFreq, numFreq,
> > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 orig, dummy_sq_fft, FFTW_FO=
RWARD, FFTW_MEASURE );
> > FFTW_plan dummyInvPlan =3D
> > =A0 =A0 FFTW_plan_dft_3d( numFreq, numFreq, numFreq,
> > =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 dummy_sq_fft, dummy_sq, FFT=
W_BACKWARD, FFTW_MEASURE );
> > for (int i=3D 0; i < numFreq3; i++) {
> > =A0 =A0 orig[ i ][ 0 ] =3D sparseProfile02[ 0 ][ i ][ 0 ];
> > =A0 =A0 //img. part =3D=3D 0
> > =A0 =A0 orig[ i ][ 1 ] =3D sparseProfile02[ 0 ][ i ] [ 1 ];
> > }
> > FFTW_execute(dummyPlan);
> > FFTW_execute(dummyInvPlan);
> > int count =3D 0;
> > for (int i=3D0; i<numFreq3; i++) {
> > =A0 =A0 double factor;
> > =A0 =A0 if (sparseProfile02[ 0 ] [ i ] [ 0 ])
> > =A0 =A0 =A0 =A0 factor =3D dummy_sq[ i ][ 0 ]/sparseProfile02[ 0 ][ i ]=
[ 0 ];
> > =A0 =A0 else
> > =A0 =A0 =A0 =A0 factor =3D dummy_sq[ i ][ 0 ];
>
> > =A0 =A0 if (factor < 0) {
> > =A0 =A0 =A0 =A0 count++;
> > =A0 =A0 }
> > }
> > std::cout<<"Count "<<count<<"\n";
> > FFTW_free(dummy_sq_fft);
> > FFTW_free(dummy_sq);
> > FFTW_destroy_plan(dummyPlan);
> > FFTW_destroy_plan(dummyInvPlan);
> > (Here sparseProfile02[0] is of type FFTW_complex*, and contains only
> > positive real data.)
> > Since we have dummy_sq =3D IFFT(FFT(sparseProfile02[ 0 ])), we must hav=
e
> > dummy_sq =3D n^3*sparseProfile02. But this is true only some of the
> > time; in fact, values on the dummy_sq grid are negative whenever
> > corresponding values on the sparseProfile02 grid are zero (but not
> > vice-versa). Does anyone know why this is happening?
>
> This is C++, I hope you realize (and this is comp.lang.c, not
> comp.lang.c++). And it's about the result of some library, that
> might be written in C but maybe not), not about the C language.
>
> But what do you expect? That you get back the exact same values
> you started with after some rather complex computations? That's
> not how it works - computations with floating point values in-
> troduce rounding errors and what you see here might be just a
> result of that. How far off are the elements that were 0 at the
> start compared to the deviations for elements that were non-
> zero? If for example an element that was 1 before the forward
> and backward transformation becomes maybe 0.9999999987 (after
> adjustment for the n^3 factor) and an element that was 0 before
> becomes -1.27e-8 than that's absolutely normal. Only if e.g.
> elements that started off with 1 remain more or less 1 but
> elements that were 0 to begin with become e.g. -0.7 then that
> really would be something to worry about. In that case (and
> after you made sure that you got everything else right) you
> probably should consider talking to the guys that developed
> that library, their email address is f...@fftw.org. If you
> do so better send them a complete, compilable program plus
> the data you were using.
> =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 =A0 Regards, Jens
> --
> =A0 \ =A0 Jens Thoms Toerring =A0___ =A0 =A0 =A0j...@toerring.de
> =A0 =A0\__________________________ =A0 =A0 =A0http://toerring.de


I agree that std::cout is actually c++, but as FFTW, the library in
question, is written in C, I thought I was better off posting to this
group. I apologize if I was off-topic.

In any case, the errors are too large to be ascribed to floating-point
computations alone.  For example, when the input value is zero, the
output value is sometimes -256 and sometimes -345.  There's definitely
a major problem here; I just wanted to know if it was an obvious one.

I'll take Jens's recommendation and email the guys at fftw.

thanks for your time,
rk

On Sep 15, 10:36=A0pm, brrk <brrk.3...@gmail.com> wrote:
> I'm trying to use FFTW to compute fast summations, and I've run into
> an issue:
>
I can tell ypu how to solve this type of problem generally.

1) Make sure your FFTW library is working in general. eg calculate a
forwards transformation and its inverse, in 1d. Are the results what
you expect. Move to 2d, then to 3d. Is the transform working as you
would expect on small datasets and does it scale up to large datasets.

2) Now try a trivial sum using your mathod. Do you get the right
answer?

If no, what's the simplest example that will produce the wrong answer?
If yes, do you still get the wrong answer on the test data? Can you
identify what is different about the test data that triggers the
error?

Usually this will lead you to the result, either uncover a bug in your
own code or, less likely, tell you unambigiously that the transform is
producing the wrong / insufficiently accurate result.