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