Merge remote-tracking branch 'origin/master' into YaSp
[scilab.git] / scilab / modules / fftw / sci_gateway / c / sci_fftw.c
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2006/2007 - INRIA - Alan LAYEC
4 * Copyright (C) 2007 - INRIA - Allan CORNET
5 * Copyright (C) 2012 - DIGITEO - Allan CORNET
6 * Copyright (C) 2012 - INRIA - Serge STEER
7 *
8 * This file must be used under the terms of the CeCILL.
9 * This source file is licensed as described in the file COPYING, which
10 * you should have received as part of this distribution.  The terms
11 * are also available at
12 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
13 *
14 */
15 /*--------------------------------------------------------------------------*/
16 #include "stack-c.h"
17 #include "fftw_utilities.h"
18 #include "callfftw.h"
19 #include "MALLOC.h"
20 #include "gw_fftw.h"
21 #include "api_scilab.h"
22 #include "localization.h"
23 #include "Scierror.h"
24 #include "BOOL.h"
25 /*--------------------------------------------------------------------------*/
26 enum Scaling
27 {
28     Divide = -1,
29     None = 0,
30     Multiply = 1,
31 };
32 /*--------------------------------------------------------------------------*/
33 extern void C2F(dscal)(int *n, double *da, double *dx, int *incx); /* blas routine */
34 /*--------------------------------------------------------------------------*/
35 static int getArrayOfDouble(void* _pvCtx, int *piAddr, int *ndims, int **dims, double **Ar, double **Ai);
36 static SciErr allocArrayOfDouble(void* _pvCtx, int _iVar, int ndims, int *dims, double **Ar);
37 static SciErr allocComplexArrayOfDouble(void* _pvCtx, int _iVar, int ndims, int *dims, double **Ar, double **Ai);
38 static SciErr getScalarIntArg(void* _pvCtx, int _iVar, char *fname, int *value);
39 static SciErr getVectorIntArg(void* _pvCtx, int _iVar, char *fname, int *pndims, int **pDim);
40 static BOOL isHyperMatrixMlist(void* _pvCtx, int *piAddressVar);
41 static int sci_fft_gen(void* _pvCtx, char *fname, int ndimsA, int *dimsA, double *Ar,  double *Ai, int isn, int iopt, guru_dim_struct gdim);
42 static int sci_fft_2args(void* _pvCtx, char *fname, int ndimsA, int *dimsA, double *Ar,  double *Ai, int isn, int iopt);
43 static int sci_fft_3args(void* _pvCtx, char *fname, int ndimsA, int *dimsA, double *Ar,  double *Ai, int isn, int iopt);
44 static int sci_fft_4args(void* _pvCtx, char *fname, int ndimsA, int *dimsA, double *Ar,  double *Ai, int isn, int iopt);
45 /*--------------------------------------------------------------------------*/
46
47 int WITHMKL = 0;
48 /* fftw function.
49 *
50 * Scilab Calling sequence :
51 *   fftw(A [,option])
52 *   fftw(A,sign [,option])
53 *   fftw(A,sel,sign [,option])
54 *   fftw(A,sign,dim,incr [,option])
55 *
56 * Input : A : a scilab double complex or real vector, matrix or hypermatrix
57 *
58 *         sign : a scilab double or integer scalar (-1 or 1): the sign
59 *                  in the exponential component
60 *
61 *         sel : a scilab double or integer vector, the selection of dimensions
62
63 *         dim : a scilab double or integer vector: the dimensions
64 *                  of the Fast Fourier Transform to perform
65 *
66 *         incr : a scilab double or integer vector: the increments
67 *                  of the Fast Fourier Transform to perform
68 *
69 * Output : a scilab double complex or real array with same shape as A that
70 *          gives the result of the transform.
71 *
72 */
73 int sci_fftw(char *fname, void* pvApiCtx)
74 {
75     SciErr sciErr;
76     int *piAddr = NULL;
77     char *option = NULL;
78     int iopt = 0; /* automatic r2c or c2r transform use decision */
79     int rhs = Rhs;
80     int iTypeOne = 0;
81
82     int ndimsA = 0;
83     int *dimsA = NULL;
84     double *Ar = NULL, *Ai = NULL;
85
86     int isn = FFTW_FORWARD;
87     WITHMKL = withMKL();
88     /****************************************
89     * Basic constraints on rhs arguments  *
90     ****************************************/
91
92     /* check min/max lhs/rhs arguments of scilab function */
93     CheckRhs(1, 5);
94     CheckLhs(1, 1);
95
96     sciErr = getVarAddressFromPosition(pvApiCtx, 1, &piAddr);
97     if (sciErr.iErr)
98     {
99         printError(&sciErr, 0);
100         Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, 1);
101         return 0;
102     }
103
104     sciErr = getVarType(pvApiCtx, piAddr, &iTypeOne);
105     if (sciErr.iErr)
106     {
107         printError(&sciErr, 0);
108         Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, 1);
109         return 0;
110     }
111
112     if ((iTypeOne == sci_list) ||
113         (iTypeOne == sci_tlist))
114     {
115         OverLoad(1);
116         return 0;
117     }
118
119     if (iTypeOne == sci_mlist)
120     {
121         /* We allow overload for not hypermatrix type */
122         if (!isHyperMatrixMlist(pvApiCtx, piAddr))
123         {
124             OverLoad(1);
125             return 0;
126         }
127     }
128
129     /* checking if last argument is a potential option argument (character string) */
130     sciErr = getVarAddressFromPosition(pvApiCtx, Rhs, &piAddr);
131     if (sciErr.iErr)
132     {
133         printError(&sciErr, 0);
134         Scierror(999, _("%s: Can not read input argument #%d.\n"), fname, Rhs);
135         return 0;
136     }
137
138     if (isStringType(pvApiCtx, piAddr))   /*  fftw(...,option); */
139     {
140         if (isScalar(pvApiCtx, piAddr))
141         {
142             if (getAllocatedSingleString(pvApiCtx, piAddr, &option) == 0)
143             {
144                 if (strcmp("symmetric", option) == 0)  iopt = 1; /*user assumes symmetry */
145                 else if (strcmp("nonsymmetric", option) == 0) iopt = 2; /*user claims full transform */
146                 else
147                 {
148                     Scierror(999, _("%s: Wrong value for input argument #%d: '%s' or '%s' expected.\n"),
149                         fname, Rhs, "\"symmetric\"", "\"nonsymmetric\"");
150                     freeAllocatedSingleString(option);
151                     option = NULL;
152                     return 0;
153                 }
154                 freeAllocatedSingleString(option);
155                 option = NULL;
156                 rhs = Rhs - 1;
157             }
158             else
159             {
160                 Scierror(999, _("%s: Wrong value for input argument #%d: '%s' or '%s' expected.\n"),
161                     fname, Rhs, "\"full\", \"same\"", "\"valid\"");
162                 return 0;
163             }
164         }
165     }
166
167
168     /********************  Checking if isn is given  ************************************************/
169     if (rhs == 1)  /*only one rhs argument: forward fft*/
170     {
171         isn = FFTW_FORWARD; /* default value */
172     }
173     else   /*get isn out of second argument*/
174     {
175         sciErr = getScalarIntArg(pvApiCtx, 2, fname, &isn);
176         if (sciErr.iErr)
177         {
178             Scierror(sciErr.iErr, getErrorMessage(sciErr));
179             return 0;
180         }
181         /* check value of second rhs argument */
182         if ((isn !=  FFTW_FORWARD) && (isn !=  FFTW_BACKWARD))
183         {
184             Scierror(53, _("%s: Wrong value for input argument #%d: %d or %d expected.\n"),
185                 fname, 2, FFTW_FORWARD, FFTW_BACKWARD);
186             return(0);
187         }
188     }
189
190     /********************  getting the array A      ************************************************/
191     getVarAddressFromPosition(pvApiCtx, 1, &piAddr);
192     if (!getArrayOfDouble(pvApiCtx, piAddr, &ndimsA, &dimsA, &Ar, &Ai))
193     {
194         Scierror(999, _("%s: Wrong type for argument #%d: Array of floating point numbers expected.\n"),
195             fname, 1);
196         return 0;
197     }
198
199
200     /********************  Select proper method     ************************************************/
201     if (rhs < 3)
202     {
203         /* fftw(A ,sign [,option])*/
204         sci_fft_2args(pvApiCtx, fname, ndimsA, dimsA, Ar, Ai, isn, iopt);
205     }
206     else if (rhs == 3)
207     {
208         /* fftw(A ,sign ,sel [,option])*/
209         sci_fft_3args(pvApiCtx, fname, ndimsA, dimsA, Ar, Ai, isn, iopt);
210     }
211     else if (rhs == 4)
212     {
213         /* fftw(A ,sign ,dim,incr [option])*/
214         sci_fft_4args(pvApiCtx, fname, ndimsA, dimsA, Ar, Ai, isn, iopt);
215     }
216
217     return 0;
218 }
219
220 /*--------------------------------------------------------------------------*/
221 int getArrayOfDouble(void* _pvCtx, int *piAddr, int *ndims, int **dims, double **Ar, double **Ai)
222 {
223     SciErr sciErr;
224     int *piAddrChild = NULL;
225     int *piOffset = NULL;
226     int *piData = NULL;
227     int nItems = 0;
228     int iRows = 0;
229     int iCols = 0;
230     int iType = 0;
231
232     sciErr = getVarType(_pvCtx, piAddr, &iType);
233     if (iType == sci_matrix)
234     {
235         *ndims = 2;
236         *dims = &(piAddr[1]);
237         if (isVarComplex(_pvCtx, piAddr))
238         {
239             getComplexMatrixOfDouble(_pvCtx, piAddr, &iRows, &iCols, Ar, Ai);
240         }
241         else
242         {
243             getMatrixOfDouble(_pvCtx, piAddr, &iRows, &iCols, Ar);
244             *Ai = NULL;
245         }
246         return 1;
247     }
248     else if (iType == sci_mlist)
249     {
250         sciErr = getListItemNumber(_pvCtx, piAddr, &nItems);
251         if (nItems != 3) return 0;
252         /*Check if first item is ["hm","dims","entries"] */
253         sciErr = getListItemAddress(_pvCtx, piAddr, 1, &piAddrChild);
254         sciErr = getVarType(_pvCtx, piAddrChild, &iType);
255         if (iType != sci_strings) return 0;
256         sciErr = getVarDimension(_pvCtx, piAddrChild, &iRows, &iCols);
257         if (iRows*iCols != 3) return 0;
258         /* Check if first entry of the first item is "hm" */
259         piOffset = piAddrChild + 4;
260         if (piOffset[1] - piOffset[0] != 2)  return 0;
261         piData = piOffset + iRows * iCols + 1;
262         if (piData[0] != 17 || piData[1] != 22) return 0; /* check "hm" */
263         /* Get second item dims */
264         sciErr = getListItemAddress(_pvCtx, piAddr, 2, &piAddrChild);
265         sciErr = getVarType(_pvCtx, piAddrChild, &iType);
266         if (iType != sci_ints) return 0;
267         sciErr = getMatrixOfInteger32(_pvCtx, piAddrChild, &iRows, &iCols, dims);
268         if (sciErr.iErr)  return 0;
269         *ndims = iRows * iCols;
270         /* Get thirds item entries */
271         sciErr = getListItemAddress(_pvCtx, piAddr, 3, &piAddrChild);
272         sciErr = getVarType(_pvCtx, piAddrChild, &iType);
273         if (iType != sci_matrix) return 0;
274         if (isVarComplex(_pvCtx, piAddrChild))
275         {
276             getComplexMatrixOfDouble(_pvCtx, piAddrChild, &iRows, &iCols, Ar, Ai);
277         }
278         else
279         {
280             getMatrixOfDouble(_pvCtx, piAddrChild, &iRows, &iCols, Ar);
281             *Ai = NULL;
282         }
283         return 1;
284     }
285     else
286     {
287         return 0;
288     }
289 }
290
291 SciErr allocComplexArrayOfDouble(void* _pvCtx, int _iVar, int ndims, int *dims, double **Ar, double **Ai)
292 {
293     SciErr sciErr;
294     int *piAddr = NULL;
295
296     if (ndims == 2)
297     {
298         sciErr = allocComplexMatrixOfDouble( _pvCtx, _iVar, dims[0], dims[1], Ar, Ai);
299         if (sciErr.iErr) return sciErr;
300     }
301     else
302     {
303         int i = 0;
304         int n = 1;
305         const char * hmType[] = {"hm", "dims", "entries"};
306
307         for (i = 0; i < ndims; i++)   n *= dims[i];
308
309         sciErr = createMList(_pvCtx, _iVar, 3, &piAddr);
310         if (sciErr.iErr) return sciErr;
311
312         sciErr = createMatrixOfStringInList(_pvCtx, _iVar, piAddr, 1, 1, 3, hmType);
313         if (sciErr.iErr) return sciErr;
314         sciErr = createMatrixOfInteger32InList(_pvCtx, _iVar, piAddr, 2, 1, ndims, dims);
315         if (sciErr.iErr) return sciErr;
316
317         sciErr = allocComplexMatrixOfDoubleInList(_pvCtx, _iVar, piAddr, 3, n, 1, Ar, Ai);
318         if (sciErr.iErr) return sciErr;
319     }
320     return sciErr;
321 }
322
323 SciErr allocArrayOfDouble(void* _pvCtx, int _iVar,  int ndims, int *dims, double **Ar)
324 {
325     SciErr sciErr;
326     int *piAddr      = NULL;
327
328
329     if (ndims == 2)
330     {
331         sciErr = allocMatrixOfDouble( _pvCtx, _iVar, dims[0], dims[1], Ar);
332         if (sciErr.iErr) return sciErr;
333     }
334     else
335     {
336         int i = 0;
337         int n = 1;
338         const char * hmType[] = {"hm", "dims", "entries"};
339
340         for (i = 0; i < ndims; i++) n *= dims[i];
341
342         sciErr = createMList(_pvCtx,  _iVar, 3, &piAddr);
343         if (sciErr.iErr) return sciErr;
344         sciErr = createMatrixOfStringInList(_pvCtx, _iVar, piAddr, 1, 1, 3, hmType);
345         if (sciErr.iErr) return sciErr;
346         sciErr = createMatrixOfInteger32InList(_pvCtx, _iVar, piAddr, 2, 1, ndims, dims);
347         if (sciErr.iErr) return sciErr;
348
349         sciErr = allocMatrixOfDoubleInList(_pvCtx, _iVar, piAddr, 3, n, 1, Ar);
350         if (sciErr.iErr) return sciErr;
351     }
352     return sciErr;
353 }
354
355 SciErr getScalarIntArg(void* _pvCtx, int _iVar, char *fname, int *value)
356 {
357     SciErr sciErr;
358     int *piAddr = NULL;
359     int iType = 0;
360     int iPrec = 0;
361     double t_d = 0.0;
362     char t_c = 0;
363     unsigned char t_uc = 0;
364     short t_s = 0;
365     unsigned short t_us = 0;
366     int t_i = 0;
367     unsigned int  t_ui = 0;
368     sciErr.iErr = 0;
369     sciErr.iMsgCount = 0;
370
371     sciErr = getVarAddressFromPosition(_pvCtx, _iVar, &piAddr);
372     if (sciErr.iErr)
373     {
374         addErrorMessage(&sciErr, API_ERROR_GET_STRING,  _("%s: Can not read input argument #%d.\n"), fname, _iVar);
375         return sciErr;
376     }
377
378     //check type
379     sciErr = getVarType(_pvCtx, piAddr, &iType);
380     if (sciErr.iErr)
381     {
382         addErrorMessage(&sciErr, API_ERROR_GET_INT,  _("%s: Can not read input argument #%d.\n"), fname, _iVar);
383         return sciErr;
384     }
385
386     if (!isScalar(_pvCtx, piAddr))
387     {
388         addErrorMessage(&sciErr, API_ERROR_GET_INT, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), fname, _iVar);
389         return sciErr;
390     }
391
392     if (iType == sci_matrix)
393     {
394         getScalarDouble(_pvCtx, piAddr, &t_d);
395         *value = (int)t_d;
396     }
397     else if (iType == sci_ints)
398     {
399         sciErr = getMatrixOfIntegerPrecision(_pvCtx, piAddr, &iPrec);
400         if (sciErr.iErr)
401         {
402             addErrorMessage(&sciErr, API_ERROR_GET_INT, _("%s: Can not read input argument #%d.\n"), fname, _iVar);
403             return sciErr;
404         }
405
406         switch (iPrec)
407         {
408         case SCI_INT8 :
409             {
410                 getScalarInteger8(_pvCtx, piAddr, &t_c);
411                 *value = (int)t_c;
412             }
413         case SCI_INT16 :
414             {
415                 getScalarInteger16(_pvCtx, piAddr, &t_s);
416                 *value = (int)t_s;
417             }
418         case SCI_INT32 :
419             {
420                 getScalarInteger32(_pvCtx, piAddr, &t_i);
421                 *value = (int)t_i;
422             }
423         case SCI_UINT8 :
424             {
425                 getScalarUnsignedInteger8(_pvCtx, piAddr, &t_uc);
426                 *value = (int)t_uc;
427             }
428         case SCI_UINT16 :
429             {
430                 getScalarUnsignedInteger16(_pvCtx, piAddr, &t_us);
431                 *value = (int)t_us;
432             }
433         case SCI_UINT32 :
434             {
435                 getScalarUnsignedInteger32(_pvCtx, piAddr, &t_ui);
436                 *value = (int)t_ui;
437             }
438         }
439     }
440     else
441     {
442         addErrorMessage(&sciErr, API_ERROR_GET_INT,
443             _("%s: Wrong type for argument #%d: An integer or a floating point number expected.\n"),
444             fname, _iVar);
445         return sciErr;
446     }
447     return sciErr;
448 }
449
450 SciErr getVectorIntArg(void* _pvCtx, int _iVar, char *fname, int *pndims, int **pDim)
451 {
452     SciErr sciErr;
453     int *piAddr = NULL;
454     int iType = 0;
455     int iPrec = 0;
456     int mDim = 0;
457     int nDim = 0;
458     int *Dim = NULL;
459     int ndims = 0;
460
461     double* p_d = NULL;
462     char* p_c = NULL;
463     unsigned char* p_uc = NULL;
464     short* p_s = NULL;
465     unsigned short* p_us = NULL;
466     int* p_i = NULL;
467     unsigned int*  p_ui = NULL;
468     int i = 0;
469
470     sciErr.iErr = 0;
471     sciErr.iMsgCount = 0;
472
473     getVarAddressFromPosition(_pvCtx, _iVar, &piAddr);
474
475     //check type
476     getVarType(_pvCtx, piAddr, &iType);
477
478     if (isVarMatrixType(_pvCtx, piAddr) == 0)
479     {
480         addErrorMessage(&sciErr, API_ERROR_GET_INT, _("%s: Wrong type for input argument #%d.\n"), fname, _iVar);
481         return sciErr;
482     }
483
484     getVarDimension(_pvCtx, piAddr, &mDim, &nDim);
485
486     ndims = mDim * nDim;
487     *pndims = ndims;
488     if (ndims <= 0)
489     {
490         addErrorMessage(&sciErr, API_ERROR_GET_INT,
491             _("%s: Wrong size for input argument #%d.\n"), fname, _iVar);
492         return sciErr;
493     }
494     if ((Dim = (int *)MALLOC(ndims * sizeof(int))) == NULL)
495     {
496         addErrorMessage(&sciErr, API_ERROR_GET_INT,
497             _("%s: Cannot allocate more memory.\n"), fname);
498         return sciErr;
499     }
500     *pDim = Dim;
501     if (iType == sci_matrix)
502     {
503         sciErr = getMatrixOfDouble(_pvCtx, piAddr, &mDim, &nDim, &p_d);
504         for (i = 0; i < ndims; i++)  Dim[i] = (int)(p_d[i]);
505     }
506     else if (iType == sci_ints)
507     {
508         getMatrixOfIntegerPrecision(_pvCtx, piAddr, &iPrec);
509         switch (iPrec)
510         {
511         case SCI_INT8 :
512             getMatrixOfInteger8(_pvCtx, piAddr, &mDim, &nDim, &p_c);
513             for (i = 0; i < ndims; i++) Dim[i]  = (int)(p_c[i]);
514             break;
515         case SCI_INT16 :
516             getMatrixOfInteger16(_pvCtx, piAddr, &mDim, &nDim, &p_s);
517             for (i = 0; i < ndims; i++) Dim[i]  = (int)(p_s[i]);
518             break;
519         case SCI_INT32 :
520             getMatrixOfInteger32(_pvCtx, piAddr, &mDim, &nDim, &p_i);
521             for (i = 0; i < ndims; i++)  Dim[i]  = (int)(p_i[i]);
522             break;
523         case SCI_UINT8 :
524             getMatrixOfUnsignedInteger8(_pvCtx, piAddr, &mDim, &nDim, &p_uc);
525             for (i = 0; i < ndims; i++) Dim[i]  = (int)(p_uc[i]);
526             break;
527         case SCI_UINT16 :
528             getMatrixOfUnsignedInteger16(_pvCtx, piAddr, &mDim, &nDim, &p_us);
529             for (i = 0; i < ndims; i++) Dim[i]  = (int) p_us[i];
530             break;
531         case SCI_UINT32 :
532             getMatrixOfUnsignedInteger32(_pvCtx, piAddr, &mDim, &nDim, &p_ui);
533             for (i = 0; i < ndims; i++) Dim[i]  = (int)(p_ui[i]);
534             break;
535         }
536     }
537     else
538     {
539         FREE(Dim);
540         Dim = NULL;
541         addErrorMessage(&sciErr, API_ERROR_GET_INT,
542             _("%s: Wrong type for argument #%d: An array of floating point or integer numbers expected.\n"), fname, _iVar);
543         return sciErr;
544     }
545     return sciErr;
546 }
547
548 int sci_fft_2args(void* _pvCtx, char *fname, int ndimsA, int *dimsA, double *Ar,  double *Ai, int isn, int iopt)
549 {
550     /*FFTW specific library variable */
551     guru_dim_struct gdim = {0, NULL, 0, NULL};
552
553     /* local variable */
554     int ndims = 0; /* number of non singleton dimensions */
555     int first_nonsingleton = -1;
556     int i = 0, j = 0;
557     int prd = 1;
558
559     /* ignore singleton dimensions */
560     first_nonsingleton = -1;
561     ndims = 0;
562     for (i = 0; i < ndimsA; i++)
563     {
564         if (dimsA[i] > 1)
565         {
566             ndims++;
567             if (first_nonsingleton < 0) first_nonsingleton = i;
568         }
569     }
570
571     /* void or scalar input gives void output or scalar*/
572     if (ndims == 0 )
573     {
574         AssignOutputVariable(_pvCtx, 1) =  1;
575         ReturnArguments(_pvCtx);
576         return(0);
577     }
578
579     gdim.rank = ndims;
580     if ((gdim.dims = (fftw_iodim *)MALLOC(sizeof(fftw_iodim) * gdim.rank)) == NULL)
581     {
582         Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
583         goto ERR;
584     }
585
586     j = 0;
587     prd = 1;
588     for (i = (first_nonsingleton); i < ndimsA; i++)
589     {
590         if (dimsA[i] > 1)
591         {
592             gdim.dims[j].n = dimsA[i];
593             gdim.dims[j].is = prd;
594             gdim.dims[j].os = prd;
595             prd *= dimsA[i];
596             j++;
597         }
598     }
599     gdim.howmany_rank = 0;
600     gdim.howmany_dims = NULL;
601
602
603     if (!sci_fft_gen(_pvCtx, fname, ndimsA, dimsA,  Ar,  Ai, isn, iopt, gdim))  goto ERR;
604
605
606     /***********************************
607     * Return results in lhs argument *
608     ***********************************/
609
610     ReturnArguments(_pvCtx);
611 ERR:
612     FREE(gdim.dims);
613     FREE(gdim.howmany_dims);
614     return(0);
615 }
616
617
618 int  sci_fft_3args(void* _pvCtx, char *fname, int ndimsA, int *dimsA, double *Ar,  double *Ai, int isn, int iopt)
619 {
620     /* API variables */
621     SciErr sciErr;
622     int *piAddr = NULL;
623
624     int *Sel = NULL;
625     int rank = 0;
626
627     /*FFTW specific library variable */
628     guru_dim_struct gdim = {0, NULL, 0, NULL};
629     /* local variable */
630     int ndims = 0;
631     int first_nonsingleton = -1;
632     int ih = 0;
633     int pd = 1; /* used to store prod(Dims(1:sel(k-1)))*/
634     int pds = 1; /* used to store prod(Dims(sel(k-1):sel(k)))*/
635     int i = 0, j = 0;
636
637     /* ignore singleton dimensions */
638     first_nonsingleton = -1;
639     ndims = 0;
640     for (i = 0; i < ndimsA; i++)
641     {
642         if (dimsA[i] > 1)
643         {
644             ndims++;
645             if (first_nonsingleton < 0) first_nonsingleton = i;
646         }
647     }
648
649     /* void or scalar input gives void output or scalar*/
650     if (ndims == 0 )
651     {
652         AssignOutputVariable(_pvCtx, 1) =  1;
653         ReturnArguments(_pvCtx);
654         return(0);
655     }
656
657
658     /******************** get and check third argument (sel) ****************************************/
659     getVarAddressFromPosition(_pvCtx, 3, &piAddr);
660     if (isVarMatrixType(_pvCtx, piAddr) == 0)
661     {
662         Scierror(999, _("%s: Wrong type for input argument #%d.\n"), fname, 3);
663         goto ERR;
664     }
665     sciErr = getVectorIntArg(_pvCtx, 3, fname, &rank, &Sel);
666     if (sciErr.iErr)
667     {
668         Scierror(sciErr.iErr, getErrorMessage(sciErr));
669         goto ERR;
670     }
671     /* size of Sel must be less than ndimsA */
672     if (rank <= 0 || rank >= ndimsA)
673     {
674         Scierror(999, _("%s: Wrong size for input argument #%d: Must be between %d and %d.\n"), fname, 3, 1, ndimsA - 1);
675         goto ERR;
676     }
677     /* check values of Sel[i] */
678     for (i = 0; i < rank; i++)
679     {
680         if (Sel[i] <= 0)
681         {
682             Scierror(999, _("%s: Wrong values for input argument #%d: Positive integers expected.\n"), fname, 3);
683             goto ERR;
684         }
685         if (Sel[i] > ndimsA)
686         {
687             Scierror(999, _("%s: Wrong values for input argument #%d: Elements must be less than %d.\n"), fname, 3, ndimsA);
688             goto ERR;
689         }
690         if (i > 0 && Sel[i] <= Sel[i - 1])
691         {
692             Scierror(999, _("%s: Wrong values for input argument #%d: Elements must be in increasing order.\n"), fname, 3);
693             goto ERR;
694         }
695     }
696
697     /* Create  gdim struct */
698     gdim.rank = rank;
699     if ((gdim.dims = (fftw_iodim *)MALLOC(sizeof(fftw_iodim) * gdim.rank)) == NULL)
700     {
701         Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
702         goto ERR;
703     }
704
705     pd = 1; /* used to store prod(Dims(1:sel(k-1)))*/
706     pds = 1; /* used to store prod(Dims(sel(k-1):sel(k)))*/
707     j = 0;
708     for (i = 0; i < ndimsA; i++)
709     {
710         if (j >= gdim.rank) break;
711         if (Sel[j] == i + 1)
712         {
713             gdim.dims[j].n = dimsA[i];
714             gdim.dims[j].is = pd;
715             gdim.dims[j].os = pd;
716             j++;
717         }
718         pd *= dimsA[i];
719     }
720     /* Compute howmany_rank based on jumps in the Sel sequence */
721     gdim.howmany_rank = 0;
722     if ((Sel[0] != 1) && (Sel[0] != ndimsA)) gdim.howmany_rank++;
723
724     for (i = 1; i <= rank - 1; i++)
725         if (Sel[i] != Sel[i - 1] + 1) gdim.howmany_rank++;
726
727     if ((Sel[rank - 1] != ndimsA) || (rank == 1)) gdim.howmany_rank++;
728
729     /* Fill the howmany_dims struct */
730     if (gdim.howmany_rank > 0)
731     {
732         /* it must be the case */
733         if ((gdim.howmany_dims = (fftw_iodim *)MALLOC(gdim.howmany_rank * sizeof(fftw_iodim))) == NULL)
734         {
735             Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
736             goto ERR;
737         }
738         pd = 1;
739         for (j = 1; j <= (Sel[0] - 1); j++) pd *= dimsA[j - 1]; /*prod(Dims(1:(sel(1)-1)))*/
740         ih = 0;
741         if ((Sel[0] != 1) && (Sel[0] != ndimsA))
742         {
743             /* First seleted dimension */
744             gdim.howmany_dims[ih].is = 1;
745             gdim.howmany_dims[ih].os = 1;
746             gdim.howmany_dims[ih].n = pd;
747             ih++;
748         }
749         pd *= dimsA[Sel[0] - 1]; /*prod(Dims(1:sel(1)))*/
750         for (i = 2; i <= rank; i++)
751         {
752             /* intermediate selected dimensions */
753             if (Sel[i - 1] != Sel[i - 2] + 1)
754             {
755                 pds = 1;
756                 for (j = (Sel[i - 2] + 1); j <= (Sel[i - 1] - 1); j++) pds *= dimsA[j - 1]; /*prod(Dims(sel(i-1)+1:(sel(i)-1)))*/
757                 gdim.howmany_dims[ih].is = pd;
758                 gdim.howmany_dims[ih].os = pd;
759                 gdim.howmany_dims[ih].n = pds;
760                 ih++;
761             }
762             pd *= pds * dimsA[Sel[i - 1] - 1]; /*prod(Dims(1:sel(i)))*/
763         }
764
765         if (Sel[rank - 1] != ndimsA)
766         {
767             /* last selected dimension*/
768             pds = 1;
769             for (j = (Sel[rank - 1] + 1); j <= ndimsA; j++) pds *= dimsA[j - 1]; /*prod(Dims(sel(i-1)+1:(sel(i)-1)))*/
770             gdim.howmany_dims[ih].is = pd;
771             gdim.howmany_dims[ih].os = pd;
772             gdim.howmany_dims[ih].n = pds;
773             ih++;
774         }
775         else if (rank == 1)
776         {
777             /* the only selected dimension is the last one */
778             gdim.howmany_dims[ih].is = 1;
779             gdim.howmany_dims[ih].os = 1;
780             gdim.howmany_dims[ih].n = pd / dimsA[Sel[0] - 1];
781             ih++;
782         }
783     }
784
785     if (!sci_fft_gen(_pvCtx, fname, ndimsA, dimsA, Ar,  Ai, isn, iopt, gdim))  goto ERR;
786     /***********************************
787     * Return results in lhs argument *
788     ***********************************/
789
790     ReturnArguments(_pvCtx);
791 ERR:
792     FREE(gdim.dims);
793     FREE(gdim.howmany_dims);
794     return(0);
795 }
796
797 int sci_fft_4args(void* _pvCtx, char *fname, int ndimsA, int *dimsA, double *Ar,  double *Ai, int isn, int iopt)
798 {
799     /* API variables */
800     SciErr sciErr;
801     int *piAddr = NULL;
802
803     /* Input  array variables */
804     int *Dim1 = NULL;
805     int ndims = 0;
806     int *Incr = NULL;
807     int nincr = 0;
808
809     /*FFTW specific library variable */
810     guru_dim_struct gdim = {0, NULL, 0, NULL};
811     /* input/output address for transform variables */
812
813     /* local variable */
814     int *Dim = NULL, *Sel = NULL;
815     int pd = 1;
816     int pds = 1;
817     int nd = 0;
818     int rank = 0;
819     int i = 0, j = 0, k = 0, lA = 1;
820
821     for (i = 0; i < ndimsA; i++)
822     {
823         lA *= dimsA[i];
824     }
825
826     /* void or scalar input gives void output or scalar*/
827     if (lA <= 1 )
828     {
829         AssignOutputVariable(_pvCtx, 1) =  1;
830         ReturnArguments(_pvCtx);
831         return(0);
832     }
833
834     /******************** get and check third argument (dim) ****************************************/
835     getVarAddressFromPosition(_pvCtx, 3, &piAddr);
836     if (isVarMatrixType(_pvCtx, piAddr) == 0)
837     {
838         Scierror(999, _("%s: Wrong type for input argument #%d.\n"), fname, 3);
839         goto ERR;
840     }
841     sciErr = getVectorIntArg(_pvCtx, 3, fname, &ndims, &Dim1);
842     if (sciErr.iErr)
843     {
844         Scierror(sciErr.iErr, getErrorMessage(sciErr));
845         goto ERR;
846     }
847     /* check values of Dim1[i} */
848     pd = 1;
849     for (i = 0; i < ndims; i++)
850     {
851         if (Dim1[i] <= 1)
852         {
853             Scierror(999, _("%s: Wrong values for input argument #%d: Elements must be greater than %d.\n"), fname, 3, 1);
854             goto ERR;
855         }
856         pd *= Dim1[i];
857     }
858     if ( pd > lA)
859     {
860         Scierror(999, _("%s: Wrong values for input argument #%d: Must be less than %d.\n"), fname, 3, lA);
861         goto ERR;
862     }
863     if (lA % pd)
864     {
865         Scierror(999, _("%s: Wrong values for input argument #%d: Must be a divisor of %d.\n"), fname, 3, lA);
866         goto ERR;
867     }
868     /******************** get and check fourth argument (incr) ****************************************/
869     sciErr = getVectorIntArg(_pvCtx, 4, fname, &nincr, &Incr);
870     if (sciErr.iErr)
871     {
872         Scierror(sciErr.iErr, getErrorMessage(sciErr));
873         goto ERR;
874     }
875     if (nincr != ndims)
876     {
877         Scierror(999, _("%s: Incompatible input arguments #%d and #%d: Same sizes expected.\n"), fname, 3, 4);
878         goto ERR;
879     }
880
881     /* check values of Incr[i] */
882     if (Incr[0] <= 0)
883     {
884         Scierror(999, _("%s: Wrong values for input argument #%d: Positive integers expected.\n"), fname, 4);
885         goto ERR;
886     }
887     for (i = 0; i < ndims; i++)
888     {
889         if (lA % Incr[i])
890         {
891             Scierror(999, _("%s: Wrong values for input argument #%d: Elements must be divisors of %d.\n"), fname, 3, lA);
892             goto ERR;
893         }
894         if (i > 0 && (Incr[i] <= Incr[i - 1]))
895         {
896             Scierror(999, _("%s: Wrong values for input argument #%d: Elements must be in increasing ""order.\n"), fname, 4);
897             goto ERR;
898         }
899     }
900     if ((Dim = (int *)MALLOC((2 * ndims + 1) * sizeof(int))) == NULL)
901     {
902         Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
903         goto ERR;
904     }
905     if ((Sel = (int *)MALLOC((ndims) * sizeof(int))) == NULL)
906     {
907         Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
908         goto ERR;
909     }
910
911
912     /*Transform  Dim1 and Incr into Dim and Sel and check validity*/
913
914     nd = 0;
915     pd = 1;
916     if (Incr[0] != 1)
917     {
918         Dim[nd++] = Incr[0];
919         pd *= Incr[0];
920     }
921     Dim[nd++] = Dim1[0];
922     pd *= Dim1[0];
923     Sel[0] = nd;
924
925     for (k = 1; k < ndims; k++)
926     {
927         if (Incr[k] % pd != 0)
928         {
929             Scierror(999, _("%s: Incompatible input arguments #%d and #%d.\n"), fname, 3, 4);
930             goto ERR;
931         }
932         if (Incr[k] != pd)
933         {
934             Dim[nd++] = (int)(Incr[k] / pd);
935             pd = Incr[k];
936         }
937         Dim[nd++] = Dim1[k];
938         pd *= Dim1[k];
939         Sel[k] = nd;
940     }
941     if (pd < lA)
942     {
943         if (lA % pd != 0)
944         {
945             Scierror(999, _("%s: Incompatible input arguments #%d and #%d.\n"), fname, 3, 4);
946             goto ERR;
947         }
948         Dim[nd++] = (int)(lA / pd);
949     }
950
951     rank = ndims;
952     ndims = nd;
953
954     /* now one  same algorithm than sci_fft_3args applies */
955     /* Create  gdim struct */
956     gdim.rank = rank;
957     if ((gdim.dims = (fftw_iodim *)MALLOC(sizeof(fftw_iodim) * gdim.rank)) == NULL)
958     {
959         Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
960         goto ERR;
961     }
962
963     pd = 1; /* used to store prod(Dims(1:sel(k-1)))*/
964     pds = 1; /* used to store prod(Dims(sel(k-1):sel(k)))*/
965     j = 0;
966     for (i = 0; i < ndims; i++)
967     {
968         if (j >= gdim.rank) break;
969         if (Sel[j] == i + 1)
970         {
971             gdim.dims[j].n = Dim[i];
972             gdim.dims[j].is = pd;
973             gdim.dims[j].os = pd;
974             j++;
975         }
976         pd *= Dim[i];
977     }
978     /* Compute howmany_rank based on jumps in the Sel sequence */
979     gdim.howmany_rank = 0;
980     if ((Sel[0] != 1) && (Sel[0] != ndims)) gdim.howmany_rank++;
981
982     for (i = 1; i <= rank - 1; i++)
983     {
984         if (Sel[i] != Sel[i - 1] + 1) gdim.howmany_rank++;
985     }
986     if ((Sel[rank - 1] != ndims) || (rank == 1)) gdim.howmany_rank++;
987     /* Fill the howmany_dims struct */
988     if (gdim.howmany_rank > 0)
989     {
990         /* it must be the case */
991         int ih = 0;
992
993         if ((gdim.howmany_dims = (fftw_iodim *)MALLOC(gdim.howmany_rank * sizeof(fftw_iodim))) == NULL)
994         {
995             Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
996             goto ERR;
997         }
998         pd = 1;
999         for (j = 1; j <= (Sel[0] - 1); j++) pd *= Dim[j - 1]; /*prod(Dims(1:(sel(1)-1)))*/
1000         ih = 0;
1001         if ((Sel[0] != 1) && (Sel[0] != ndims))
1002         {
1003             /* First seleted dimension */
1004             gdim.howmany_dims[ih].is = 1;
1005             gdim.howmany_dims[ih].os = 1;
1006             gdim.howmany_dims[ih].n = pd;
1007             ih++;
1008         }
1009         pd *= Dim[Sel[0] - 1]; /*prod(Dims(1:sel(1)))*/
1010         for (i = 2; i <= rank; i++)
1011         {
1012             /* intermediate selected dimensions */
1013             if (Sel[i - 1] != Sel[i - 2] + 1)
1014             {
1015                 pds = 1;
1016                 for (j = (Sel[i - 2] + 1); j <= (Sel[i - 1] - 1); j++) pds *= Dim[j - 1]; /*prod(Dims(sel(i-1)+1:(sel(i)-1)))*/
1017                 gdim.howmany_dims[ih].is = pd;
1018                 gdim.howmany_dims[ih].os = pd;
1019                 gdim.howmany_dims[ih].n = pds;
1020                 ih++;
1021             }
1022             pd *= pds * Dim[Sel[i - 1] - 1]; /*prod(Dims(1:sel(i)))*/
1023         }
1024
1025         if (Sel[rank - 1] != ndims)
1026         {
1027             /* last selected dimension*/
1028             pds = 1;
1029             for (j = (Sel[rank - 1] + 1); j <= ndims; j++) pds *= Dim[j - 1]; /*prod(Dims(sel(i-1)+1:(sel(i)-1)))*/
1030             gdim.howmany_dims[ih].is = pd;
1031             gdim.howmany_dims[ih].os = pd;
1032             gdim.howmany_dims[ih].n = pds;
1033             ih++;
1034         }
1035         else if (rank == 1) /* the only selected dimension is the last one */
1036         {
1037             gdim.howmany_dims[ih].is = 1;
1038             gdim.howmany_dims[ih].os = 1;
1039             gdim.howmany_dims[ih].n = pd / Dim[Sel[0] - 1];
1040             ih++;
1041         }
1042     }
1043     if (!sci_fft_gen(_pvCtx, fname, ndimsA, dimsA, Ar,  Ai, isn, iopt, gdim))  goto ERR;
1044
1045     /***********************************
1046     * Return results in lhs argument *
1047     ***********************************/
1048
1049     ReturnArguments(_pvCtx);
1050 ERR:
1051     FREE(Dim1);
1052     FREE(Incr);
1053     FREE(Dim);
1054     FREE(Sel);
1055     FREE(gdim.dims);
1056     FREE(gdim.howmany_dims);
1057     return(0);
1058 }
1059 /*--------------------------------------------------------------------------*/
1060 BOOL isHyperMatrixMlist(void* _pvCtx, int *piAddressVar)
1061 {
1062     char **fields = NULL;
1063     SciErr sciErr;
1064     int iType = 0;
1065     int m = 0, n = 0;
1066
1067     if (piAddressVar == NULL)
1068     {
1069         return FALSE;
1070     }
1071
1072     sciErr = getVarType(_pvCtx, piAddressVar, &iType);
1073     if (sciErr.iErr)
1074     {
1075         return FALSE;
1076     }
1077
1078     if (iType == sci_mlist)
1079     {
1080         int* piAddrChild  = NULL;
1081         int iItem   = 0;
1082
1083         sciErr = getListItemNumber(_pvCtx, piAddressVar, &iItem);
1084         if (sciErr.iErr)
1085         {
1086             return FALSE;
1087         }
1088
1089         sciErr = getListItemAddress(_pvCtx, piAddressVar, 1, &piAddrChild);
1090         if (sciErr.iErr)
1091         {
1092             return FALSE;
1093         }
1094
1095         if (!isStringType(_pvCtx, piAddrChild))
1096         {
1097             return FALSE;
1098         }
1099
1100         if (getAllocatedMatrixOfString(_pvCtx, piAddrChild, &m, &n , &fields) == 0)
1101         {
1102             if (strcmp(fields[0], "hm") != 0)
1103             {
1104                 freeAllocatedMatrixOfString(m, n, fields);
1105                 fields = NULL;
1106                 return FALSE;
1107             }
1108             freeAllocatedMatrixOfString(m, n, fields);
1109             fields = NULL;
1110         }
1111         else
1112         {
1113             return FALSE;
1114         }
1115         return TRUE;
1116     }
1117     return FALSE;
1118 }
1119 /*--------------------------------------------------------------------------*/
1120 int sci_fft_gen(void* _pvCtx, char *fname, int ndimsA, int *dimsA, double *Ar,  double *Ai, int isn, int iopt, guru_dim_struct gdim)
1121 {
1122     /* API variables */
1123     SciErr sciErr;
1124
1125     /* Input  array variables */
1126     int  isrealA = (Ai == NULL), issymA = 1, lA = 1;
1127     /*for MKL*/
1128     int isrealA_save = isrealA ;
1129
1130     /*FFTW specific library variable */
1131     enum Scaling scale = None;
1132     enum Plan_Type type;
1133     fftw_plan p;
1134
1135     /* input/output address for transform variables */
1136     double *ri = NULL, *ii = NULL, *ro = NULL, *io = NULL;
1137
1138
1139     /* local variable */
1140     int one = 1;
1141     int i = 0;
1142
1143
1144     for (i = 0; i < ndimsA; i++)
1145     {
1146         lA *= dimsA[i];
1147     }
1148
1149
1150     if (iopt == 0)
1151     {
1152         /* automatically selected algorithm*/
1153         issymA =  check_array_symmetry(Ar, Ai, gdim);
1154         if (issymA < 0 )
1155         {
1156             Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
1157             goto ERR;
1158         }
1159     }
1160
1161     else if (iopt == 1)
1162     {
1163         issymA = 1; /* user forces symmetry */
1164     }
1165     else
1166     {
1167         issymA = 0;
1168     }
1169
1170     AssignOutputVariable(_pvCtx, 1) =  1;/* assume inplace transform*/
1171
1172     if (WITHMKL)
1173     {
1174         double dzero = 0.0;
1175         if (isrealA)
1176         {
1177             /*MKL does not implement the r2c nor r2r guru split methods, make A complex */
1178             if (issymA)
1179             {
1180                 /* result will be real, the imaginary part of A can be allocated alone */
1181                 sciErr = allocMatrixOfDouble(_pvCtx, *getNbInputArgument(_pvCtx) + 1, 1, lA, &Ai);
1182                 if (sciErr.iErr)
1183                 {
1184                     Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
1185                     goto ERR;
1186                 }
1187                 C2F(dset)(&lA, &dzero, Ai, &one);
1188             }
1189             else
1190             {
1191                 /* result will be complex, realloc A for inplace computation */
1192                 sciErr = allocComplexArrayOfDouble(_pvCtx, *getNbInputArgument(_pvCtx) + 1, ndimsA, dimsA, &ri, &Ai);
1193                 if (sciErr.iErr)
1194                 {
1195                     Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
1196                     goto ERR;
1197                 }
1198                 C2F(dcopy)(&lA, Ar, &one, ri, &one);
1199                 Ar = ri;
1200                 C2F(dset)(&lA, &dzero, Ai, &one);
1201                 AssignOutputVariable(_pvCtx, 1) =  nbInputArgument(_pvCtx) + 1;
1202                 isrealA = 0;
1203             }
1204         }
1205     }
1206
1207     if (!isrealA && issymA) /* A is complex but result is real */
1208     {
1209         /* result will be complex, realloc real part of A for real part inplace computation */
1210         sciErr = allocArrayOfDouble(_pvCtx, *getNbInputArgument(_pvCtx) + 1, ndimsA, dimsA, &ri);
1211         if (sciErr.iErr)
1212         {
1213             Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
1214             goto ERR;
1215         }
1216         C2F(dcopy)(&lA, Ar, &one, ri, &one);
1217         Ar = ri;
1218         AssignOutputVariable(_pvCtx, 1) = nbInputArgument(_pvCtx) + 1;
1219     }
1220
1221     /* Set pointers on real and imaginary part of the input */
1222     ri = Ar;
1223     ii = Ai;
1224
1225     scale = None; /*no scaling needed */
1226     if (isn == FFTW_BACKWARD) scale = Divide;
1227     if (isrealA)
1228     {
1229         /*A is real */
1230         if (issymA)
1231         {
1232             /*r2r =  isrealA &&  issymA*/
1233             /* there is no general plan able to compute r2r transform so it is tranformed into
1234             a R2c plan. The computed imaginary part will be zero*/
1235             sciErr = allocMatrixOfDouble(_pvCtx, *getNbInputArgument(_pvCtx) + 1, 1, lA,  &io);
1236             if (sciErr.iErr)
1237             {
1238                 Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
1239                 goto ERR;
1240             }
1241             type = R2C_PLAN;
1242             ro = Ar;
1243         }
1244         else
1245         {
1246             /*r2c =  isrealA && ~issymA;*/
1247             /* transform cannot be done in place */
1248             sciErr = allocComplexArrayOfDouble(_pvCtx, *getNbInputArgument(_pvCtx) + 1, ndimsA, dimsA, &ro, &io);
1249             if (sciErr.iErr)
1250             {
1251                 Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
1252                 goto ERR;
1253             }
1254             AssignOutputVariable(_pvCtx, 1) = nbInputArgument(_pvCtx) + 1;
1255             type = R2C_PLAN; /* fftw_plan_guru_split_dft_r2c plans for an FFTW_FORWARD transform*/
1256             if (isn == FFTW_BACKWARD)
1257             {
1258                 /*transform problem into a FORWARD fft*/
1259                 /*ifft(A)=conj(fft(A/N)) cas vect*/
1260                 /* pre traitement A must be  divided by N cas vect*/
1261                 /* post treatment result must conjugated */
1262             }
1263         }
1264     }
1265     else
1266     {
1267         /* A is complex */
1268         if (!WITHMKL && issymA) /*result is real*/
1269         {
1270             /*c2r =  ~isrealA &&  issymA*/
1271             ro = ri;
1272             io = NULL;
1273
1274             type = C2R_PLAN; /*fftw_plan_guru_split_dft_c2r plans for an FFTW_BACKWARD transform*/
1275             if (isn == FFTW_FORWARD)
1276             {
1277                 /*transform problem into a BACKWARD fft : fft(A)=ifft(conj(A))*/
1278                 double minusone = -1.0;
1279                 C2F(dscal)(&lA, &minusone, ii, &one);
1280             }
1281         }
1282         else
1283         {
1284             /*c2c =  ~isrealA && ~issymA;*/
1285             /* use inplace transform*/
1286             type = C2C_PLAN; /*  fftw_plan_guru_split_dft plans for an FFTW_FORWARD transform*/
1287             if (isn == FFTW_BACKWARD)
1288             {
1289                 /*transform problem into a FORWARD fft*/
1290                 /* ifft(A) = %i*conj(fft(%i*conj(A)/N) */
1291                 /* reverse input */
1292                 ri = Ai;
1293                 ii = Ar;
1294                 /* reverse output */
1295                 ro = Ai;
1296                 io = Ar;
1297             }
1298             else
1299             {
1300                 ro = ri;
1301                 io = ii;
1302             }
1303         }
1304     }
1305     /* Set Plan */
1306     p = GetFFTWPlan(type, &gdim, ri, ii, ro, io, getCurrentFftwFlags(), isn);
1307     if (p == NULL)
1308     {
1309         Scierror(999, _("%s: No more memory.\n"), fname);
1310         goto ERR;
1311     }
1312     /* pre-treatment */
1313     if (scale != None)
1314     {
1315         double ak = 1.0;
1316         for (i = 0; i < gdim.rank; i++) ak = ak * ((double)(gdim.dims[i].n));
1317         if (scale == Divide) ak = 1.0 / ak;
1318         C2F(dscal)(&lA, &ak, ri, &one);
1319         if (isrealA == 0) C2F(dscal)(&lA, &ak, ii, &one);
1320     }
1321     /* execute FFTW plan */
1322     ExecuteFFTWPlan(type, p, ri, ii, ro, io);
1323     /* Post treatment */
1324     switch (type)
1325     {
1326     case R2R_PLAN:
1327         if (complete_array(ro, NULL, gdim) == -1)
1328         {
1329             Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
1330             goto ERR;
1331         }
1332         break;
1333     case C2R_PLAN:
1334         break;
1335     case R2C_PLAN:
1336         if (issymA)
1337         {
1338             /*R2C has been used to solve an r2r problem*/
1339             if (complete_array(ro, NULL, gdim) == -1)
1340             {
1341                 Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
1342                 goto ERR;
1343             }
1344         }
1345         else
1346         {
1347             if (complete_array(ro, io, gdim) == -1)
1348             {
1349                 Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
1350                 goto ERR;
1351             }
1352             if (isn == FFTW_BACKWARD)
1353             {
1354                 /*conjugate result */
1355                 double ak = -1.0;
1356                 C2F(dscal)(&lA, &ak, io, &one);
1357             }
1358         }
1359         break;
1360     case C2C_PLAN:
1361         if (WITHMKL && isrealA_save)
1362         {
1363             if (isn == FFTW_FORWARD)
1364             {
1365                 if (complete_array(ro, io, gdim) == -1)
1366                 {
1367                     Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
1368                     goto ERR;
1369                 }
1370             }
1371             else
1372             {
1373                 if (complete_array(io, ro, gdim) == -1)
1374                 {
1375                     Scierror(999, _("%s: Cannot allocate more memory.\n"), fname);
1376                     goto ERR;
1377                 }
1378             }
1379         }
1380         break;
1381     }
1382
1383     return(1);
1384 ERR:
1385     return(0);
1386 }