scilab/modules/fftw/sci_gateway/cpp/sci_fftw.cpp

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