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 - INRIA - Serge STEER
6  *
7  * Copyright (C) 2012 - 2016 - Scilab Enterprises
8  *
9  * This file is hereby licensed under the terms of the GNU GPL v2.0,
10  * pursuant to article 5.3.4 of the CeCILL v.2.1.
11  * This file was originally licensed under the terms of the CeCILL v2.1,
12  * and continues to be available under such terms.
13  * For more information, see the COPYING file which you should have received
14  * along with this program.
15  *
16  */
17 #include <math.h>
18 #include "fftw_utilities.h"
19 #include "sci_malloc.h"
20 #include "callfftw.h"
21 int check_1D_symmetry(double *Ar, double *Ai, int nA, int iA);
22 int check_2D_symmetry(double *Ar, double *Ai, int mA, int iA, int nA, int jA);
23 int check_ND_symmetry(double *Ar, double *Ai, int ndims, int *dims, int *incr);
24 
25 void complete_1D_array(double *Ar, double *Ai, int nA, int iA);
26 void complete_2D_array(double *Ar, double *Ai, int mA, int iA, int nA, int jA);
27 int complete_ND_array(double *Ar, double *Ai, int ndims, int *dims, int *incr);
28 void dct_scale_1D_array(double *Ar, double *Ai, int nA, int iA, int isn, double fact);
29 void dct_scale_2D_array(double *Ar, double *Ai, int mA, int iA, int nA, int jA, int isn, double fact);
30 int dct_scale_ND_array(double *Ar, double *Ai, int ndims, int *dims, int *incr, int isn, double fact);
31 
32 void dst_scale_1D_array(double *Ar, double *Ai, int nA, int iA, int isn, double fact);
33 void dst_scale_2D_array(double *Ar, double *Ai, int mA, int iA, int nA, int jA, int isn, double fact);
34 int dst_scale_ND_array(double *Ar, double *Ai, int ndims, int *dims, int *incr, int isn, double fact);
35 
36 /*--------------------------------------------------------------------------*/
37 /* definition of structures to store parameters
38  * of FFTW planners - set here default value -
39  */
40 FFTW_Plan_struct Sci_Forward_Plan =
41 {
42     0,                  /* int plan_type            */
43     NULL,               /* fftw_plan p              */
44     {0, NULL, 0, NULL}, /* guru_dim_struct gdim     */
45     FFTW_ESTIMATE,      /* unsigned flags           */
46     NULL                /* kind                     */
47 };
48 
49 FFTW_Plan_struct Sci_Backward_Plan =
50 {
51     0,                  /* int plan_type            */
52     NULL,               /* fftw_plan p              */
53     {0, NULL, 0, NULL}, /* guru_dim_struct gdim     */
54     FFTW_ESTIMATE,      /* unsigned flags           */
55     NULL                /* kind                     */
56 };
57 /*--------------------------------------------------------------------------*/
getSci_Backward_Plan(void)58 FFTW_Plan_struct *getSci_Backward_Plan(void)
59 {
60     return &Sci_Backward_Plan;
61 }
62 /*--------------------------------------------------------------------------*/
getSci_Forward_Plan(void)63 FFTW_Plan_struct *getSci_Forward_Plan(void)
64 {
65     return &Sci_Forward_Plan;
66 }
67 /*--------------------------------------------------------------------------*/
68 unsigned cur_fftw_flags = FFTW_ESTIMATE;
69 /*--------------------------------------------------------------------------*/
getCurrentFftwFlags(void)70 unsigned int getCurrentFftwFlags(void)
71 {
72     return cur_fftw_flags;
73 }
74 /*--------------------------------------------------------------------------*/
setCurrentFftwFlags(unsigned int newFftwFlags)75 void setCurrentFftwFlags(unsigned int newFftwFlags)
76 {
77     cur_fftw_flags = newFftwFlags;
78 }
79 /*--------------------------------------------------------------------------*/
80 /* Free a FFTW_Plan_struct
81  *
82  * Input : FFTW_Plan_struct *Sci_Plan
83  *
84  * Output : int, return 1.
85  *
86  */
FreeFFTWPlan(FFTW_Plan_struct * Sci_Plan)87 int FreeFFTWPlan(FFTW_Plan_struct *Sci_Plan)
88 {
89     if (Sci_Plan->p != NULL)
90     {
91         call_fftw_destroy_plan(Sci_Plan->p);
92         Sci_Plan->p = NULL;
93     }
94     if (Sci_Plan->gdim.rank != 0)
95     {
96         Sci_Plan->gdim.rank = 0;
97         FREE(Sci_Plan->gdim.dims);
98         Sci_Plan->gdim.dims = NULL;
99         FREE(Sci_Plan->kind);
100         Sci_Plan->kind = NULL;
101     }
102     if (Sci_Plan->gdim.howmany_rank != 0)
103     {
104         Sci_Plan->gdim.howmany_rank = 0;
105         FREE(Sci_Plan->gdim.howmany_dims);
106         Sci_Plan->gdim.howmany_dims = NULL;
107     }
108 
109     return (1);
110 }
111 /*--------------------------------------------------------------------------*/
112 /* Return a valid plan ptr.
113  * This function search in the Sci_xx_Plan structures if
114  * the given input parameters follows an already stored
115  * set of parameters.
116  * If found then return an already stored plan ptr.
117  * If not found, then returns
118  * a new set of parameters (and a new plan)
119  * with fftw_plan_guru_split_dft
120  * and store it in Sci_xx_Plan structures
121  *
122  * Input : guru_dim_struct *gdim
123  *         double *ri, double *ii
124  *         double *ro, double *io
125  *         unsigned flags, int isn
126  *
127  * Output : fftw_plan
128  *
129  *
130  */
GetFFTWPlan(enum Plan_Type type,guru_dim_struct * gdim,double * ri,double * ii,double * ro,double * io,unsigned flags,int isn,fftw_r2r_kind * kind,int * errflag)131 fftw_plan GetFFTWPlan(enum Plan_Type type, guru_dim_struct *gdim,
132                       double *ri, double *ii,
133                       double *ro, double *io,
134                       unsigned flags, int isn, fftw_r2r_kind *kind, int *errflag)
135 {
136     FFTW_Plan_struct *Sci_Plan;
137     int i = 0;
138 
139     *errflag = 0;
140 
141     if (isn == -1)
142     {
143         Sci_Plan = &Sci_Backward_Plan;
144     }
145     else
146     {
147         Sci_Plan = &Sci_Forward_Plan;
148     }
149 
150     /* plan must be changed */
151     FreeFFTWPlan(Sci_Plan);
152 
153     Sci_Plan->plan_type = type;
154     if (gdim->rank != 0)
155     {
156         Sci_Plan->gdim.rank = gdim->rank;
157         if ((Sci_Plan->gdim.dims = (fftw_iodim *) MALLOC(sizeof(fftw_iodim) * (gdim->rank))) == NULL)
158         {
159             *errflag = 1;
160             return (NULL);
161         }
162         for (i = 0; i < gdim->rank; i++)
163         {
164             Sci_Plan->gdim.dims[i].n  = gdim->dims[i].n;
165             Sci_Plan->gdim.dims[i].is = gdim->dims[i].is;
166             Sci_Plan->gdim.dims[i].os = gdim->dims[i].os;
167         }
168 
169         if (kind != NULL)
170         {
171             if ((Sci_Plan->kind = (fftw_r2r_kind *) MALLOC(sizeof(fftw_r2r_kind) * (gdim->rank))) == NULL)
172             {
173                 FREE(Sci_Plan->gdim.dims);
174                 *errflag = 1;
175                 return (NULL);
176             }
177             for (i = 0; i < gdim->rank; i++)
178             {
179                 Sci_Plan->kind[i]  = kind[i];
180             }
181         }
182     }
183 
184     if (gdim->howmany_rank != 0)
185     {
186         Sci_Plan->gdim.howmany_rank = gdim->howmany_rank;
187         if ((Sci_Plan->gdim.howmany_dims = (fftw_iodim *) MALLOC(sizeof(fftw_iodim) * (gdim->howmany_rank))) == NULL)
188         {
189             FREE(Sci_Plan->gdim.dims);
190             *errflag = 1;
191             return (NULL);
192         }
193         for (i = 0; i < gdim->howmany_rank; i++)
194         {
195             Sci_Plan->gdim.howmany_dims[i].n  = gdim->howmany_dims[i].n;
196             Sci_Plan->gdim.howmany_dims[i].is = gdim->howmany_dims[i].is;
197             Sci_Plan->gdim.howmany_dims[i].os = gdim->howmany_dims[i].os;
198         }
199     }
200 
201     Sci_Plan->flags = cur_fftw_flags;
202 
203     switch (type)
204     {
205         case C2C_PLAN:
206             Sci_Plan->p = call_fftw_plan_guru_split_dft(Sci_Plan->gdim.rank,
207                           Sci_Plan->gdim.dims,
208                           Sci_Plan->gdim.howmany_rank,
209                           Sci_Plan->gdim.howmany_dims,
210                           ri, ii, ro, io,
211                           Sci_Plan->flags);
212             break;
213         case C2R_PLAN:
214             Sci_Plan->p = call_fftw_plan_guru_split_dft_c2r(Sci_Plan->gdim.rank,
215                           Sci_Plan->gdim.dims,
216                           Sci_Plan->gdim.howmany_rank,
217                           Sci_Plan->gdim.howmany_dims,
218                           ri, ii, ro, flags);
219             break;
220         case R2C_PLAN:
221             Sci_Plan->p = call_fftw_plan_guru_split_dft_r2c(Sci_Plan->gdim.rank,
222                           Sci_Plan->gdim.dims,
223                           Sci_Plan->gdim.howmany_rank,
224                           Sci_Plan->gdim.howmany_dims,
225                           ri, ro, io, flags);
226 
227             break;
228         case R2R_PLAN:
229             Sci_Plan->p = call_fftw_plan_guru_split_dft_r2r(Sci_Plan->gdim.rank,
230                           Sci_Plan->gdim.dims,
231                           Sci_Plan->gdim.howmany_rank,
232                           Sci_Plan->gdim.howmany_dims,
233                           ri, ro, kind, flags);
234             break;
235     }
236 
237     if (Sci_Plan->p == NULL)
238     {
239         *errflag = 2;
240     }
241 
242     return (Sci_Plan->p);
243 }
244 /*--------------------------------------------------------------------------*/
245 /* Check if two guru_dim structures are equal
246  *
247  * Input : guru_dim_struct *gdim1
248  *         guru_dim_struct *gdim2
249  *
250  * Output : int, return 0 if False, else 1.
251  *
252  */
CheckGuruDims(guru_dim_struct * gdim1,guru_dim_struct * gdim2)253 int CheckGuruDims(guru_dim_struct *gdim1, guru_dim_struct *gdim2)
254 {
255     int i;
256 
257     if ( (gdim1->rank == gdim2->rank) &&
258             (gdim1->howmany_rank == gdim2->howmany_rank))
259     {
260         for (i = 0; i < gdim1->rank; i++)
261         {
262             if (gdim1->dims[i].n  != gdim2->dims[i].n)
263             {
264                 return (0);
265             }
266             if (gdim1->dims[i].is != gdim2->dims[i].is)
267             {
268                 return (0);
269             }
270             if (gdim1->dims[i].os != gdim2->dims[i].os)
271             {
272                 return (0);
273             }
274         }
275         for (i = 0; i < gdim1->howmany_rank; i++)
276         {
277             if (gdim1->howmany_dims[i].n  != gdim2->howmany_dims[i].n)
278             {
279                 return (0);
280             }
281             if (gdim1->howmany_dims[i].is != gdim2->howmany_dims[i].is)
282             {
283                 return (0);
284             }
285             if (gdim1->howmany_dims[i].os != gdim2->howmany_dims[i].os)
286             {
287                 return (0);
288             }
289         }
290         return (1);
291     }
292     else
293     {
294         return (0);
295     }
296 }
297 /*--------------------------------------------------------------------------*/
298 /* Check if kind array  are equal
299  *
300  * Input : fftw_r2r_kind *kind1
301  *         gfftw_r2r_kind *kind2
302  *         int rank
303  *
304  * Output : int, return 0 if False, else 1.
305  *
306  */
CheckKindArray(fftw_r2r_kind * kind1,fftw_r2r_kind * kind2,int rank)307 int CheckKindArray(fftw_r2r_kind *kind1, fftw_r2r_kind *kind2, int rank)
308 {
309     int i;
310     if ((kind1 == NULL) && (kind2 == NULL))
311     {
312         return (1);
313     }
314     if ((kind1 == NULL) || (kind2 == NULL))
315     {
316         return (0);
317     }
318     for (i = 0; i < rank; i++)
319     {
320         if (kind1[i]  != kind2[i])
321         {
322             return (0);
323         }
324     }
325     return (1);
326 }
327 /*--------------------------------------------------------------------------*/
328 /* call different fftw_execute_split_dft_xxx procedures according to type input
329  *
330  * Input :  Plan_Type type
331  *          fftw_plan p
332  *          double *ri
333  *          double *ii
334  * Output : double *ro
335  *          double *io
336  */
ExecuteFFTWPlan(enum Plan_Type type,const fftw_plan p,double * ri,double * ii,double * ro,double * io)337 void ExecuteFFTWPlan(enum Plan_Type type, const fftw_plan p, double *ri, double *ii, double *ro, double *io)
338 {
339     switch (type)
340     {
341         case C2C_PLAN:
342             call_fftw_execute_split_dft(p, ri, ii, ro, io);
343             break;
344         case C2R_PLAN:
345             call_fftw_execute_split_dft_c2r(p, ri, ii, ro);
346             break;
347         case R2C_PLAN:
348             call_fftw_execute_split_dft_r2c(p, ri, ro, io);
349             break;
350         case R2R_PLAN:
351             call_fftw_execute_split_dft_r2r(p, ri, ro);
352             break;
353     }
354 }
355 /*--------------------------------------------------------------------------*/
is_real(double * Ar,double * Ai,int ndims,int * dims)356 int is_real(double *Ar, double *Ai, int ndims, int *dims)
357 {
358     double zero = 0.0;
359     int t = 1;
360     int i = 0;
361     int lA = 1;
362 
363 
364     for (i = 0; i < ndims; i++)
365     {
366         lA = lA * dims[i];
367     }
368 
369     /*Check if A is real*/
370     if (Ai != NULL)
371     {
372         for (i = 0; i < lA; i++)
373         {
374             if (Ai[i] != zero)
375             {
376                 t = 0;
377                 break;
378             }
379         }
380     }
381     return t;
382 }
383 
384 /*--------------------------------------------------------------------------
385  * Check if a 1D array A  is "symmetric" or hermitian symmetric for fft.
386  * A==conj(A([1 $:-1:2]))
387  * Ar : pointer on real part array
388  * Ai : pointer on imaginary part array or NULL
389  * nA : number of elements
390  * iA : increment between 2 consecutive element of the array
391  */
392 
check_1D_symmetry(double * Ar,double * Ai,int nA,int iA)393 int check_1D_symmetry(double *Ar, double *Ai, int nA, int iA)
394 {
395     int i = 0;
396     int nas2 = (int)(nA / 2);
397     double zero = 0.0;
398     //Checking real part
399     if ( nA % 2 == 0)
400     {
401         /* A length is even */
402         for (i = 1; i < nas2; i++)
403         {
404             if (Ar[iA * i] != Ar[iA * (nA - i)])
405             {
406                 return 0;
407             }
408         }
409     }
410     else
411     {
412         /* A length is odd */
413         for (i = 1; i <= nas2; i++)
414         {
415             if (Ar[iA * i] != Ar[iA * (nA - i)])
416             {
417                 return 0;
418             }
419         }
420     }
421     if (Ai == NULL)
422     {
423         return 1;
424     }
425     //Checking imaginary part
426     if ( nA % 2 == 0)
427     {
428         /* A length is even */
429         if (Ai[0] != zero || Ai[iA * (nA / 2)] != zero)
430         {
431             return 0;
432         }
433         for (i = 1; i < nas2; i++)
434         {
435             if (Ai[iA * i] != -Ai[iA * (nA - i)])
436             {
437                 return 0;
438             }
439         }
440     }
441     else
442     {
443         /* A length is odd */
444         if (Ai[0] != zero)
445         {
446             return 0;
447         }
448         for (i = 1; i <= nas2; i++)
449         {
450             if (Ai[iA * i] != -Ai[iA * (nA - i)])
451             {
452                 return 0;
453             }
454         }
455     }
456     return 1;
457 }
458 /*--------------------------------------------------------------------------
459  * Check if a 2D array A  is "symmetric" or hermitian symmetric for fft.
460  * A==conj(A([1 $:-1:2],[1 $:-1:2])
461  * Ar : pointer on real part array
462  * Ai : pointer on imaginary part array or NULL
463  * mA : number of rows
464  * nA : number of columns
465  * iA : increment between 2 consecutive element of a row
466  * jA : increment between 2 consecutive element of a column
467  */
468 
check_2D_symmetry(double * Ar,double * Ai,int mA,int iA,int nA,int jA)469 int check_2D_symmetry(double *Ar, double *Ai, int mA, int iA, int nA, int jA)
470 {
471     int l1 = 0;
472     int l2 = 0;
473     int k = 0;
474     int l = 0;
475     int nAs2 = (int)(nA / 2) + 1; /* A VERIFIER */
476 
477     /* Check first column */
478     if (!check_1D_symmetry(Ar, Ai, mA, iA))
479     {
480         return 0;
481     }
482     /* Check first row */
483     if (!check_1D_symmetry(Ar, Ai, nA, jA))
484     {
485         return 0;
486     }
487 
488     /* Check A(2:$,2:$) block */
489     if (Ai == NULL)
490     {
491         for (k = 1; k < nAs2; k++)
492         {
493             l1 = jA * k + iA ;
494             l2 = jA * (nA - k) + iA * (mA - 1);
495             for (l = 1; l < mA; l++)
496             {
497                 if (Ar[l1] != Ar[l2])
498                 {
499                     return 0;
500                 }
501                 l1 += iA;
502                 l2 -= iA;
503             }
504         }
505     }
506     else
507     {
508         for (k = 1; k < nAs2; k++)
509         {
510             l1 = jA * k + iA ;
511             l2 = jA * (nA - k) + iA * (mA - 1);
512             for (l = 1; l < mA; l++)
513             {
514                 if ((Ar[l1] != Ar[l2]) || (Ai[l1] != -Ai[l2]))
515                 {
516                     return 0;
517                 }
518                 l1 += iA;
519                 l2 -= iA;
520             }
521         }
522     }
523     return 1;
524 }
525 
526 /*--------------------------------------------------------------------------
527  * Check if a N-D array A  is "symmetric" or hermitian symmetric for fft
528  * A==conj(A([1 $:-1:2],...,[1 $:-1:2])
529  * Ar : pointer on real part array
530  * Ai : pointer on imaginary part array or NULL
531  * mA : number of rows
532  * nA : number of columns
533  * iA : increment between 2 consecutive element of a row
534  * jA : increment between 2 consecutive element of a column
535  */
536 
check_ND_symmetry(double * Ar,double * Ai,int ndims,int * dims,int * incr)537 int check_ND_symmetry(double *Ar, double *Ai, int ndims, int *dims, int *incr)
538 {
539     int i = 0, j = 0, l = 0;
540     int r = 0;
541     int l1 = 0;/* current 1-D index in array*/
542     int l2 = 0;/* associated symmetric value 1-D index */
543     int *temp = NULL;
544     int *dims1 = NULL;
545     int *incr1 = NULL;
546     int nSub = 0, nSubs2 = 0;
547     int k = 0, step = 0;
548 
549     if (ndims == 2)
550     {
551         r = check_2D_symmetry(Ar, Ai, dims[0], incr[0], dims[1], incr[1]);
552         return r;
553     }
554     else if  (ndims == 1)
555     {
556         r = check_1D_symmetry(Ar, Ai, dims[0], incr[0]);
557         return r;
558     }
559 
560     if  ((temp =  (int *)MALLOC(sizeof(int) * (2 * ndims))) == NULL)
561     {
562         return -1;
563     }
564     dims1 = temp;
565     incr1 = temp + ndims;
566 
567     for (i = 0; i < ndims; i++)
568     {
569         /* remove current dimension and increment out of  dims ans incr */
570         l = 0;
571         for   (j = 0; j < ndims; j++)
572         {
573             if (j != i)
574             {
575                 dims1[l] = dims[j];
576                 incr1[l] = incr[j];
577                 l++;
578             }
579         }
580         r = check_ND_symmetry(Ar, Ai, ndims - 1, dims1, incr1);
581         if (r != 1)
582         {
583             dims1 = NULL;
584             incr1 = NULL;
585             FREE(temp);
586             return r;
587         }
588     }
589 
590     /* check bloc A(2:$,....,2:$)*/
591     /*A(2,...,2) index*/
592     l1 = 0;
593     for (i = 0; i < ndims; i++)
594     {
595         l1 += incr[i];
596     }
597     /*A($,...,$) index*/
598     l2 = 0;
599     for (i = 0; i < ndims; i++)
600     {
601         l2 += (dims[i] - 1) * incr[i];
602     }
603 
604 
605     /* cumprod(size(A(2:$,....,2:$)))*/
606     incr1[0] = dims[0] - 1;
607     for (i = 1; i < (ndims - 1); i++)
608     {
609         incr1[i] = incr1[i - 1] * (dims[i] - 1) ;
610     }
611     /* steps*/
612     dims1[0] = (dims[0] - 2) * incr[0];
613     for (i = 1; i < (ndims - 1); i++)
614     {
615         dims1[i] = dims1[i - 1] + (dims[i] - 2) * incr[i];
616     }
617 
618     /*  A(2:$,....,2:$) block number of elements*/
619     nSub = 1;
620     for (i = 0; i < ndims; i++)
621     {
622         nSub *= (dims[i] - 1);
623     }
624 
625     nSubs2 = (int)(nSub / 2);
626 
627 
628     if (Ai == NULL)
629     {
630         /* Real case */
631         for (i = 0; i < nSubs2; i++)
632         {
633 
634             if (Ar[l1] != Ar[l2])
635             {
636                 return 0;
637             }
638             step = incr[0];
639             for (j = ndims - 2; j >= 0; j--)
640             {
641                 if ((i + 1) % incr1[j] == 0)
642                 {
643                     step = -dims1[j] + incr[j + 1] ;
644                     break;
645                 }
646             }
647             l1 += step;
648             l2 -= step;
649         }
650     }
651     else    /* Complex case */
652     {
653         for (i = 0; i < nSubs2; i++)
654         {
655             if (Ar[l1] != Ar[l2] || Ai[l1] != -Ai[l2])
656             {
657                 return 0;
658             }
659             step = incr[0];
660             for (j = ndims - 2; j >= 0; j--)
661             {
662                 if ((i + 1) % incr1[j] == 0)
663                 {
664                     step = -dims1[j] + incr[j + 1] ;
665                     break;
666                 }
667             }
668             l1 += step;
669             l2 -= step;
670         }
671     }
672     dims1 = NULL;
673     incr1 = NULL;
674     FREE(temp);
675     return 1;
676 }
677 
678 
679 
check_array_symmetry(double * Ar,double * Ai,guru_dim_struct gdim)680 int check_array_symmetry(double *Ar, double *Ai, guru_dim_struct gdim)
681 {
682     int ndims = gdim.rank;
683     int * dims = NULL;
684     int * incr = NULL;
685     int r = -1;
686     int i = 0, j = 0, k = 0;
687 
688     if (gdim.howmany_rank == 0)
689     {
690         switch (gdim.rank)
691         {
692             case 1:
693                 return check_1D_symmetry(Ar, Ai, gdim.dims[0].n, gdim.dims[0].is);
694             case 2:
695                 return check_2D_symmetry(Ar, Ai, gdim.dims[0].n, gdim.dims[0].is, gdim.dims[1].n, gdim.dims[1].is);
696             default: /*general N-D case*/
697                 if ((dims = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
698                 {
699                     return -1;
700                 }
701                 if ((incr = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
702                 {
703                     FREE(dims);
704                     return -1;
705                 }
706                 for (i = 0; i < ndims; i++)
707                 {
708                     dims[i] = gdim.dims[i].n;
709                     incr[i] = gdim.dims[i].is;
710                 }
711                 r = check_ND_symmetry(Ar, Ai, ndims, dims, incr);
712                 FREE(dims);
713                 FREE(incr);
714                 return r;
715         }
716     }
717     else
718     {
719         int m = 0;
720         int p = 1;
721         int *dims1 = NULL;
722         int *incr1 = NULL;
723         int ir = 0;
724 
725         if ((dims1 = (int *)MALLOC(sizeof(int) * gdim.howmany_rank)) == NULL)
726         {
727             return -1;
728         }
729         dims1[0] = gdim.howmany_dims[0].n;
730         for (i = 1; i < gdim.howmany_rank; i++)
731         {
732             dims1[i] = dims1[i - 1] * gdim.howmany_dims[i].n;
733         }
734         m = dims1[gdim.howmany_rank - 1];
735 
736         if ((incr1 = (int *)MALLOC(sizeof(int) * gdim.howmany_rank)) == NULL)
737         {
738             FREE(dims1);
739             return -1;
740         }
741         p = 1;
742         for (i = 0; i < gdim.howmany_rank; i++)
743         {
744             p += (gdim.howmany_dims[i].n - 1) * gdim.howmany_dims[i].is;
745             incr1[i] = p;
746         }
747         switch (gdim.rank)
748         {
749             case 1:
750                 if (Ai == NULL)
751                 {
752                     /* multiple 1D fft */
753                     for (ir = 0; ir < gdim.howmany_rank; ir++)
754                     {
755                         j = 0;
756                         for (i = 1; i <= m; i++)
757                         {
758                             if ((r = check_1D_symmetry(Ar + j, NULL, gdim.dims[0].n, gdim.dims[0].is)) != 1 )
759                             {
760                                 FREE(dims1);
761                                 FREE(incr1);
762                                 return r;
763                             }
764                             j += gdim.howmany_dims[0].is;
765                             for (k = gdim.howmany_rank - 2; k >= 0; k--)
766                             {
767                                 if (i % dims1[k] == 0)
768                                 {
769                                     j += -incr1[k] + gdim.howmany_dims[k + 1].is;
770                                     break;
771                                 }
772                             }
773                         }
774                     }
775                 }
776                 else
777                 {
778                     for (ir = 0; ir < gdim.howmany_rank; ir++)
779                     {
780                         j = 0;
781                         for (i = 1; i <= m; i++)
782                         {
783                             if ((r = check_1D_symmetry(Ar + j, Ai + j, gdim.dims[0].n, gdim.dims[0].is)) != 1 )
784                             {
785                                 FREE(dims1);
786                                 FREE(incr1);
787                                 return r;
788                             }
789                             j += gdim.howmany_dims[0].is;
790                             for (k = gdim.howmany_rank - 2; k >= 0; k--)
791                             {
792                                 if (i % dims1[k] == 0)
793                                 {
794                                     j += -incr1[k] + gdim.howmany_dims[k + 1].is;
795                                     break;
796                                 }
797                             }
798                         }
799                     }
800                 }
801                 FREE(dims1);
802                 FREE(incr1);
803                 return 1;
804             case 2:  /* multiple 2D fft */
805                 if (Ai == NULL)
806                 {
807                     for (ir = 0; ir < gdim.howmany_rank; ir++)
808                     {
809                         j = 0;
810                         for (i = 1; i <= m; i++)
811                         {
812                             if ((r = check_2D_symmetry(Ar + j, NULL, gdim.dims[0].n, gdim.dims[0].is,
813                                                        gdim.dims[1].n, gdim.dims[1].is)) != 1 )
814                             {
815                                 FREE(dims1);
816                                 FREE(incr1);
817                                 return r;
818                             }
819                             j += gdim.howmany_dims[0].is;
820 
821                             for (k = gdim.howmany_rank - 2; k >= 0; k--)
822                             {
823                                 if (i % dims1[k] == 0)
824                                 {
825                                     j += -incr1[k] + gdim.howmany_dims[k + 1].is;
826                                     break;
827                                 }
828                             }
829                         }
830                     }
831                 }
832                 else
833                 {
834                     for (ir = 0; ir < gdim.howmany_rank; ir++)
835                     {
836                         j = 0;
837                         for (i = 1; i <= m; i++)
838                         {
839                             if ((r = check_2D_symmetry(Ar + j, Ai + j, gdim.dims[0].n, gdim.dims[0].is,
840                                                        gdim.dims[1].n, gdim.dims[1].is)) != 1 )
841                             {
842                                 FREE(dims1);
843                                 FREE(incr1);
844                                 return r;
845                             }
846                             j += gdim.howmany_dims[0].is;
847                             for (k = gdim.howmany_rank - 2; k >= 0; k--)
848                             {
849                                 if (i % dims1[k] == 0)
850                                 {
851                                     j += -incr1[k] + gdim.howmany_dims[k + 1].is;
852                                     break;
853                                 }
854                             }
855                         }
856                     }
857                 }
858                 FREE(dims1);
859                 FREE(incr1);
860                 return 1;
861             default: /*general N-D case*/
862                 if ((dims = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
863                 {
864                     FREE(dims1);
865                     FREE(incr1);
866                     return -1;
867                 }
868                 if ((incr = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
869                 {
870                     FREE(dims);
871                     FREE(dims1);
872                     FREE(incr1);
873                     return -1;
874                 }
875                 for (i = 0; i < ndims; i++)
876                 {
877                     dims[i] = gdim.dims[i].n;
878                     incr[i] = gdim.dims[i].is;
879                 }
880                 for (ir = 0; ir < gdim.howmany_rank; ir++)
881                 {
882                     j = 0;
883                     for (i = 1; i <= m; i++)
884                     {
885                         if (Ai == NULL)
886                         {
887                             r = check_ND_symmetry(Ar + j, NULL, ndims, dims, incr);
888                         }
889                         else
890                         {
891                             r = check_ND_symmetry(Ar + j, Ai + j, ndims, dims, incr);
892                         }
893                         if (r <= 0)
894                         {
895                             FREE(dims);
896                             FREE(dims1);
897                             FREE(incr);
898                             FREE(incr1);
899                             return r;
900                         }
901                         j += gdim.howmany_dims[0].is;
902                         for (k = gdim.howmany_rank - 2; k >= 0; k--)
903                         {
904                             if (i % dims1[k] == 0)
905                             {
906                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
907                                 break;
908                             }
909                         }
910                     }
911                 }
912                 FREE(dims);
913                 FREE(dims1);
914                 FREE(incr);
915                 FREE(incr1);
916                 return 1;
917         }
918     }
919     return 1;
920 }
921 /*--------------------------------------------------------------------------
922  * "symmetrizing" a vector A of length nA modifying the second half part of the vector
923  * nA even: A=[a0 A1 conj(A1($:-1:1))]
924  * nA odd : A=[a0 A1 am conj(A1($:-1:1))]
925  */
926 
complete_1D_array(double * Ar,double * Ai,int nA,int iA)927 void complete_1D_array(double *Ar, double *Ai, int nA, int iA)
928 {
929 
930     if (nA > 2)
931     {
932         int nAs2 = (int)(nA / 2);
933         int n = (nA % 2 == 0) ? nAs2 - 1 : nAs2;
934         int l1 = iA; /* ignore first element */
935         int l2 = (nA - 1) * iA;
936         int i = 0;
937         if (Ai == NULL)
938         {
939             for (i = 0; i < n; i++)
940             {
941                 Ar[l2] = Ar[l1];
942                 l1 += iA;
943                 l2 -= iA;
944             }
945         }
946         else
947         {
948             for (i = 0; i < n; i++)
949             {
950                 Ar[l2] = Ar[l1];
951                 Ai[l2] = -Ai[l1];
952                 l1 += iA;
953                 l2 -= iA;
954             }
955         }
956     }
957 }
958 
959 /*--------------------------------------------------------------------------
960  * "symmetrizing" a mA x nA array modifying the second half part of the columns
961  * nA even: A=[a11 A12 conj(A12($:-1:1))
962  *             A21 A22 conj(A22($:-1:1,$:-1:1))]
963  *
964  * nA odd : A=[a11 A12 am  conj(A12($:-1:1))
965  A21 A22 A2m conj(A22($:-1:1,$:-1:1))]]
966 */
967 
complete_2D_array(double * Ar,double * Ai,int mA,int iA,int nA,int jA)968 void complete_2D_array(double *Ar, double *Ai, int mA, int iA, int nA, int jA)
969 {
970     if (nA > 2)
971     {
972         int n = (nA % 2 == 0) ? (int)(nA / 2) - 1 : (int)(nA / 2);
973         int i = 0, j = 0; /* loop variables */
974         int l1 = jA + iA; /* the index of the first element of the A22 block A(2,2)*/
975         int l2 = (mA - 1) * iA + (nA - 1) * jA; /* the index of the last element of the A22 block A(mA,nA)*/
976         int step = 0;
977         /* first column  */
978         /*could not be useful because fftw only skip half of the rightmost dimension but it may be not exactly symmetric */
979 
980         complete_1D_array(Ar, Ai, mA, iA);
981 
982         /* first row */
983         complete_1D_array(Ar, Ai, nA, jA);
984 
985         /* A22 block */
986         if (Ai == NULL)
987         {
988             for (j = 0; j < n; j++)
989             {
990                 for (i = 1; i < mA; i++)
991                 {
992                     Ar[l2] = Ar[l1];
993                     l1 += iA;
994                     l2 -= iA;
995                 }
996                 step = -(mA - 1) * iA + jA;
997                 l1 += step;
998                 l2 -= step;
999             }
1000         }
1001         else
1002         {
1003             for (j = 0; j < n; j++)
1004             {
1005                 for (i = 1; i < mA; i++)
1006                 {
1007                     Ar[l2] = Ar[l1];
1008                     Ai[l2] = -Ai[l1];
1009                     l1 += iA;
1010                     l2 -= iA;
1011                 }
1012                 step = -(mA - 1) * iA + jA;
1013                 l1 += step;
1014                 l2 -= step;
1015             }
1016         }
1017     }
1018 }
1019 
complete_ND_array(double * Ar,double * Ai,int ndims,int * dims,int * incr)1020 int complete_ND_array(double *Ar, double *Ai, int ndims, int *dims, int *incr)
1021 {
1022     int i = 0, j = 0, l = 0;
1023     int r = 0;
1024     int l1 = 0;/* current 1-D index in array*/
1025     int l2 = 0;/* associated symmetric value 1-D index */
1026 
1027     int *temp = NULL;
1028     int *dims1 = NULL;
1029     int *incr1 = NULL;
1030     int nSub = 0, nSubs2 = 0, step = 0, k = 0;
1031 
1032     if (ndims == 2)
1033     {
1034         complete_2D_array(Ar, Ai, dims[0], incr[0], dims[1], incr[1]);
1035         return 0;
1036     }
1037     else if (ndims == 1)
1038     {
1039         complete_1D_array(Ar, Ai, dims[0], incr[0]);
1040         return 0;
1041     }
1042     if  ((temp =  (int *)MALLOC(sizeof(int) * (2 * ndims))) == NULL)
1043     {
1044         return -1;
1045     }
1046     dims1 = temp;
1047     incr1 = temp + ndims;
1048 
1049     for (i = 0; i < ndims; i++)
1050     {
1051         /* remove current dimension and increment out of  dims ans incr */
1052         l = 0;
1053         for   (j = 0; j < ndims; j++)
1054         {
1055             if (j != i)
1056             {
1057                 dims1[l] = dims[j];
1058                 incr1[l] = incr[j];
1059                 l++;
1060             }
1061         }
1062         r = complete_ND_array(Ar, Ai, ndims - 1, dims1, incr1);
1063         if (r < 0)
1064         {
1065             dims1 = NULL;
1066             incr1 = NULL;
1067             FREE(temp);
1068             return r;
1069         }
1070     }
1071 
1072     /* check bloc A(2:$,....,2:$)*/
1073     l1 = 0;
1074     for (i = 0; i < ndims; i++)
1075     {
1076         l1 += incr[i];
1077     }
1078     /*A($,...,$) index*/
1079     l2 = 0;
1080     for (i = 0; i < ndims; i++)
1081     {
1082         l2 += (dims[i] - 1) * incr[i];
1083     }
1084 
1085 
1086     /* cumprod(size(A(2:$,....,2:$)))*/
1087     incr1[0] = dims[0] - 1;
1088     for (i = 1; i < (ndims - 1); i++)
1089     {
1090         incr1[i] = incr1[i - 1] * (dims[i] - 1) ;
1091     }
1092     /* steps*/
1093     dims1[0] = (dims[0] - 2) * incr[0];
1094     for (i = 1; i < (ndims - 1); i++)
1095     {
1096         dims1[i] = dims1[i - 1] + (dims[i] - 2) * incr[i];
1097     }
1098 
1099     /*  A(2:$,....,2:$) block number of elements*/
1100     nSub = 1;
1101     for (i = 0; i < ndims; i++)
1102     {
1103         nSub *= (dims[i] - 1);
1104     }
1105 
1106     nSubs2 = (int)(nSub / 2);
1107 
1108     if (Ai == 0)
1109     {
1110         /* Real case */
1111         for (i = 0; i < nSubs2; i++)
1112         {
1113             Ar[l2] = Ar[l1];
1114             step = incr[0];
1115             for (j = ndims - 2; j >= 0; j--)
1116             {
1117                 if ((i + 1) % incr1[j] == 0)
1118                 {
1119                     step = -dims1[j] + incr[j + 1] ;
1120                     break;
1121                 }
1122             }
1123             l1 += step;
1124             l2 -= step;
1125         }
1126     }
1127     else    /* Complex case */
1128     {
1129         for (i = 0; i < nSubs2; i++)
1130         {
1131             Ar[l2] = Ar[l1];
1132             Ai[l2] = -Ai[l1];
1133             step = incr[0];
1134             for (j = ndims - 2; j >= 0; j--)
1135             {
1136                 if ((i + 1) % incr1[j] == 0)
1137                 {
1138                     step = -dims1[j] + incr[j + 1] ;
1139                     break;
1140                 }
1141             }
1142             l1 += step;
1143             l2 -= step;
1144         }
1145     }
1146     dims1 = NULL;
1147     incr1 = NULL;
1148     FREE(temp);
1149     return 1;
1150 }
1151 
complete_array(double * Ar,double * Ai,guru_dim_struct gdim)1152 int complete_array(double *Ar, double *Ai, guru_dim_struct gdim)
1153 {
1154     int ndims = gdim.rank;
1155     int * dims = NULL;
1156     int * incr = NULL;
1157     int r = -1;
1158     int i = 0, j = 0, k = 0;
1159     if (gdim.howmany_rank == 0)
1160     {
1161         switch (gdim.rank)
1162         {
1163             case 1:
1164                 complete_1D_array(Ar, Ai, gdim.dims[0].n, gdim.dims[0].is);
1165                 return 0;
1166             case 2:
1167                 complete_2D_array(Ar, Ai, gdim.dims[0].n, gdim.dims[0].is, gdim.dims[1].n, gdim.dims[1].is);
1168                 return 0;
1169             default: /*general N-D case*/
1170                 if ((dims = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1171                 {
1172                     return -1;
1173                 }
1174                 if ((incr = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1175                 {
1176                     FREE(dims);
1177                     return -1;
1178                 }
1179                 for (i = 0; i < ndims; i++)
1180                 {
1181                     dims[i] = gdim.dims[i].n;
1182                     incr[i] = gdim.dims[i].is;
1183                 }
1184                 r = complete_ND_array(Ar, Ai, ndims, dims, incr);
1185                 FREE(dims);
1186                 FREE(incr);
1187                 return r;
1188         }
1189     }
1190     else
1191     {
1192         int m = 0;
1193         int *dims1 = NULL;
1194         int *incr1 = NULL;
1195         int hrank = gdim.howmany_rank;
1196 
1197         if ((dims1 = (int *)MALLOC(sizeof(int) * hrank)) == NULL)
1198         {
1199             return -1;
1200         }
1201         dims1[0] = gdim.howmany_dims[0].n;
1202         for (i = 1; i < hrank; i++)
1203         {
1204             dims1[i] = dims1[i - 1] * gdim.howmany_dims[i].n;
1205         }
1206         m = dims1[gdim.howmany_rank - 1];
1207 
1208         if ((incr1 = (int *)MALLOC(sizeof(int) * hrank)) == NULL)
1209         {
1210             FREE(dims1);
1211             return -1;
1212         }
1213         incr1[0] = gdim.howmany_dims[0].n * gdim.howmany_dims[0].is;
1214         for (i = 1; i < hrank; i++)
1215         {
1216             incr1[i] = incr1[i - 1] + (gdim.howmany_dims[i].n - 1) * gdim.howmany_dims[i].is;;
1217         }
1218         switch (gdim.rank)
1219         {
1220             case 1: /* multiple 1D fft */
1221                 if (Ai == NULL)
1222                 {
1223                     j = 0;
1224                     for (i = 1; i <= m; i++)
1225                     {
1226                         complete_1D_array(Ar + j, NULL, gdim.dims[0].n, gdim.dims[0].is);
1227                         j += gdim.howmany_dims[0].is;
1228                         for (k = hrank - 2; k >= 0; k--)
1229                         {
1230                             if (i % dims1[k] == 0)
1231                             {
1232                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1233                                 break;
1234                             }
1235                         }
1236                     }
1237                 }
1238                 else
1239                 {
1240                     j = 0;
1241                     for (i = 1; i <= m; i++)
1242                     {
1243                         complete_1D_array(Ar + j, Ai + j, gdim.dims[0].n, gdim.dims[0].is);
1244                         j += gdim.howmany_dims[0].is;
1245                         for (k = hrank - 2; k >= 0; k--)
1246                         {
1247                             if (i % dims1[k] == 0)
1248                             {
1249                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1250                                 break;
1251                             }
1252                         }
1253                     }
1254                 }
1255                 FREE(dims1);
1256                 FREE(incr1);
1257                 return 0;
1258             case 2: /* multiple 2D fft */
1259                 if (Ai == NULL)
1260                 {
1261                     j = 0;
1262                     for (i = 1; i <= m; i++)
1263                     {
1264                         complete_2D_array(Ar + j, NULL, gdim.dims[0].n, gdim.dims[0].is, gdim.dims[1].n, gdim.dims[1].is);
1265                         j += gdim.howmany_dims[0].is;
1266                         for (k = hrank - 2; k >= 0; k--)
1267                         {
1268                             if (i % dims1[k] == 0)
1269                             {
1270                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1271                                 break;
1272                             }
1273                         }
1274                     }
1275                 }
1276                 else
1277                 {
1278                     j = 0;
1279                     for (i = 1; i <= m; i++)
1280                     {
1281                         complete_2D_array(Ar + j, Ai + j, gdim.dims[0].n, gdim.dims[0].is, gdim.dims[1].n, gdim.dims[1].is);
1282 
1283                         j += gdim.howmany_dims[0].is;
1284                         for (k = hrank - 2; k >= 0; k--)
1285                         {
1286                             if (i % dims1[k] == 0)
1287                             {
1288                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1289                                 break;
1290                             }
1291                         }
1292                     }
1293                 }
1294                 FREE(dims1);
1295                 FREE(incr1);
1296                 return 0;
1297             default:  /* multiple ND fft */
1298                 if ((dims = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1299                 {
1300                     FREE(dims1);
1301                     FREE(incr1);
1302                     return -1;
1303                 }
1304                 if ((incr = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1305                 {
1306                     FREE(dims);
1307                     FREE(dims1);
1308                     FREE(incr1);
1309                     return -1;
1310                 }
1311                 for (i = 0; i < ndims; i++)
1312                 {
1313                     dims[i] = gdim.dims[i].n;
1314                     incr[i] = gdim.dims[i].is;
1315                 }
1316                 j = 0;
1317                 for (i = 1; i <= m; i++)
1318                 {
1319                     if (Ai == NULL)
1320                     {
1321                         r = complete_ND_array(Ar + j, NULL, ndims, dims, incr);
1322                     }
1323                     else
1324                     {
1325                         r = complete_ND_array(Ar + j, Ai + j, ndims, dims, incr);
1326                     }
1327                     if (r < 0)
1328                     {
1329                         FREE(dims);
1330                         FREE(dims1);
1331                         FREE(incr);
1332                         FREE(incr1);
1333                         return r;
1334                     }
1335                     j += gdim.howmany_dims[0].is;
1336                     for (k = hrank - 2; k >= 0; k--)
1337                     {
1338                         if (i % dims1[k] == 0)
1339                         {
1340                             j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1341                             break;
1342                         }
1343                     }
1344                 }
1345                 FREE(dims);
1346                 FREE(dims1);
1347                 FREE(incr);
1348                 FREE(incr1);
1349         }
1350     }
1351     return 0;
1352 }
1353 /*--------------------------------------------------------------------------
1354  * Check if Scilab is linked with MKL library * Some fftw functions
1355  * are not yet implemented in MKL in particular wisdom; guru_split real case
1356  * functions and  guru_split complex with homany_rank>1
1357  */
1358 
withMKL(void)1359 int withMKL(void)
1360 {
1361     static int iWithMKL = -1;
1362     if (iWithMKL == -1)
1363     {
1364         char* str = NULL;
1365         iWithMKL = 1;
1366         str = call_fftw_export_wisdom_to_string();
1367         if (str)
1368         {
1369             iWithMKL = 0;
1370             // According to the FFTW documentation we should free str
1371             // string but doing makes Scilab crash!?
1372             //free(str);
1373         }
1374     }
1375 
1376     return iWithMKL;
1377 }/*--------------------------------------------------------------------------*/
1378 
1379 
1380 
1381 
dct_scale_1D_array(double * Ar,double * Ai,int nA,int iA,int isn,double fact)1382 void dct_scale_1D_array(double *Ar, double *Ai, int nA, int iA, int isn, double fact)
1383 {
1384     /* fact: multiplication factor for all terms but the first one*/
1385     double s, s0;
1386     int i = 0;
1387 
1388     if (isn == -1)
1389     {
1390         s0 = fact * 0.5 / sqrt(nA);
1391     }
1392     else
1393     {
1394         s0 = fact / sqrt(nA); /* 2.0*sqrt(nA)/(2*nA) */
1395     }
1396     s = fact / sqrt(2.0 * nA); /* sqrt(2.0*nA)/(2*nA) */
1397     if (Ai == NULL)
1398     {
1399         Ar[0] *= s0;
1400         for (i = 1; i < nA; i++)
1401         {
1402             Ar[i * iA] *= s;
1403         }
1404     }
1405     else
1406     {
1407         Ar[0] *= s0;
1408         Ai[0] *= s0;
1409         for (i = 1; i < nA; i++)
1410         {
1411             Ar[i * iA] *= s;
1412             Ai[i * iA] *= s;
1413         }
1414 
1415     }
1416 }
1417 
1418 
1419 
dct_scale_2D_array(double * Ar,double * Ai,int mA,int iA,int nA,int jA,int isn,double fact)1420 void dct_scale_2D_array(double *Ar, double *Ai, int mA, int iA, int nA, int jA, int isn, double fact)
1421 {
1422     int j = 0; /* loop variables */
1423     double s, s0;
1424     s = fact / sqrt(2 * nA);
1425     s0 = fact / sqrt(nA);
1426     if (isn == -1)
1427     {
1428         s0 *= 0.5;
1429     }
1430 
1431     /* first column  */
1432     dct_scale_1D_array(Ar, Ai, mA, iA, isn, s0);
1433     /* other columns */
1434     if (Ai == NULL)
1435     {
1436         for (j = 1; j < nA; j++)
1437         {
1438             dct_scale_1D_array(&Ar[j * jA], NULL, mA, iA, isn, s);
1439         }
1440     }
1441     else
1442     {
1443         for (j = 1; j < nA; j++)
1444         {
1445             dct_scale_1D_array(&Ar[j * jA], &Ai[j * jA], mA, iA, isn, s);
1446         }
1447     }
1448 }
1449 
dct_scale_ND_array(double * Ar,double * Ai,int ndims,int * dims,int * incr,int isn,double fact)1450 int dct_scale_ND_array(double *Ar, double *Ai, int ndims, int *dims, int *incr, int isn, double fact)
1451 {
1452     int i = 0;
1453     double s = 1.0, s0 = 1.0;
1454 
1455     if (ndims == 2)
1456     {
1457         dct_scale_2D_array(Ar, Ai, dims[0], incr[0], dims[1], incr[1], isn, fact);
1458     }
1459     else if (ndims == 1)
1460     {
1461         dct_scale_1D_array(Ar, Ai, dims[0], incr[0], isn, fact);
1462     }
1463     else
1464     {
1465         /* Decompose recursively along the first array dimension
1466            A_scaled(i,:,...,:)=s1(i)*scale(A(i,:,...,:))
1467            with
1468            s1(1) = 1/(2*sqrt(n1) and  s1(i>1) = 1/(sqrt(2*n1)
1469         */
1470         s = fact / sqrt(2.0 * dims[0]);
1471         s0 = fact / sqrt(dims[0]);
1472         if (isn == -1)
1473         {
1474             s0 *= 0.5;
1475         }
1476 
1477         if (Ai == NULL)
1478         {
1479             /* first index: s1(1)*/
1480             dct_scale_ND_array(Ar, Ai, ndims - 1, dims + 1, incr + 1, isn, s0);
1481             /* next indexes: s1(i>1)*/
1482             for (i = 1; i < dims[0]; i++)
1483             {
1484                 dct_scale_ND_array(&Ar[i * incr[0]], NULL, ndims - 1, dims + 1, incr + 1, isn, s);
1485             }
1486         }
1487         else
1488         {
1489             dct_scale_ND_array(Ar, Ai, ndims - 1, dims + 1, incr + 1, isn, s0);
1490 
1491             for (i = 1; i < dims[0]; i++)
1492             {
1493                 dct_scale_ND_array(&Ar[i * incr[0]], &Ai[i * incr[0]], ndims - 1, dims + 1, incr + 1, isn, s);
1494             }
1495         }
1496     }
1497     return 0;
1498 
1499 }
1500 
1501 
dct_scale_array(double * Ar,double * Ai,guru_dim_struct gdim,int isn)1502 int dct_scale_array(double *Ar, double *Ai, guru_dim_struct gdim, int isn)
1503 {
1504     int * dims = NULL;
1505     int * incr = NULL;
1506     int *dims1 = NULL;
1507     int *incr1 = NULL;
1508 
1509     int i = 0, j = 0, k = 0;
1510     if (gdim.howmany_rank == 0)
1511     {
1512         switch (gdim.rank)
1513         {
1514             case 1:
1515                 dct_scale_1D_array(Ar, Ai, gdim.dims[0].n, gdim.dims[0].is, isn, (double)1.0);
1516                 return 0;
1517             case 2:
1518                 dct_scale_2D_array(Ar, Ai, gdim.dims[0].n, gdim.dims[0].is, gdim.dims[1].n, gdim.dims[1].is, isn, (double)1.0);
1519                 return 0;
1520             default: /*general N-D case*/
1521                 if ((dims = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1522                 {
1523                     goto ERR;
1524                 }
1525                 if ((incr = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1526                 {
1527                     goto ERR;
1528                 }
1529                 for (i = 0; i < gdim.rank; i++)
1530                 {
1531                     dims[i] = gdim.dims[i].n;
1532                     incr[i] = gdim.dims[i].is;
1533                 }
1534                 dct_scale_ND_array(Ar, Ai, gdim.rank, dims, incr, isn, (double)1.0);
1535         }
1536     }
1537     else
1538     {
1539         int m = 0;
1540         int hrank = gdim.howmany_rank;
1541         if ((dims1 = (int *)MALLOC(sizeof(int) * hrank)) == NULL)
1542         {
1543             goto ERR;
1544         }
1545         dims1[0] = gdim.howmany_dims[0].n;
1546         for (i = 1; i < hrank; i++)
1547         {
1548             dims1[i] = dims1[i - 1] * gdim.howmany_dims[i].n;
1549         }
1550         m = dims1[gdim.howmany_rank - 1];
1551 
1552         if ((incr1 = (int *)MALLOC(sizeof(int) * hrank)) == NULL)
1553         {
1554             goto ERR;
1555         }
1556 
1557         incr1[0] = gdim.howmany_dims[0].n * gdim.howmany_dims[0].is;
1558         for (i = 1; i < hrank; i++)
1559         {
1560             incr1[i] = incr1[i - 1] + (gdim.howmany_dims[i].n - 1) * gdim.howmany_dims[i].is;;
1561         }
1562         switch (gdim.rank)
1563         {
1564             case 1: /* multiple 1D dct */
1565                 if (Ai == NULL)
1566                 {
1567                     j = 0;
1568                     for (i = 1; i <= m; i++)
1569                     {
1570                         dct_scale_1D_array(Ar + j, NULL, gdim.dims[0].n, gdim.dims[0].is, isn, (double)1.0);
1571                         j += gdim.howmany_dims[0].is;
1572                         for (k = hrank - 2; k >= 0; k--)
1573                         {
1574                             if (i % dims1[k] == 0)
1575                             {
1576                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1577                                 break;
1578                             }
1579                         }
1580                     }
1581                 }
1582                 else
1583                 {
1584                     j = 0;
1585                     for (i = 1; i <= m; i++)
1586                     {
1587                         dct_scale_1D_array(Ar + j, Ai + j, gdim.dims[0].n, gdim.dims[0].is, isn, (double)1.0);
1588                         j += gdim.howmany_dims[0].is;
1589                         for (k = hrank - 2; k >= 0; k--)
1590                         {
1591                             if (i % dims1[k] == 0)
1592                             {
1593                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1594                                 break;
1595                             }
1596                         }
1597                     }
1598                 }
1599                 break;
1600             case 2: /* multiple 2D dct */
1601                 if (Ai == NULL)
1602                 {
1603                     j = 0;
1604                     for (i = 1; i <= m; i++)
1605                     {
1606                         dct_scale_2D_array(&Ar[j], NULL, gdim.dims[0].n, gdim.dims[0].is, gdim.dims[1].n, gdim.dims[1].is, isn, (double)1.0);
1607                         j += gdim.howmany_dims[0].is;
1608                         for (k = hrank - 2; k >= 0; k--)
1609                         {
1610                             if (i % dims1[k] == 0)
1611                             {
1612                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1613                                 break;
1614                             }
1615                         }
1616                     }
1617                 }
1618 
1619                 else
1620                 {
1621                     j = 0;
1622                     for (i = 1; i <= m; i++)
1623                     {
1624                         dct_scale_2D_array(&Ar[j], &Ai[j], gdim.dims[0].n, gdim.dims[0].is, gdim.dims[1].n, gdim.dims[1].is, isn, (double)1.0);
1625 
1626                         j += gdim.howmany_dims[0].is;
1627                         for (k = hrank - 2; k >= 0; k--)
1628                         {
1629                             if (i % dims1[k] == 0)
1630                             {
1631                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1632                                 break;
1633                             }
1634                         }
1635                     }
1636                 }
1637                 break;
1638             default:  /* multiple ND dct */
1639                 if ((dims = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1640                 {
1641                     goto ERR;
1642                 }
1643                 if ((incr = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1644                 {
1645                     goto ERR;
1646                 }
1647 
1648                 for (i = 0; i < gdim.rank; i++)
1649                 {
1650                     dims[i] = gdim.dims[i].n;
1651                     incr[i] = gdim.dims[i].is;
1652                 }
1653                 j = 0;
1654                 for (i = 1; i <= m; i++)
1655                 {
1656                     if (Ai == NULL)
1657                     {
1658                         dct_scale_ND_array(Ar + j, NULL, gdim.rank, dims, incr, isn, (double)1.0);
1659                     }
1660                     else
1661                     {
1662                         dct_scale_ND_array(Ar + j, Ai + j, gdim.rank, dims, incr, isn, (double)1.0);
1663                     }
1664 
1665                     j += gdim.howmany_dims[0].is;
1666                     for (k = hrank - 2; k >= 0; k--)
1667                     {
1668                         if (i % dims1[k] == 0)
1669                         {
1670                             j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1671                             break;
1672                         }
1673                     }
1674                 }
1675         }
1676 
1677     }
1678     FREE(dims);
1679     FREE(incr);
1680     FREE(dims1);
1681     FREE(incr1);
1682     return 0;
1683 
1684 ERR:
1685     FREE(dims);
1686     FREE(incr);
1687     FREE(dims1);
1688     FREE(incr1);
1689     return -1;
1690 }
1691 
dst_scale_1D_array(double * Ar,double * Ai,int nA,int iA,int isn,double fact)1692 void dst_scale_1D_array(double *Ar, double *Ai, int nA, int iA, int isn, double fact)
1693 {
1694     /* fact: multiplication factor for all terms but the first one*/
1695     double s = fact / (1.0 + nA);
1696     int i = 0;
1697 
1698     if (Ai == NULL)
1699     {
1700         for (i = 0; i < nA; i++)
1701         {
1702             Ar[i * iA] *= s;
1703         }
1704     }
1705     else
1706     {
1707         for (i = 0; i < nA; i++)
1708         {
1709             Ar[i * iA] *= s;
1710             Ai[i * iA] *= s;
1711         }
1712 
1713     }
1714 }
1715 
1716 
dst_scale_2D_array(double * Ar,double * Ai,int mA,int iA,int nA,int jA,int isn,double fact)1717 void dst_scale_2D_array(double *Ar, double *Ai, int mA, int iA, int nA, int jA, int isn, double fact)
1718 {
1719     int j = 0; /* loop variables */
1720     double s = fact / (1.0 + nA);
1721 
1722     if (Ai == NULL)
1723     {
1724         for (j = 0; j < nA; j++)
1725         {
1726             dst_scale_1D_array(&Ar[j * jA], NULL, mA, iA, isn, s);
1727         }
1728     }
1729     else
1730     {
1731         for (j = 0; j < nA; j++)
1732         {
1733             dst_scale_1D_array(&Ar[j * jA], &Ai[j * jA], mA, iA, isn, s);
1734         }
1735     }
1736 }
1737 
dst_scale_ND_array(double * Ar,double * Ai,int ndims,int * dims,int * incr,int isn,double fact)1738 int dst_scale_ND_array(double *Ar, double *Ai, int ndims, int *dims, int *incr, int isn, double fact)
1739 {
1740     int i = 0;
1741     double s = 1.0;
1742 
1743     if (ndims == 2)
1744     {
1745         dst_scale_2D_array(Ar, Ai, dims[0], incr[0], dims[1], incr[1], isn, fact);
1746     }
1747     else if (ndims == 1)
1748     {
1749         dst_scale_1D_array(Ar, Ai, dims[0], incr[0], isn, fact);
1750     }
1751     else
1752     {
1753         /* Decompose recursively along the first array dimension
1754            A_scaled(i,:,...,:)=s1*scale(A(i,:,...,:))
1755            with
1756            s1 = 1/(n+1)
1757         */
1758 
1759         s = fact / (1.0 + dims[0]);
1760 
1761         if (Ai == NULL)
1762         {
1763             /* next indexes: s1(i>1)*/
1764             for (i = 0; i < dims[0]; i++)
1765             {
1766                 dst_scale_ND_array(&Ar[i * incr[0]], NULL, ndims - 1, dims + 1, incr + 1, isn, s);
1767             }
1768         }
1769         else
1770         {
1771             for (i = 0; i < dims[0]; i++)
1772             {
1773                 dst_scale_ND_array(&Ar[i * incr[0]], &Ai[i * incr[0]], ndims - 1, dims + 1, incr + 1, isn, s);
1774             }
1775         }
1776     }
1777     return 0;
1778 }
1779 
1780 
dst_scale_array(double * Ar,double * Ai,guru_dim_struct gdim,int isn)1781 int dst_scale_array(double *Ar, double *Ai, guru_dim_struct gdim, int isn)
1782 {
1783     int * dims = NULL;
1784     int * incr = NULL;
1785     int *dims1 = NULL;
1786     int *incr1 = NULL;
1787 
1788     int i = 0, j = 0, k = 0;
1789     if (gdim.howmany_rank == 0)
1790     {
1791         switch (gdim.rank)
1792         {
1793             case 1:
1794                 dst_scale_1D_array(Ar, Ai, gdim.dims[0].n, gdim.dims[0].is, isn, (double)1.0);
1795                 return 0;
1796             case 2:
1797                 dst_scale_2D_array(Ar, Ai, gdim.dims[0].n, gdim.dims[0].is, gdim.dims[1].n, gdim.dims[1].is, isn, (double)1.0);
1798                 return 0;
1799             default: /*general N-D case*/
1800                 if ((dims = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1801                 {
1802                     goto ERR;
1803                 }
1804                 if ((incr = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1805                 {
1806                     goto ERR;
1807                 }
1808                 for (i = 0; i < gdim.rank; i++)
1809                 {
1810                     dims[i] = gdim.dims[i].n;
1811                     incr[i] = gdim.dims[i].is;
1812                 }
1813                 dst_scale_ND_array(Ar, Ai, gdim.rank, dims, incr, isn, (double)1.0);
1814         }
1815     }
1816     else
1817     {
1818         int m = 0;
1819         int hrank = gdim.howmany_rank;
1820 
1821         if ((dims1 = (int *)MALLOC(sizeof(int) * hrank)) == NULL)
1822         {
1823             goto ERR;
1824         }
1825         dims1[0] = gdim.howmany_dims[0].n;
1826         for (i = 1; i < hrank; i++)
1827         {
1828             dims1[i] = dims1[i - 1] * gdim.howmany_dims[i].n;
1829         }
1830         m = dims1[gdim.howmany_rank - 1];
1831 
1832         if ((incr1 = (int *)MALLOC(sizeof(int) * hrank)) == NULL)
1833         {
1834             goto ERR;
1835         }
1836 
1837         incr1[0] = gdim.howmany_dims[0].n * gdim.howmany_dims[0].is;
1838         for (i = 1; i < hrank; i++)
1839         {
1840             incr1[i] = incr1[i - 1] + (gdim.howmany_dims[i].n - 1) * gdim.howmany_dims[i].is;;
1841         }
1842 
1843         switch (gdim.rank)
1844         {
1845             case 1: /* multiple 1D dst */
1846                 if (Ai == NULL)
1847                 {
1848                     j = 0;
1849                     for (i = 1; i <= m; i++)
1850                     {
1851                         dst_scale_1D_array(Ar + j, NULL, gdim.dims[0].n, gdim.dims[0].is, isn, (double)1.0);
1852                         j += gdim.howmany_dims[0].is;
1853                         for (k = hrank - 2; k >= 0; k--)
1854                         {
1855                             if (i % dims1[k] == 0)
1856                             {
1857                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1858                                 break;
1859                             }
1860                         }
1861                     }
1862                 }
1863                 else
1864                 {
1865                     j = 0;
1866                     for (i = 1; i <= m; i++)
1867                     {
1868                         dst_scale_1D_array(Ar + j, Ai + j, gdim.dims[0].n, gdim.dims[0].is, isn, (double)1.0);
1869                         j += gdim.howmany_dims[0].is;
1870                         for (k = hrank - 2; k >= 0; k--)
1871                         {
1872                             if (i % dims1[k] == 0)
1873                             {
1874                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1875                                 break;
1876                             }
1877                         }
1878                     }
1879                 }
1880                 break;
1881             case 2: /* multiple 2D dst */
1882                 if (Ai == NULL)
1883                 {
1884                     j = 0;
1885                     for (i = 1; i <= m; i++)
1886                     {
1887                         dst_scale_2D_array(&Ar[j], NULL, gdim.dims[0].n, gdim.dims[0].is, gdim.dims[1].n, gdim.dims[1].is, isn, (double)1.0);
1888                         j += gdim.howmany_dims[0].is;
1889                         for (k = hrank - 2; k >= 0; k--)
1890                         {
1891                             if (i % dims1[k] == 0)
1892                             {
1893                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1894                                 break;
1895                             }
1896                         }
1897                     }
1898                 }
1899 
1900                 else
1901                 {
1902                     j = 0;
1903                     for (i = 1; i <= m; i++)
1904                     {
1905                         dst_scale_2D_array(&Ar[j], &Ai[j], gdim.dims[0].n, gdim.dims[0].is, gdim.dims[1].n, gdim.dims[1].is, isn, (double)1.0);
1906 
1907                         j += gdim.howmany_dims[0].is;
1908                         for (k = hrank - 2; k >= 0; k--)
1909                         {
1910                             if (i % dims1[k] == 0)
1911                             {
1912                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1913                                 break;
1914                             }
1915                         }
1916                     }
1917                 }
1918                 break;
1919             default:  /* multiple ND dst */
1920                 if ((dims = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1921                 {
1922                     goto ERR;
1923                 }
1924                 if ((incr = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1925                 {
1926                     goto ERR;
1927                 }
1928 
1929                 for (i = 0; i < gdim.rank; i++)
1930                 {
1931                     dims[i] = gdim.dims[i].n;
1932                     incr[i] = gdim.dims[i].is;
1933                 }
1934                 j = 0;
1935                 for (i = 1; i <= m; i++)
1936                 {
1937                     if (Ai == NULL)
1938                     {
1939                         dst_scale_ND_array(Ar + j, NULL, gdim.rank, dims, incr, isn, (double)1.0);
1940                     }
1941                     else
1942                     {
1943                         dst_scale_ND_array(Ar + j, Ai + j, gdim.rank, dims, incr, isn, (double)1.0);
1944                     }
1945 
1946                     j += gdim.howmany_dims[0].is;
1947                     for (k = hrank - 2; k >= 0; k--)
1948                     {
1949                         if (i % dims1[k] == 0)
1950                         {
1951                             j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1952                             break;
1953                         }
1954                     }
1955                 }
1956         }
1957 
1958     }
1959     FREE(dims);
1960     FREE(incr);
1961     FREE(dims1);
1962     FREE(incr1);
1963     return 0;
1964 
1965 ERR:
1966     FREE(dims);
1967     FREE(incr);
1968     FREE(dims1);
1969     FREE(incr1);
1970     return -1;
1971 }
1972