[bug_14225] Scilab binary return 1 when an error occured in scilab before exit. Retur...
[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     FREE(str);
1337     return iWithMKL;
1338 }
1339 /*--------------------------------------------------------------------------*/
1340
1341
1342
1343
1344 void dct_scale_1D_array(double *Ar, double *Ai, int nA, int iA, int isn, double fact)
1345 {
1346     /* fact: multiplication factor for all terms but the first one*/
1347     double s, s0;
1348     int i = 0;
1349
1350     if (isn == -1)
1351     {
1352         s0 = fact * 0.5 / sqrt(nA);
1353     }
1354     else
1355     {
1356         s0 = fact / sqrt(nA); /* 2.0*sqrt(nA)/(2*nA) */
1357     }
1358     s = fact / sqrt(2.0 * nA); /* sqrt(2.0*nA)/(2*nA) */
1359     if (Ai == NULL)
1360     {
1361         Ar[0] *= s0;
1362         for (i = 1; i < nA; i++)
1363         {
1364             Ar[i * iA] *= s;
1365         }
1366     }
1367     else
1368     {
1369         Ar[0] *= s0;
1370         Ai[0] *= s0;
1371         for (i = 1; i < nA; i++)
1372         {
1373             Ar[i * iA] *= s;
1374             Ai[i * iA] *= s;
1375         }
1376
1377     }
1378 }
1379
1380
1381
1382 void dct_scale_2D_array(double *Ar, double *Ai, int mA, int iA, int nA, int jA, int isn, double fact)
1383 {
1384     int j = 0; /* loop variables */
1385     double s, s0;
1386     s = fact / sqrt(2 * nA);
1387     s0 = fact / sqrt(nA);
1388     if (isn == -1)
1389     {
1390         s0 *= 0.5;
1391     }
1392
1393     /* first column  */
1394     dct_scale_1D_array(Ar, Ai, mA, iA, isn, s0);
1395     /* other columns */
1396     if (Ai == NULL)
1397     {
1398         for (j = 1; j < nA; j++)
1399         {
1400             dct_scale_1D_array(&Ar[j * jA], NULL, mA, iA, isn, s);
1401         }
1402     }
1403     else
1404     {
1405         for (j = 1; j < nA; j++)
1406         {
1407             dct_scale_1D_array(&Ar[j * jA], &Ai[j * jA], mA, iA, isn, s);
1408         }
1409     }
1410 }
1411
1412 int dct_scale_ND_array(double *Ar, double *Ai, int ndims, int *dims, int *incr, int isn, double fact)
1413 {
1414     int i = 0;
1415     double s = 1.0, s0 = 1.0;
1416
1417     if (ndims == 2)
1418     {
1419         dct_scale_2D_array(Ar, Ai, dims[0], incr[0], dims[1], incr[1], isn, fact);
1420     }
1421     else if (ndims == 1)
1422     {
1423         dct_scale_1D_array(Ar, Ai, dims[0], incr[0], isn, fact);
1424     }
1425     else
1426     {
1427         /* Decompose recursively along the first array dimension
1428            A_scaled(i,:,...,:)=s1(i)*scale(A(i,:,...,:))
1429            with
1430            s1(1) = 1/(2*sqrt(n1) and  s1(i>1) = 1/(sqrt(2*n1)
1431         */
1432         s = fact / sqrt(2.0 * dims[0]);
1433         s0 = fact / sqrt(dims[0]);
1434         if (isn == -1)
1435         {
1436             s0 *= 0.5;
1437         }
1438
1439         if (Ai == NULL)
1440         {
1441             /* first index: s1(1)*/
1442             dct_scale_ND_array(Ar, Ai, ndims - 1, dims + 1, incr + 1, isn, s0);
1443             /* next indexes: s1(i>1)*/
1444             for (i = 1; i < dims[0]; i++)
1445             {
1446                 dct_scale_ND_array(&Ar[i * incr[0]], NULL, ndims - 1, dims + 1, incr + 1, isn, s);
1447             }
1448         }
1449         else
1450         {
1451             dct_scale_ND_array(Ar, Ai, ndims - 1, dims + 1, incr + 1, isn, s0);
1452
1453             for (i = 1; i < dims[0]; i++)
1454             {
1455                 dct_scale_ND_array(&Ar[i * incr[0]], &Ai[i * incr[0]], ndims - 1, dims + 1, incr + 1, isn, s);
1456             }
1457         }
1458     }
1459     return 0;
1460
1461 }
1462
1463
1464 int dct_scale_array(double *Ar, double *Ai, guru_dim_struct gdim, int isn)
1465 {
1466     int * dims = NULL;
1467     int * incr = NULL;
1468     int *dims1 = NULL;
1469     int *incr1 = NULL;
1470
1471     int i = 0, j = 0, k = 0;
1472     if (gdim.howmany_rank == 0)
1473     {
1474         switch (gdim.rank)
1475         {
1476             case 1:
1477                 dct_scale_1D_array(Ar, Ai, gdim.dims[0].n, gdim.dims[0].is, isn, (double)1.0);
1478                 return 0;
1479             case 2:
1480                 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);
1481                 return 0;
1482             default: /*general N-D case*/
1483                 if ((dims = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1484                 {
1485                     goto ERR;
1486                 }
1487                 if ((incr = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1488                 {
1489                     goto ERR;
1490                 }
1491                 for (i = 0; i < gdim.rank; i++)
1492                 {
1493                     dims[i] = gdim.dims[i].n;
1494                     incr[i] = gdim.dims[i].is;
1495                 }
1496                 dct_scale_ND_array(Ar, Ai, gdim.rank, dims, incr, isn, (double)1.0);
1497         }
1498     }
1499     else
1500     {
1501         int m = 0;
1502         int hrank = gdim.howmany_rank;
1503         if ((dims1 = (int *)MALLOC(sizeof(int) * hrank)) == NULL)
1504         {
1505             goto ERR;
1506         }
1507         dims1[0] = gdim.howmany_dims[0].n;
1508         for (i = 1; i < hrank; i++)
1509         {
1510             dims1[i] = dims1[i - 1] * gdim.howmany_dims[i].n;
1511         }
1512         m = dims1[gdim.howmany_rank - 1];
1513
1514         if ((incr1 = (int *)MALLOC(sizeof(int) * hrank)) == NULL)
1515         {
1516             goto ERR;
1517         }
1518
1519         incr1[0] = gdim.howmany_dims[0].n * gdim.howmany_dims[0].is;
1520         for (i = 1; i < hrank; i++)
1521         {
1522             incr1[i] = incr1[i - 1] + (gdim.howmany_dims[i].n - 1) * gdim.howmany_dims[i].is;;
1523         }
1524         switch (gdim.rank)
1525         {
1526             case 1: /* multiple 1D dct */
1527                 if (Ai == NULL)
1528                 {
1529                     j = 0;
1530                     for (i = 1; i <= m; i++)
1531                     {
1532                         dct_scale_1D_array(Ar + j, NULL, gdim.dims[0].n, gdim.dims[0].is, isn, (double)1.0);
1533                         j += gdim.howmany_dims[0].is;
1534                         for (k = hrank - 2; k >= 0; k--)
1535                         {
1536                             if (i % dims1[k] == 0)
1537                             {
1538                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1539                                 break;
1540                             }
1541                         }
1542                     }
1543                 }
1544                 else
1545                 {
1546                     j = 0;
1547                     for (i = 1; i <= m; i++)
1548                     {
1549                         dct_scale_1D_array(Ar + j, Ai + j, gdim.dims[0].n, gdim.dims[0].is, isn, (double)1.0);
1550                         j += gdim.howmany_dims[0].is;
1551                         for (k = hrank - 2; k >= 0; k--)
1552                         {
1553                             if (i % dims1[k] == 0)
1554                             {
1555                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1556                                 break;
1557                             }
1558                         }
1559                     }
1560                 }
1561                 break;
1562             case 2: /* multiple 2D dct */
1563                 if (Ai == NULL)
1564                 {
1565                     j = 0;
1566                     for (i = 1; i <= m; i++)
1567                     {
1568                         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);
1569                         j += gdim.howmany_dims[0].is;
1570                         for (k = hrank - 2; k >= 0; k--)
1571                         {
1572                             if (i % dims1[k] == 0)
1573                             {
1574                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1575                                 break;
1576                             }
1577                         }
1578                     }
1579                 }
1580
1581                 else
1582                 {
1583                     j = 0;
1584                     for (i = 1; i <= m; i++)
1585                     {
1586                         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);
1587
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             default:  /* multiple ND dct */
1601                 if ((dims = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1602                 {
1603                     goto ERR;
1604                 }
1605                 if ((incr = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1606                 {
1607                     goto ERR;
1608                 }
1609
1610                 for (i = 0; i < gdim.rank; i++)
1611                 {
1612                     dims[i] = gdim.dims[i].n;
1613                     incr[i] = gdim.dims[i].is;
1614                 }
1615                 j = 0;
1616                 for (i = 1; i <= m; i++)
1617                 {
1618                     if (Ai == NULL)
1619                     {
1620                         dct_scale_ND_array(Ar + j, NULL, gdim.rank, dims, incr, isn, (double)1.0);
1621                     }
1622                     else
1623                     {
1624                         dct_scale_ND_array(Ar + j, Ai + j, gdim.rank, dims, incr, isn, (double)1.0);
1625                     }
1626
1627                     j += gdim.howmany_dims[0].is;
1628                     for (k = hrank - 2; k >= 0; k--)
1629                     {
1630                         if (i % dims1[k] == 0)
1631                         {
1632                             j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1633                             break;
1634                         }
1635                     }
1636                 }
1637         }
1638
1639     }
1640     FREE(dims);
1641     FREE(incr);
1642     FREE(dims1);
1643     FREE(incr1);
1644     return 0;
1645
1646 ERR:
1647     FREE(dims);
1648     FREE(incr);
1649     FREE(dims1);
1650     FREE(incr1);
1651     return -1;
1652 }
1653
1654 void dst_scale_1D_array(double *Ar, double *Ai, int nA, int iA, int isn, double fact)
1655 {
1656     /* fact: multiplication factor for all terms but the first one*/
1657     double s = fact / (1.0 + nA);
1658     int i = 0;
1659
1660     if (Ai == NULL)
1661     {
1662         for (i = 0; i < nA; i++)
1663         {
1664             Ar[i * iA] *= s;
1665         }
1666     }
1667     else
1668     {
1669         for (i = 0; i < nA; i++)
1670         {
1671             Ar[i * iA] *= s;
1672             Ai[i * iA] *= s;
1673         }
1674
1675     }
1676 }
1677
1678
1679 void dst_scale_2D_array(double *Ar, double *Ai, int mA, int iA, int nA, int jA, int isn, double fact)
1680 {
1681     int j = 0; /* loop variables */
1682     double s = fact / (1.0 + nA);
1683
1684     if (Ai == NULL)
1685     {
1686         for (j = 0; j < nA; j++)
1687         {
1688             dst_scale_1D_array(&Ar[j * jA], NULL, mA, iA, isn, s);
1689         }
1690     }
1691     else
1692     {
1693         for (j = 0; j < nA; j++)
1694         {
1695             dst_scale_1D_array(&Ar[j * jA], &Ai[j * jA], mA, iA, isn, s);
1696         }
1697     }
1698 }
1699
1700 int dst_scale_ND_array(double *Ar, double *Ai, int ndims, int *dims, int *incr, int isn, double fact)
1701 {
1702     int i = 0;
1703     double s = 1.0;
1704
1705     if (ndims == 2)
1706     {
1707         dst_scale_2D_array(Ar, Ai, dims[0], incr[0], dims[1], incr[1], isn, fact);
1708     }
1709     else if (ndims == 1)
1710     {
1711         dst_scale_1D_array(Ar, Ai, dims[0], incr[0], isn, fact);
1712     }
1713     else
1714     {
1715         /* Decompose recursively along the first array dimension
1716            A_scaled(i,:,...,:)=s1*scale(A(i,:,...,:))
1717            with
1718            s1 = 1/(n+1)
1719         */
1720
1721         s = fact / (1.0 + dims[0]);
1722
1723         if (Ai == NULL)
1724         {
1725             /* next indexes: s1(i>1)*/
1726             for (i = 0; i < dims[0]; i++)
1727             {
1728                 dst_scale_ND_array(&Ar[i * incr[0]], NULL, ndims - 1, dims + 1, incr + 1, isn, s);
1729             }
1730         }
1731         else
1732         {
1733             for (i = 0; i < dims[0]; i++)
1734             {
1735                 dst_scale_ND_array(&Ar[i * incr[0]], &Ai[i * incr[0]], ndims - 1, dims + 1, incr + 1, isn, s);
1736             }
1737         }
1738     }
1739     return 0;
1740 }
1741
1742
1743 int dst_scale_array(double *Ar, double *Ai, guru_dim_struct gdim, int isn)
1744 {
1745     int * dims = NULL;
1746     int * incr = NULL;
1747     int *dims1 = NULL;
1748     int *incr1 = NULL;
1749
1750     int i = 0, j = 0, k = 0;
1751     if (gdim.howmany_rank == 0)
1752     {
1753         switch (gdim.rank)
1754         {
1755             case 1:
1756                 dst_scale_1D_array(Ar, Ai, gdim.dims[0].n, gdim.dims[0].is, isn, (double)1.0);
1757                 return 0;
1758             case 2:
1759                 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);
1760                 return 0;
1761             default: /*general N-D case*/
1762                 if ((dims = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1763                 {
1764                     goto ERR;
1765                 }
1766                 if ((incr = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1767                 {
1768                     goto ERR;
1769                 }
1770                 for (i = 0; i < gdim.rank; i++)
1771                 {
1772                     dims[i] = gdim.dims[i].n;
1773                     incr[i] = gdim.dims[i].is;
1774                 }
1775                 dst_scale_ND_array(Ar, Ai, gdim.rank, dims, incr, isn, (double)1.0);
1776         }
1777     }
1778     else
1779     {
1780         int m = 0;
1781         int hrank = gdim.howmany_rank;
1782
1783         if ((dims1 = (int *)MALLOC(sizeof(int) * hrank)) == NULL)
1784         {
1785             goto ERR;
1786         }
1787         dims1[0] = gdim.howmany_dims[0].n;
1788         for (i = 1; i < hrank; i++)
1789         {
1790             dims1[i] = dims1[i - 1] * gdim.howmany_dims[i].n;
1791         }
1792         m = dims1[gdim.howmany_rank - 1];
1793
1794         if ((incr1 = (int *)MALLOC(sizeof(int) * hrank)) == NULL)
1795         {
1796             goto ERR;
1797         }
1798
1799         incr1[0] = gdim.howmany_dims[0].n * gdim.howmany_dims[0].is;
1800         for (i = 1; i < hrank; i++)
1801         {
1802             incr1[i] = incr1[i - 1] + (gdim.howmany_dims[i].n - 1) * gdim.howmany_dims[i].is;;
1803         }
1804
1805         switch (gdim.rank)
1806         {
1807             case 1: /* multiple 1D dst */
1808                 if (Ai == NULL)
1809                 {
1810                     j = 0;
1811                     for (i = 1; i <= m; i++)
1812                     {
1813                         dst_scale_1D_array(Ar + j, NULL, gdim.dims[0].n, gdim.dims[0].is, isn, (double)1.0);
1814                         j += gdim.howmany_dims[0].is;
1815                         for (k = hrank - 2; k >= 0; k--)
1816                         {
1817                             if (i % dims1[k] == 0)
1818                             {
1819                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1820                                 break;
1821                             }
1822                         }
1823                     }
1824                 }
1825                 else
1826                 {
1827                     j = 0;
1828                     for (i = 1; i <= m; i++)
1829                     {
1830                         dst_scale_1D_array(Ar + j, Ai + j, gdim.dims[0].n, gdim.dims[0].is, isn, (double)1.0);
1831                         j += gdim.howmany_dims[0].is;
1832                         for (k = hrank - 2; k >= 0; k--)
1833                         {
1834                             if (i % dims1[k] == 0)
1835                             {
1836                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1837                                 break;
1838                             }
1839                         }
1840                     }
1841                 }
1842                 break;
1843             case 2: /* multiple 2D dst */
1844                 if (Ai == NULL)
1845                 {
1846                     j = 0;
1847                     for (i = 1; i <= m; i++)
1848                     {
1849                         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);
1850                         j += gdim.howmany_dims[0].is;
1851                         for (k = hrank - 2; k >= 0; k--)
1852                         {
1853                             if (i % dims1[k] == 0)
1854                             {
1855                                 j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1856                                 break;
1857                             }
1858                         }
1859                     }
1860                 }
1861
1862                 else
1863                 {
1864                     j = 0;
1865                     for (i = 1; i <= m; i++)
1866                     {
1867                         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);
1868
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             default:  /* multiple ND dst */
1882                 if ((dims = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1883                 {
1884                     goto ERR;
1885                 }
1886                 if ((incr = (int *)MALLOC(sizeof(int) * gdim.rank)) == NULL)
1887                 {
1888                     goto ERR;
1889                 }
1890
1891                 for (i = 0; i < gdim.rank; i++)
1892                 {
1893                     dims[i] = gdim.dims[i].n;
1894                     incr[i] = gdim.dims[i].is;
1895                 }
1896                 j = 0;
1897                 for (i = 1; i <= m; i++)
1898                 {
1899                     if (Ai == NULL)
1900                     {
1901                         dst_scale_ND_array(Ar + j, NULL, gdim.rank, dims, incr, isn, (double)1.0);
1902                     }
1903                     else
1904                     {
1905                         dst_scale_ND_array(Ar + j, Ai + j, gdim.rank, dims, incr, isn, (double)1.0);
1906                     }
1907
1908                     j += gdim.howmany_dims[0].is;
1909                     for (k = hrank - 2; k >= 0; k--)
1910                     {
1911                         if (i % dims1[k] == 0)
1912                         {
1913                             j += -incr1[k] + gdim.howmany_dims[k + 1].is;
1914                             break;
1915                         }
1916                     }
1917                 }
1918         }
1919
1920     }
1921     FREE(dims);
1922     FREE(incr);
1923     FREE(dims1);
1924     FREE(incr1);
1925     return 0;
1926
1927 ERR:
1928     FREE(dims);
1929     FREE(incr);
1930     FREE(dims1);
1931     FREE(incr1);
1932     return -1;
1933 }