Coverity: fftw module memory errors fixed
[scilab.git] / scilab / modules / fftw / sci_gateway / cpp / sci_fftw.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2012 - INRIA - Serge STEER
4 * Copyright (C) 2015 - Scilab Enterprises - Antoine ELIAS
5 *
6  * Copyright (C) 2012 - 2016 - Scilab Enterprises
7  *
8  * This file is hereby licensed under the terms of the GNU GPL v2.0,
9  * pursuant to article 5.3.4 of the CeCILL v.2.1.
10  * This file was originally licensed under the terms of the CeCILL v2.1,
11  * and continues to be available under such terms.
12  * For more information, see the COPYING file which you should have received
13  * along with this program.
14 *
15 */
16
17 #include "fftw_gw.hxx"
18 #include "function.hxx"
19 #include "fftw_common.hxx"
20
21 extern "C"
22 {
23 #include "callfftw.h"
24 #include "localization.h"
25 #include "charEncoding.h"
26 #include "Scierror.h"
27
28     int WITHMKL = 0;
29     extern void C2F(dscal)(int *n, double *da, double *dx, int *incx); /* blas routine */
30     extern void C2F(dset)(int *n, double *da, double *dx, int *incx); /* blas routine */
31 }
32 /*--------------------------------------------------------------------------*/
33 enum Scaling
34 {
35     Divide = -1,
36     None = 0,
37     Multiply = 1,
38 };
39 /*--------------------------------------------------------------------------*/
40 static int sci_fft_gen(const char *fname, types::Double* A, types::Double** O, int isn, guru_dim_struct gdim, int iopt);
41 /*--------------------------------------------------------------------------*/
42 types::Function::ReturnValue sci_fftw(types::typed_list &in, int _iRetCount, types::typed_list &out)
43 {
44     std::wstring name(L"fftw");
45     return fftw_common(name, in, _iRetCount, out, sci_fft_gen);
46 }
47 /*--------------------------------------------------------------------------*/
48 int sci_fft_gen(const char *fname, types::Double* A, types::Double** O, int isn, guru_dim_struct gdim, int iopt)
49 {
50     types::Double* tmp = A->clone()->getAs<types::Double>();
51     int ndimsA = tmp->getDims();
52     int *dimsA = tmp->getDimsArray();
53     double *Ar = tmp->get();
54     double *Ai = tmp->getImg();
55
56     /* Input  array variables */
57     int  isrealA = (Ai == NULL), issymA = 1, lA = 1;
58     /*for MKL*/
59     int isrealA_save = isrealA;
60
61     /*FFTW specific library variable */
62     enum Scaling scale = None;
63     enum Plan_Type type;
64     fftw_plan p;
65
66     /* input/output address for transform variables */
67     double *ri = NULL, *ii = NULL, *ro = NULL, *io = NULL;
68
69     /* for MKL special cases */
70     int * dims1 = NULL;
71     int * incr1 = NULL;
72
73     /* local variable */
74     int one = 1;
75     int i = 0;
76     int errflag = 0;
77
78     bool bFreeAi = false;
79
80     for (i = 0; i < ndimsA; i++)
81     {
82         lA *= dimsA[i];
83     }
84
85     if (iopt == 0)
86     {
87         /* automatically selected algorithm*/
88         issymA = check_array_symmetry(Ar, Ai, gdim);
89         if (issymA < 0)
90         {
91             Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
92             return 0;
93         }
94     }
95     else if (iopt == 1)
96     {
97         issymA = 1; /* user forces symmetry */
98     }
99     else
100     {
101         issymA = 0;
102     }
103
104     if (WITHMKL)
105     {
106         double dzero = 0.0;
107         if (isrealA)
108         {
109             /*MKL does not implement the r2c nor r2r guru split methods, make A complex */
110             if (issymA)
111             {
112                 /* result will be real, the imaginary part of A can be allocated alone */
113                 Ai = (double*)MALLOC(sizeof(double) * lA);
114                 C2F(dset)(&lA, &dzero, Ai, &one);
115                 bFreeAi = true;
116             }
117             else
118             {
119                 /* result will be complex, set clone of A complex for inplace computation */
120                 tmp->setComplex(true);
121                 Ai = tmp->getImg();
122                 C2F(dset)(&lA, &dzero, Ai, &one);
123                 isrealA = 0;
124             }
125         }
126     }
127
128     /* Set pointers on real and imaginary part of the input */
129     ri = Ar;
130     ii = Ai;
131
132     scale = None; /*no scaling needed */
133     if (isn == FFTW_BACKWARD)
134     {
135         scale = Divide;
136     }
137     if (isrealA & !WITHMKL) /* To have type = C2C_PLAN*/
138     {
139         /*A is real */
140         if (issymA)
141         {
142             ///*r2r =  isrealA &&  issymA*/
143             ///* there is no general plan able to compute r2r transform so it is tranformed into
144             //a R2c plan. The computed imaginary part will be zero*/
145             double dzero = 0.0;
146             *O = tmp;
147             tmp = NULL;
148             (*O)->setComplex(false);
149             ro = (*O)->get();
150             Ai = (double*)MALLOC(sizeof(double) * lA);
151             io = Ai;
152             bFreeAi = true;
153             C2F(dset)(&lA, &dzero, io, &one);
154             type = R2C_PLAN;
155         }
156         else
157         {
158             /*r2c =  isrealA && ~issymA;*/
159             /* transform cannot be done in place */
160             *O = new types::Double(ndimsA, dimsA, true);
161             ro = (*O)->get();
162             io = (*O)->getImg();
163             type = R2C_PLAN; /* fftw_plan_guru_split_dft_r2c plans for an FFTW_FORWARD transform*/
164             if (isn == FFTW_BACKWARD)
165             {
166                 /*transform problem into a FORWARD fft*/
167                 /*ifft(A)=conj(fft(A/N)) cas vect*/
168                 /* pre traitement A must be  divided by N cas vect*/
169                 /* post treatment result must conjugated */
170             }
171         }
172     }
173     else
174     {
175         *O = tmp;
176         tmp = NULL;
177         /* A is complex */
178         if (!WITHMKL && issymA) /*result is real*/
179         {
180             /*c2r =  ~isrealA &&  issymA*/
181             ro = ri;
182             io = NULL;
183
184             type = C2R_PLAN; /*fftw_plan_guru_split_dft_c2r plans for an FFTW_BACKWARD transform*/
185             if (isn == FFTW_FORWARD)
186             {
187                 /*transform problem into a BACKWARD fft : fft(A)=ifft(conj(A))*/
188                 double minusone = -1.0;
189                 C2F(dscal)(&lA, &minusone, ii, &one);
190             }
191         }
192         else
193         {
194             /*c2c =  ~isrealA && ~issymA;*/
195             /* use inplace transform*/
196             isrealA = 0;
197             type = C2C_PLAN; /*  fftw_plan_guru_split_dft plans for an FFTW_FORWARD transform*/
198             if (isn == FFTW_BACKWARD)
199             {
200                 /*transform problem into a FORWARD fft*/
201                 /* ifft(A) = %i*conj(fft(%i*conj(A)/N) */
202                 /* reverse input */
203                 ri = Ai;
204                 ii = Ar;
205                 /* reverse output */
206                 ro = Ai;
207                 io = Ar;
208             }
209             else
210             {
211                 ro = ri;
212                 io = ii;
213             }
214         }
215     }
216
217     /* pre-treatment */
218     if (scale != None)
219     {
220         double ak = 1.0;
221         for (i = 0; i < gdim.rank; i++)
222         {
223             ak = ak * ((double)(gdim.dims[i].n));
224         }
225         if (scale == Divide)
226         {
227             ak = 1.0 / ak;
228         }
229         C2F(dscal)(&lA, &ak, ri, &one);
230         if (isrealA == 0)
231         {
232             C2F(dscal)(&lA, &ak, ii, &one);
233         }
234     }
235
236     if (!WITHMKL || gdim.howmany_rank <= 1)
237     {
238         /* Set Plan */
239         p = GetFFTWPlan(type, &gdim, ri, ii, ro, io, getCurrentFftwFlags(), isn, (fftw_r2r_kind *)NULL, &errflag);
240         if (errflag)
241         {
242             if (errflag == 1)
243             {
244                 Scierror(999, _("%s: No more memory.\n"), fname);
245             }
246             else if (errflag == 2)
247             {
248                 Scierror(999, _("%s: Creation of requested fftw plan failed.\n"), fname);
249             }
250
251             delete (*O);
252             (*O) = NULL;
253
254             if (tmp)
255             {
256                 delete tmp;
257             }
258
259             if (bFreeAi)
260             {
261                 FREE(Ai);
262             }
263
264             return 0;
265         }
266
267         /* execute FFTW plan */
268         ExecuteFFTWPlan(type, p, ri, ii, ro, io);
269     }
270     else
271     {
272         /*FFTW MKL does not implement yet guru plan with howmany_rank>1             */
273         /*   associated loops described in gdim.howmany_rank and  gdim.howmany_dims */
274         /*   are implemented here by a set of call with howmany_rank==1             */
275         fftw_iodim *howmany_dims = gdim.howmany_dims;
276         int howmany_rank = gdim.howmany_rank;
277         int i1 = 0, i2 = 0;
278         int nloop = 0;
279         int t = 0;
280
281
282         gdim.howmany_rank = 0;
283         gdim.howmany_dims = NULL;
284
285         p = GetFFTWPlan(type, &gdim, ri, ii, ro, io, getCurrentFftwFlags(), isn, (fftw_r2r_kind *)NULL, &errflag);
286         if (errflag)
287         {
288             if (errflag == 1)
289             {
290                 Scierror(999, _("%s: No more memory.\n"), fname);
291             }
292             else if (errflag == 2)
293             {
294                 Scierror(999, _("%s: Creation of requested fftw plan failed.\n"), fname);
295             }
296
297             FREE(dims1);
298             FREE(incr1);
299
300             delete (*O);
301             (*O) = NULL;
302
303             if (tmp)
304             {
305                 delete tmp;
306             }
307
308             if (bFreeAi)
309             {
310                 FREE(Ai);
311             }
312
313             return 0;
314         }
315
316         /* flatten  nested loops: replace howmany_rank nested loops by a single one*/
317         /* Build temporary arrays used by flatened loop */
318         if ((dims1 = (int *)MALLOC(sizeof(int) * howmany_rank)) == NULL)
319         {
320             Scierror(999, _("%s: No more memory.\n"), fname);
321             FREE(dims1);
322             FREE(incr1);
323
324             delete (*O);
325             (*O) = NULL;
326
327             if (tmp)
328             {
329                 delete tmp;
330             }
331
332             if (bFreeAi)
333             {
334                 FREE(Ai);
335             }
336
337             return 0;
338         }
339         dims1[0] = howmany_dims[0].n;
340         for (i = 1; i < howmany_rank; i++)
341         {
342             dims1[i] = dims1[i - 1] * howmany_dims[i].n;
343         }
344         nloop = dims1[howmany_rank - 1];
345
346         if ((incr1 = (int *)MALLOC(sizeof(int) * howmany_rank)) == NULL)
347         {
348             Scierror(999, _("%s: No more memory.\n"), fname);
349             FREE(dims1);
350             FREE(incr1);
351
352             delete (*O);
353             (*O) = NULL;
354
355             if (tmp)
356             {
357                 delete tmp;
358             }
359
360             if (bFreeAi)
361             {
362                 FREE(Ai);
363             }
364
365             return 0;
366         }
367         t = 1;
368         for (i = 0; i < howmany_rank; i++)
369         {
370             t += (howmany_dims[i].n - 1) * howmany_dims[i].is;
371             incr1[i] = t;
372         }
373         /*loop on each "plan" */
374         i = 0; /*index on the first plan entry */
375         for (i1 = 1; i1 <= nloop; i1++)
376         {
377             /* the input and output are assumed to be complex because
378             within MKL real cases are transformed to complex ones in
379             previous steps of sci_fft_gen*/
380             ExecuteFFTWPlan(type, p, &ri[i], &ii[i], &ro[i], &io[i]);
381             i += howmany_dims[0].is;
382             /* check if  a loop ends*/
383             for (i2 = howmany_rank - 2; i2 >= 0; i2--)
384             {
385                 if ((i1 % dims1[i2]) == 0)
386                 {
387                     /*loop on dimension i2 ends, compute jump on the first plan entry index*/
388                     i += howmany_dims[i2 + 1].is - incr1[i2];
389                     break;
390                 }
391             }
392         }
393         /* free temporary arrays */
394         FREE(dims1);
395         FREE(incr1);
396         /* reset initial value of gdim for post treatment*/
397         gdim.howmany_rank = howmany_rank;
398         gdim.howmany_dims = howmany_dims;
399
400     }
401     /* Post treatment */
402     int iErr = 0;
403     switch (type)
404     {
405         case R2R_PLAN:
406             iErr = complete_array(ro, NULL, gdim);
407             break;
408         case C2R_PLAN:
409             break;
410         case R2C_PLAN:
411             if (issymA)
412             {
413                 /*R2C has been used to solve an r2r problem*/
414                 iErr = complete_array(ro, NULL, gdim);
415             }
416             else
417             {
418                 iErr = complete_array(ro, io, gdim);
419                 if (iErr != -1 && isn == FFTW_BACKWARD)
420                 {
421                     /*conjugate result */
422                     double ak = -1.0;
423                     C2F(dscal)(&lA, &ak, io, &one);
424                 }
425             }
426             break;
427         case C2C_PLAN:
428             if (WITHMKL && isrealA_save)
429             {
430                 if (isn == FFTW_FORWARD)
431                 {
432                     iErr = complete_array(ro, io, gdim);
433                 }
434                 else
435                 {
436                     iErr = complete_array(io, ro, gdim);
437                 }
438             }
439             break;
440     }
441
442     if (tmp)
443     {
444         delete tmp;
445     }
446
447     if (bFreeAi)
448     {
449         FREE(Ai);
450     }
451
452     if (iErr == -1)
453     {
454         Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
455         delete (*O);
456         (*O) = NULL;
457         return 0;
458     }
459
460     if (io == NULL)
461     {
462         (*O)->setComplex(false);
463     }
464
465     return 1;
466 }