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 }