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