Coverity: fftw module memory errors fixed
[scilab.git] / scilab / modules / fftw / src / cpp / fftw_common.cpp
1 /*
2 *  Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 *  Copyright (C) 2015 - Scilab Enterprises - Antoine ELIAS
4 *
5  * Copyright (C) 2012 - 2016 - Scilab Enterprises
6  *
7  * This file is hereby licensed under the terms of the GNU GPL v2.0,
8  * pursuant to article 5.3.4 of the CeCILL v.2.1.
9  * This file was originally licensed under the terms of the CeCILL v2.1,
10  * and continues to be available under such terms.
11  * For more information, see the COPYING file which you should have received
12  * along with this program.
13 *
14 */
15
16 #include "fftw_common.hxx"
17 #include "overload.hxx"
18 #include "int.hxx"
19 #include "string.hxx"
20
21 extern "C"
22 {
23 #include "callfftw.h"
24 #include "localization.h"
25 #include "charEncoding.h"
26 #include "Scierror.h"
27
28     extern int WITHMKL;
29 }
30 types::Function::ReturnValue fftw_common(std::wstring& name, types::typed_list &in, int _iRetCount, types::typed_list &out, fftw_gen func)
31 {
32     int iRhs = static_cast<int>(in.size());
33
34     wchar_t* option = nullptr;
35     int isn = FFTW_FORWARD;
36     WITHMKL = withMKL();
37     int iopt = 0;
38
39
40     char* s = wide_string_to_UTF8(name.data());
41     std::string cname(s);
42     FREE(s);
43
44     int iMaxRhs = 4;
45     if (name == L"fftw")
46     {
47         iMaxRhs = 5;
48     }
49     /****************************************
50     * Basic constraints on rhs arguments  *
51     ****************************************/
52
53     /* check min/max lhs/rhs arguments of scilab function */
54     if (in.size() < 1 || in.size() > iMaxRhs)
55     {
56         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), cname.data(), 1, iMaxRhs);
57         return types::Function::Error;
58     }
59
60     if (_iRetCount > 1)
61     {
62         Scierror(78, _("%s: Wrong number of output argument(s): %d expected.\n"), cname.data(), 1);
63         return types::Function::Error;
64     }
65
66     if (in[0]->isDouble() == false)
67     {
68         return Overload::generateNameAndCall(name, in, _iRetCount, out);
69     }
70
71     //check option, last parameter
72     if (in.back()->isString())
73     {
74         types::String* pOption = in.back()->getAs<types::String>();
75         if (pOption->isScalar())
76         {
77             option = pOption->get()[0];
78             --iRhs;
79         }
80         else
81         {
82             Scierror(999, _("%s: Cannot allocate more memory.\n"), cname.data());
83             return types::Function::Error;
84         }
85     }
86
87     /********************  Checking if isn is given  ************************************************/
88     if (iRhs != 1)
89     {
90         if (in[1]->isDouble() == false)
91         {
92             Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), cname.data(), 2);
93             return types::Function::Error;
94         }
95
96         types::Double* pWay = in[1]->getAs<types::Double>();
97         if (pWay->isScalar() == false)
98         {
99             Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), cname.data(), 2);
100             return types::Function::Error;
101         }
102
103         isn = static_cast<int>(pWay->get()[0]);
104         if (isn != FFTW_FORWARD && isn != FFTW_BACKWARD)
105         {
106             Scierror(53, _("%s: Wrong value for input argument #%d: %d or %d expected.\n"), cname.data(), 2, FFTW_FORWARD, FFTW_BACKWARD);
107             return types::Function::Error;
108         }
109     }
110
111     if (option)
112     {
113         if (cname == "dct" || cname == "dst")
114         {
115             if (isn == FFTW_FORWARD)
116             {
117                 if (option == name + L"1")
118                 {
119                     iopt = 1;
120                 }
121                 else if (option == name + L"2")
122                 {
123                     iopt = 2;
124                 }
125                 else if (option == name)
126                 {
127                     iopt = 0;
128                 }
129                 else if (option == name + L"4")
130                 {
131                     iopt = 4;
132                 }
133                 else
134                 {
135                     std::string err;
136                     err += "\"" + cname + "\",";
137                     err += "\"" + cname + "1\",";
138                     err += "\"" + cname + "2\",";
139                     err += "\"" + cname + "4\"";
140                     Scierror(999, _("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"), cname.data(), in.size(), err.data());
141                     return types::Function::Error;
142                 }
143             }
144             else
145             {
146                 if (option == name + L"1")
147                 {
148                     iopt = 1;
149                 }
150                 else if (option == name + L"3")
151                 {
152                     iopt = 3;
153                 }
154                 else if (option == L"i" + name)
155                 {
156                     iopt = 0;
157                 }
158                 else if (option == name + L"4")
159                 {
160                     iopt = 4;
161                 }
162                 else
163                 {
164                     std::string err;
165                     err += "\"i" + cname + "\",";
166                     err += "\"" + cname + "1\",";
167                     err += "\"" + cname + "3\",";
168                     err += "\"" + cname + "4\"";
169                     Scierror(999, _("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"), cname.data(), in.size(), err.data());
170                     return types::Function::Error;
171                 }
172             }
173         }
174         else //fftw
175         {
176             if (option == std::wstring(L"symmetric"))
177             {
178                 iopt = 1;
179             }
180             else if (option == std::wstring(L"nonsymmetric"))
181             {
182                 iopt = 2;
183             }
184             else
185             {
186                 Scierror(999, _("%s: Wrong value for input argument #%d: '%s' or '%s' expected.\n"), name.data(), iRhs, "\"symmetric\"", "\"nonsymmetric\"");
187                 return types::Function::Error;
188             }
189         }
190     }
191
192     switch (iRhs)
193     {
194         case 4:
195         {
196             /* dct(A ,sign ,dim, incr)*/
197             return common_4args(cname, in, _iRetCount, out, func, isn, iopt);
198             break;
199         }
200         case 3:
201         {
202             /* dct(A ,sign ,sel)*/
203             return common_3args(cname, in, _iRetCount, out, func, isn, iopt);
204             break;
205         }
206         default:
207         {
208             /* dct(A ,sign)*/
209             return common_2args(cname, in, _iRetCount, out, func, isn, iopt);
210             break;
211         }
212     }
213 }
214
215 types::Function::ReturnValue common_2args(std::string& name, types::typed_list &in, int _iRetCount, types::typed_list &out, fftw_gen func, int way, int opt)
216 {
217     //getting the array A
218     types::Double* A = in[0]->getAs<types::Double>();
219     int ndimsA = A->getDims();
220     int *dimsA = A->getDimsArray();
221
222     /*FFTW specific library variable */
223     guru_dim_struct gdim = {0, NULL, 0, NULL};
224
225     /* local variable */
226     int ndims = 0; /* number of non singleton dimensions */
227     int first_nonsingleton = -1;
228     int i = 0, j = 0;
229     int prd = 1;
230
231     /* ignore singleton dimensions */
232     first_nonsingleton = -1;
233     ndims = 0;
234     for (i = 0; i < ndimsA; i++)
235     {
236         if (dimsA[i] > 1)
237         {
238             ndims++;
239             if (first_nonsingleton < 0)
240             {
241                 first_nonsingleton = i;
242             }
243         }
244     }
245
246     /* void or scalar input gives void output or scalar*/
247     if (ndims == 0)
248     {
249         out.push_back(A);
250         return types::Function::OK;
251     }
252
253     gdim.rank = ndims;
254     if ((gdim.dims = (fftw_iodim *)MALLOC(sizeof(fftw_iodim) * gdim.rank)) == nullptr)
255     {
256         Scierror(999, _("%s: Cannot allocate more memory.\n"), name.data());
257         FREE(gdim.dims);
258         FREE(gdim.howmany_dims);
259         return types::Function::Error;
260     }
261
262     j = 0;
263     prd = 1;
264     for (i = (first_nonsingleton); i < ndimsA; i++)
265     {
266         if (dimsA[i] > 1)
267         {
268             gdim.dims[j].n = dimsA[i];
269             gdim.dims[j].is = prd;
270             gdim.dims[j].os = prd;
271             prd *= dimsA[i];
272             j++;
273         }
274     }
275     gdim.howmany_rank = 0;
276     gdim.howmany_dims = nullptr;
277
278
279     //fftw functions work "in place", so we need to clone input data.
280     types::Double* D = nullptr;// A->clone()->getAs<types::Double>();
281     if (!func(name.data(), A, &D, way, gdim, opt))
282     {
283         FREE(gdim.dims);
284         FREE(gdim.howmany_dims);
285         return types::Function::Error;
286     }
287     FREE(gdim.dims);
288     FREE(gdim.howmany_dims);
289     out.push_back(D);
290     return types::Function::OK;
291 }
292
293 types::Function::ReturnValue common_3args(std::string& name, types::typed_list &in, int _iRetCount, types::typed_list &out, fftw_gen func, int way, int opt)
294 {
295     //getting the array A
296     types::Double* A = in[0]->getAs<types::Double>();
297     int ndimsA = A->getDims();
298     int *dimsA = A->getDimsArray();
299
300     int *Sel = nullptr;
301     int rank = 0;
302
303     /*FFTW specific library variable */
304     guru_dim_struct gdim = {0, NULL, 0, NULL};
305     /* local variable */
306     int ndims = 0;
307     int first_nonsingleton = -1;
308     int ih = 0;
309     int pd = 1; /* used to store prod(Dims(1:sel(k-1)))*/
310     int pds = 1; /* used to store prod(Dims(sel(k-1):sel(k)))*/
311     int i = 0, j = 0;
312
313     /* ignore singleton dimensions */
314     first_nonsingleton = -1;
315     ndims = 0;
316     for (i = 0; i < ndimsA; i++)
317     {
318         if (dimsA[i] > 1)
319         {
320             ndims++;
321             if (first_nonsingleton < 0)
322             {
323                 first_nonsingleton = i;
324             }
325         }
326     }
327
328     /* void or scalar input gives void output or scalar*/
329     if (ndims == 0)
330     {
331         out.push_back(A);
332         return types::Function::OK;
333     }
334
335     if (in[2]->isGenericType() == false)
336     {
337         Scierror(999, _("%s: Wrong type for input argument #%d.\n"), name.data(), 3);
338         FREE(gdim.dims);
339         FREE(gdim.howmany_dims);
340         return types::Function::Error;
341     }
342
343     getVarAsDims(in[2], rank, Sel);
344
345     /* size of Sel must be less than ndimsA */
346     if (rank <= 0 || rank >= ndimsA)
347     {
348         Scierror(999, _("%s: Wrong size for input argument #%d: Must be between %d and %d.\n"), name.data(), 3, 1, ndimsA - 1);
349         FREE(gdim.dims);
350         FREE(gdim.howmany_dims);
351         FREE(Sel);
352         return types::Function::Error;
353     }
354     /* check values of Sel[i] */
355     for (i = 0; i < rank; i++)
356     {
357         if (Sel[i] <= 0)
358         {
359             Scierror(999, _("%s: Wrong values for input argument #%d: Positive integers expected.\n"), name.data(), 3);
360             FREE(gdim.dims);
361             FREE(gdim.howmany_dims);
362             FREE(Sel);
363             return types::Function::Error;
364         }
365         if (Sel[i] > ndimsA)
366         {
367             Scierror(999, _("%s: Wrong values for input argument #%d: Elements must be less than %d.\n"), name.data(), 3, ndimsA);
368             FREE(gdim.dims);
369             FREE(gdim.howmany_dims);
370             FREE(Sel);
371             return types::Function::Error;
372         }
373         if (i > 0 && Sel[i] <= Sel[i - 1])
374         {
375             Scierror(999, _("%s: Wrong values for input argument #%d: Elements must be in increasing order.\n"), name.data(), 3);
376             FREE(gdim.dims);
377             FREE(gdim.howmany_dims);
378             FREE(Sel);
379             return types::Function::Error;
380         }
381     }
382
383     /* Create  gdim struct */
384     gdim.rank = rank;
385     if ((gdim.dims = (fftw_iodim *)MALLOC(sizeof(fftw_iodim) * gdim.rank)) == nullptr)
386     {
387         Scierror(999, _("%s: Cannot allocate more memory.\n"), name.data());
388         FREE(gdim.dims);
389         FREE(gdim.howmany_dims);
390         FREE(Sel);
391         return types::Function::Error;
392     }
393
394     pd = 1; /* used to store prod(Dims(1:sel(k-1)))*/
395     pds = 1; /* used to store prod(Dims(sel(k-1):sel(k)))*/
396     j = 0;
397     for (i = 0; i < ndimsA; i++)
398     {
399         if (j >= gdim.rank)
400         {
401             break;
402         }
403         if (Sel[j] == i + 1)
404         {
405             gdim.dims[j].n = dimsA[i];
406             gdim.dims[j].is = pd;
407             gdim.dims[j].os = pd;
408             j++;
409         }
410         pd *= dimsA[i];
411     }
412     /* Compute howmany_rank based on jumps in the Sel sequence */
413     gdim.howmany_rank = 0;
414     if ((Sel[0] != 1) && (Sel[0] != ndimsA))
415     {
416         gdim.howmany_rank++;
417     }
418     for (i = 1; i <= rank - 1; i++)
419     {
420         if (Sel[i] != Sel[i - 1] + 1)
421         {
422             /*check if all dimensions between Sel[i-1]+1 and Sel[i]-1 are
423             equal to one, in this case they can be ignored and there is
424             no jump*/
425             for (j = Sel[i - 1] + 1; j <= Sel[i] - 1; j++)
426             {
427                 if (dimsA[j - 1] != 1)
428                 {
429                     gdim.howmany_rank++;
430                     break;
431                 }
432             }
433         }
434     }
435
436     if ((Sel[rank - 1] != ndimsA) || (rank == 1))
437     {
438         gdim.howmany_rank++;
439     }
440     /* Fill the howmany_dims struct */
441     if (gdim.howmany_rank > 0)
442     {
443         /* it must be the case */
444         if ((gdim.howmany_dims = (fftw_iodim *)MALLOC(gdim.howmany_rank * sizeof(fftw_iodim))) == nullptr)
445         {
446             Scierror(999, _("%s: Cannot allocate more memory.\n"), name.data());
447             FREE(gdim.dims);
448             FREE(gdim.howmany_dims);
449             FREE(Sel);
450             return types::Function::Error;
451         }
452         pd = 1;
453         for (j = 1; j <= (Sel[0] - 1); j++)
454         {
455             pd *= dimsA[j - 1];    /*prod(Dims(1:(sel(1)-1)))*/
456         }
457         ih = 0;
458         if ((Sel[0] != 1) && (Sel[0] != ndimsA))
459         {
460             /* First seleted dimension */
461             gdim.howmany_dims[ih].is = 1;
462             gdim.howmany_dims[ih].os = 1;
463             gdim.howmany_dims[ih].n = pd;
464             ih++;
465         }
466         pd *= dimsA[Sel[0] - 1]; /*prod(Dims(1:sel(1)))*/
467         for (i = 1; i <= rank - 1; i++)
468         {
469             /* intermediate selected dimensions */
470             if (Sel[i] != Sel[i - 1] + 1)
471             {
472                 pds = 1;
473                 for (j = (Sel[i - 1] + 1); j <= (Sel[i] - 1); j++)
474                 {
475                     pds *= dimsA[j - 1];    /*prod(Dims(sel(i-1)+1:(sel(i)-1)))*/
476                 }
477                 /*check again if all dimensions between Sel[i-1]+1 and
478                 Sel[i]-1 are equal to one, in this case they can be
479                 ignored and there is no jump*/
480                 for (j = (Sel[i - 1] + 1); j <= (Sel[i] - 1); j++)
481                 {
482                     if (dimsA[j - 1] != 1)
483                     {
484                         gdim.howmany_dims[ih].is = pd;
485                         gdim.howmany_dims[ih].os = pd;
486                         gdim.howmany_dims[ih].n = pds;
487                         ih++;
488                         break;
489                     }
490                 }
491             }
492             pd *= pds * dimsA[Sel[i] - 1]; /*prod(Dims(1:sel(i)))*/
493         }
494
495         if (Sel[rank - 1] != ndimsA)
496         {
497             /* last selected dimension*/
498             pds = 1;
499             for (j = (Sel[rank - 1] + 1); j <= ndimsA; j++)
500             {
501                 pds *= dimsA[j - 1];    /*prod(Dims(sel(i-1)+1:(sel(i)-1)))*/
502             }
503             gdim.howmany_dims[ih].is = pd;
504             gdim.howmany_dims[ih].os = pd;
505             gdim.howmany_dims[ih].n = pds;
506             ih++;
507         }
508         else if (rank == 1)
509         {
510             /* the only selected dimension is the last one */
511             gdim.howmany_dims[ih].is = 1;
512             gdim.howmany_dims[ih].os = 1;
513             gdim.howmany_dims[ih].n = pd / dimsA[Sel[0] - 1];
514             ih++;
515         }
516     }
517
518     //fftw functions work "in place", so we need to clone input data.
519     types::Double* D = nullptr;// A->clone()->getAs<types::Double>();
520     if (!func(name.data(), A, &D, way, gdim, opt))
521     {
522         FREE(gdim.dims);
523         FREE(gdim.howmany_dims);
524         FREE(Sel);
525         return types::Function::Error;
526     }
527
528     FREE(gdim.dims);
529     FREE(gdim.howmany_dims);
530     FREE(Sel);
531     out.push_back(D);
532     return types::Function::OK;
533 }
534
535 types::Function::ReturnValue common_4args(std::string& name, types::typed_list &in, int _iRetCount, types::typed_list &out, fftw_gen func, int way, int opt)
536 {
537     //getting the array A
538     types::Double* A = in[0]->getAs<types::Double>();
539     int ndimsA = A->getDims();
540     int *dimsA = A->getDimsArray();
541
542     /* Input  array variables */
543     int *Dim1 = nullptr;
544     int ndims = 0;
545     int *Incr = nullptr;
546     int nincr = 0;
547
548     /*FFTW specific library variable */
549     guru_dim_struct gdim = {0, NULL, 0, NULL};
550     /* input/output address for transform variables */
551
552     /* local variable */
553     int *Dim = nullptr, *Sel = nullptr;
554     int pd = 1;
555     int pds = 1;
556     int nd = 0;
557     int rank = 0;
558     int i = 0, j = 0, k = 0, lA = 1;
559
560     for (i = 0; i < ndimsA; i++)
561     {
562         lA *= dimsA[i];
563     }
564
565     /* void or scalar input gives void output or scalar*/
566     if (lA == 0)
567     {
568         out.push_back(A);
569         return types::Function::OK;
570     }
571
572     if (in[2]->isGenericType() == false)
573     {
574         Scierror(999, _("%s: Wrong type for input argument #%d.\n"), name.data(), 3);
575         FREE(gdim.dims);
576         FREE(gdim.howmany_dims);
577         return types::Function::Error;
578     }
579
580     getVarAsDims(in[2], ndims, Dim1);
581     /* check values of Dim1[i] */
582     pd = 1;
583     for (i = 0; i < ndims; i++)
584     {
585         if (Dim1[i] <= 1)
586         {
587             Scierror(999, _("%s: Wrong values for input argument #%d: Elements must be greater than %d.\n"), name.data(), 3, 1);
588             FREE(Dim1);
589             FREE(gdim.dims);
590             FREE(gdim.howmany_dims);
591             return types::Function::Error;
592         }
593         pd *= Dim1[i];
594     }
595     if (pd > lA)
596     {
597         Scierror(999, _("%s: Wrong values for input argument #%d: Must be less than %d.\n"), name.data(), 3, lA);
598         FREE(Dim1);
599         FREE(gdim.dims);
600         FREE(gdim.howmany_dims);
601         return types::Function::Error;
602     }
603     if (lA % pd)
604     {
605         Scierror(999, _("%s: Wrong values for input argument #%d: Must be a divisor of %d.\n"), name.data(), 3, lA);
606         FREE(Dim1);
607         FREE(gdim.dims);
608         FREE(gdim.howmany_dims);
609         return types::Function::Error;
610     }
611
612     if (in[3]->isGenericType() == false)
613     {
614         Scierror(999, _("%s: Wrong type for input argument #%d.\n"), name.data(), 4);
615         FREE(Dim1);
616         FREE(gdim.dims);
617         FREE(gdim.howmany_dims);
618         return types::Function::Error;
619     }
620
621     getVarAsDims(in[3], nincr, Incr);
622
623     if (nincr != ndims)
624     {
625         Scierror(999, _("%s: Incompatible input arguments #%d and #%d: Same sizes expected.\n"), name.data(), 3, 4);
626         FREE(Dim1);
627         FREE(Incr);
628         FREE(Dim);
629         FREE(Sel);
630         FREE(gdim.dims);
631         FREE(gdim.howmany_dims);
632         return types::Function::Error;
633     }
634
635     /* check values of Incr[i] */
636     if (Incr[0] <= 0)
637     {
638         Scierror(999, _("%s: Wrong values for input argument #%d: Positive integers expected.\n"), name.data(), 4);
639         FREE(Dim1);
640         FREE(Incr);
641         FREE(gdim.dims);
642         FREE(gdim.howmany_dims);
643         return types::Function::Error;
644     }
645     for (i = 0; i < ndims; i++)
646     {
647         if (lA % Incr[i])
648         {
649             Scierror(999, _("%s: Wrong values for input argument #%d: Elements must be divisors of %d.\n"), name.data(), 3, lA);
650             FREE(Dim1);
651             FREE(Incr);
652             FREE(gdim.dims);
653             FREE(gdim.howmany_dims);
654             return types::Function::Error;
655         }
656         if (i > 0 && (Incr[i] <= Incr[i - 1]))
657         {
658             Scierror(999, _("%s: Wrong values for input argument #%d: Elements must be in increasing ""order.\n"), name.data(), 4);
659             FREE(Dim1);
660             FREE(Incr);
661             FREE(gdim.dims);
662             FREE(gdim.howmany_dims);
663             return types::Function::Error;
664         }
665     }
666     if ((Dim = (int *)MALLOC((2 * ndims + 1) * sizeof(int))) == nullptr)
667     {
668         Scierror(999, _("%s: Cannot allocate more memory.\n"), name.data());
669         FREE(Dim1);
670         FREE(Incr);
671         FREE(Dim);
672         FREE(gdim.dims);
673         FREE(gdim.howmany_dims);
674         return types::Function::Error;
675     }
676     if ((Sel = (int *)MALLOC((ndims) * sizeof(int))) == nullptr)
677     {
678         Scierror(999, _("%s: Cannot allocate more memory.\n"), name.data());
679         FREE(Dim1);
680         FREE(Incr);
681         FREE(Dim);
682         FREE(Sel);
683         FREE(gdim.dims);
684         FREE(gdim.howmany_dims);
685         return types::Function::Error;
686     }
687
688
689     /*Transform  Dim1 and Incr into Dim and Sel and check validity*/
690
691     nd = 0;
692     pd = 1;
693     if (Incr[0] != 1)
694     {
695         Dim[nd++] = Incr[0];
696         pd *= Incr[0];
697     }
698     Dim[nd++] = Dim1[0];
699     pd *= Dim1[0];
700     Sel[0] = nd;
701
702     for (k = 1; k < ndims; k++)
703     {
704         if (Incr[k] % pd != 0)
705         {
706             Scierror(999, _("%s: Incompatible input arguments #%d and #%d.\n"), name.data(), 3, 4);
707             FREE(Dim1);
708             FREE(Incr);
709             FREE(Dim);
710             FREE(Sel);
711             FREE(gdim.dims);
712             FREE(gdim.howmany_dims);
713             return types::Function::Error;
714         }
715         if (Incr[k] != pd)
716         {
717             Dim[nd++] = (int)(Incr[k] / pd);
718             pd = Incr[k];
719         }
720         Dim[nd++] = Dim1[k];
721         pd *= Dim1[k];
722         Sel[k] = nd;
723     }
724     if (pd < lA)
725     {
726         if (lA % pd != 0)
727         {
728             Scierror(999, _("%s: Incompatible input arguments #%d and #%d.\n"), name.data(), 3, 4);
729             FREE(Dim1);
730             FREE(Incr);
731             FREE(Dim);
732             FREE(Sel);
733             FREE(gdim.dims);
734             FREE(gdim.howmany_dims);
735             return types::Function::Error;
736         }
737         Dim[nd++] = (int)(lA / pd);
738     }
739
740     rank = ndims;
741     ndims = nd;
742     /* now  same algorithm than sci_dct_3args applies */
743     /* Create  gdim struct */
744     gdim.rank = rank;
745     if ((gdim.dims = (fftw_iodim *)MALLOC(sizeof(fftw_iodim) * gdim.rank)) == nullptr)
746     {
747         Scierror(999, _("%s: Cannot allocate more memory.\n"), name.data());
748         FREE(Dim1);
749         FREE(Incr);
750         FREE(Dim);
751         FREE(Sel);
752         FREE(gdim.dims);
753         FREE(gdim.howmany_dims);
754         return types::Function::Error;
755     }
756
757     pd = 1; /* used to store prod(Dims(1:sel(k-1)))*/
758     pds = 1; /* used to store prod(Dims(sel(k-1):sel(k)))*/
759     j = 0;
760     for (i = 0; i < ndims; i++)
761     {
762         if (j >= gdim.rank)
763         {
764             break;
765         }
766         if (Sel[j] == i + 1)
767         {
768             gdim.dims[j].n = Dim[i];
769             gdim.dims[j].is = pd;
770             gdim.dims[j].os = pd;
771             j++;
772         }
773         pd *= Dim[i];
774     }
775     /* Compute howmany_rank based on jumps in the Sel sequence */
776     gdim.howmany_rank = 0;
777     if ((Sel[0] != 1) && (Sel[0] != ndims))
778     {
779         gdim.howmany_rank++;
780     }
781
782     for (i = 1; i <= rank - 1; i++)
783     {
784         if (Sel[i] != Sel[i - 1] + 1)
785         {
786             /*check if all dimensions between Sel[i-1]+1 and Sel[i]-1 are
787             equal to one, in this case they can be ignored and there is
788             no jump*/
789             for (j = Sel[i - 1] + 1; j <= Sel[i] - 1; j++)
790             {
791                 if (Dim[j - 1] != 1)
792                 {
793                     gdim.howmany_rank++;
794                     break;
795                 }
796             }
797         }
798     }
799     if ((Sel[rank - 1] != ndims) || (rank == 1))
800     {
801         gdim.howmany_rank++;
802     }
803     /* Fill the howmany_dims struct */
804     if (gdim.howmany_rank > 0)
805     {
806         /* it must be the case */
807         int ih = 0;
808
809         if ((gdim.howmany_dims = (fftw_iodim *)MALLOC(gdim.howmany_rank * sizeof(fftw_iodim))) == nullptr)
810         {
811             Scierror(999, _("%s: Cannot allocate more memory.\n"), name.data());
812             FREE(Dim1);
813             FREE(Incr);
814             FREE(Dim);
815             FREE(Sel);
816             FREE(gdim.dims);
817             FREE(gdim.howmany_dims);
818             return types::Function::Error;
819         }
820         pd = 1;
821         for (j = 1; j <= (Sel[0] - 1); j++)
822         {
823             pd *= Dim[j - 1];    /*prod(Dims(1:(sel(1)-1)))*/
824         }
825         ih = 0;
826         if ((Sel[0] != 1) && (Sel[0] != ndims))
827         {
828             /* First seleted dimension */
829             gdim.howmany_dims[ih].is = 1;
830             gdim.howmany_dims[ih].os = 1;
831             gdim.howmany_dims[ih].n = pd;
832             ih++;
833         }
834         pd *= Dim[Sel[0] - 1]; /*prod(Dims(1:sel(1)))*/
835
836         for (i = 1; i <= rank - 1; i++)
837         {
838             /* intermediate selected dimensions */
839             if (Sel[i] != Sel[i - 1] + 1)
840             {
841                 pds = 1;
842                 for (j = (Sel[i - 1] + 1); j <= (Sel[i] - 1); j++)
843                 {
844                     pds *= Dim[j - 1];    /*prod(Dims(sel(i-1)+1:(sel(i)-1)))*/
845                 }
846                 /*check again if all dimensions between Sel[i-1]+1 and
847                 Sel[i]-1 are equal to one, in this case they can be
848                 ignored and there is no jump*/
849                 for (j = (Sel[i - 1] + 1); j <= (Sel[i] - 1); j++)
850                 {
851                     if (Dim[j - 1] != 1)
852                     {
853                         gdim.howmany_dims[ih].is = pd;
854                         gdim.howmany_dims[ih].os = pd;
855                         gdim.howmany_dims[ih].n = pds;
856                         ih++;
857                         break;
858                     }
859                 }
860             }
861
862             pd *= pds * Dim[Sel[i] - 1]; /*prod(Dims(1:sel(i)))*/
863         }
864
865         if (Sel[rank - 1] != ndims)
866         {
867             /* last selected dimension*/
868             pds = 1;
869             for (j = (Sel[rank - 1] + 1); j <= ndims; j++)
870             {
871                 pds *= Dim[j - 1];    /*prod(Dims(sel(i-1)+1:(sel(i)-1)))*/
872             }
873             gdim.howmany_dims[ih].is = pd;
874             gdim.howmany_dims[ih].os = pd;
875             gdim.howmany_dims[ih].n = pds;
876             ih++;
877         }
878         else if (rank == 1) /* the only selected dimension is the last one */
879         {
880             gdim.howmany_dims[ih].is = 1;
881             gdim.howmany_dims[ih].os = 1;
882             gdim.howmany_dims[ih].n = pd / Dim[Sel[0] - 1];
883             ih++;
884         }
885     }
886
887     //fftw functions work "in place", so we need to clone input data.
888     types::Double* D = nullptr;// A->clone()->getAs<types::Double>();
889     if (!func(name.data(), A, &D, way, gdim, opt))
890     {
891         FREE(Dim1);
892         FREE(Incr);
893         FREE(Dim);
894         FREE(Sel);
895         FREE(gdim.dims);
896         FREE(gdim.howmany_dims);
897         return types::Function::Error;
898     }
899
900     FREE(Dim1);
901     FREE(Incr);
902     FREE(Dim);
903     FREE(Sel);
904     FREE(gdim.dims);
905     FREE(gdim.howmany_dims);
906
907     out.push_back(D);
908     return types::Function::OK;
909 }
910
911 void getVarAsDims(types::InternalType* t, int& rank, int*& Sel)
912 {
913     switch (t->getType())
914     {
915         case types::InternalType::ScilabDouble:
916         {
917             getVarAsDims(t->getAs<types::Double>(), rank, Sel);
918             break;
919         }
920         case types::InternalType::ScilabInt8:
921         {
922             getVarAsDims(t->getAs<types::Int8>(), rank, Sel);
923             break;
924         }
925         case types::InternalType::ScilabInt16:
926         {
927             getVarAsDims(t->getAs<types::Int16>(), rank, Sel);
928             break;
929         }
930         case types::InternalType::ScilabInt32:
931         {
932             getVarAsDims(t->getAs<types::Int32>(), rank, Sel);
933             break;
934         }
935         case types::InternalType::ScilabInt64:
936         {
937             getVarAsDims(t->getAs<types::Int64>(), rank, Sel);
938             break;
939         }
940         case types::InternalType::ScilabUInt8:
941         {
942             getVarAsDims(t->getAs<types::UInt8>(), rank, Sel);
943             break;
944         }
945         case types::InternalType::ScilabUInt16:
946         {
947             getVarAsDims(t->getAs<types::UInt16>(), rank, Sel);
948             break;
949         }
950         case types::InternalType::ScilabUInt32:
951         {
952             getVarAsDims(t->getAs<types::UInt32>(), rank, Sel);
953             break;
954         }
955         case types::InternalType::ScilabUInt64:
956         {
957             getVarAsDims(t->getAs<types::UInt64>(), rank, Sel);
958             break;
959         }
960     }
961 }
962
963 template<class T>
964 void getVarAsDims(T* t, int& dims, int*& pdims)
965 {
966     dims = t->getSize();
967     pdims = new int[dims];
968     for (int i = 0; i < dims; ++i)
969     {
970         pdims[i] = static_cast<int>(t->get(i));
971     }
972 }