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