Chris@82: /*
Chris@82:  * Copyright (c) 2003, 2007-14 Matteo Frigo
Chris@82:  * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
Chris@82:  *
Chris@82:  * This program is free software; you can redistribute it and/or modify
Chris@82:  * it under the terms of the GNU General Public License as published by
Chris@82:  * the Free Software Foundation; either version 2 of the License, or
Chris@82:  * (at your option) any later version.
Chris@82:  *
Chris@82:  * This program is distributed in the hope that it will be useful,
Chris@82:  * but WITHOUT ANY WARRANTY; without even the implied warranty of
Chris@82:  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Chris@82:  * GNU General Public License for more details.
Chris@82:  *
Chris@82:  * You should have received a copy of the GNU General Public License
Chris@82:  * along with this program; if not, write to the Free Software
Chris@82:  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
Chris@82:  *
Chris@82:  */
Chris@82: 
Chris@82: 
Chris@82: #include "rdft/rdft.h"
Chris@82: #include <stddef.h>
Chris@82: 
Chris@82: static void destroy(problem *ego_)
Chris@82: {
Chris@82:      problem_rdft *ego = (problem_rdft *) ego_;
Chris@82: #if !defined(STRUCT_HACK_C99) && !defined(STRUCT_HACK_KR)
Chris@82:      X(ifree0)(ego->kind);
Chris@82: #endif
Chris@82:      X(tensor_destroy2)(ego->vecsz, ego->sz);
Chris@82:      X(ifree)(ego_);
Chris@82: }
Chris@82: 
Chris@82: static void kind_hash(md5 *m, const rdft_kind *kind, int rnk)
Chris@82: {
Chris@82:      int i;
Chris@82:      for (i = 0; i < rnk; ++i)
Chris@82: 	  X(md5int)(m, kind[i]);
Chris@82: }
Chris@82: 
Chris@82: static void hash(const problem *p_, md5 *m)
Chris@82: {
Chris@82:      const problem_rdft *p = (const problem_rdft *) p_;
Chris@82:      X(md5puts)(m, "rdft");
Chris@82:      X(md5int)(m, p->I == p->O);
Chris@82:      kind_hash(m, p->kind, p->sz->rnk);
Chris@82:      X(md5int)(m, X(ialignment_of)(p->I));
Chris@82:      X(md5int)(m, X(ialignment_of)(p->O));
Chris@82:      X(tensor_md5)(m, p->sz);
Chris@82:      X(tensor_md5)(m, p->vecsz);
Chris@82: }
Chris@82: 
Chris@82: static void recur(const iodim *dims, int rnk, R *I)
Chris@82: {
Chris@82:      if (rnk == RNK_MINFTY)
Chris@82:           return;
Chris@82:      else if (rnk == 0)
Chris@82:           I[0] = K(0.0);
Chris@82:      else if (rnk > 0) {
Chris@82:           INT i, n = dims[0].n, is = dims[0].is;
Chris@82: 
Chris@82: 	  if (rnk == 1) {
Chris@82: 	       /* this case is redundant but faster */
Chris@82: 	       for (i = 0; i < n; ++i)
Chris@82: 		    I[i * is] = K(0.0);
Chris@82: 	  } else {
Chris@82: 	       for (i = 0; i < n; ++i)
Chris@82: 		    recur(dims + 1, rnk - 1, I + i * is);
Chris@82: 	  }
Chris@82:      }
Chris@82: }
Chris@82: 
Chris@82: void X(rdft_zerotens)(tensor *sz, R *I)
Chris@82: {
Chris@82:      recur(sz->dims, sz->rnk, I);
Chris@82: }
Chris@82: 
Chris@82: #define KSTR_LEN 8
Chris@82: 
Chris@82: const char *X(rdft_kind_str)(rdft_kind kind)
Chris@82: {
Chris@82:      static const char kstr[][KSTR_LEN] = {
Chris@82: 	  "r2hc", "r2hc01", "r2hc10", "r2hc11",
Chris@82: 	  "hc2r", "hc2r01", "hc2r10", "hc2r11",
Chris@82: 	  "dht",
Chris@82: 	  "redft00", "redft01", "redft10", "redft11",
Chris@82: 	  "rodft00", "rodft01", "rodft10", "rodft11"
Chris@82:      };
Chris@82:      A(kind >= 0 && kind < sizeof(kstr) / KSTR_LEN);
Chris@82:      return kstr[kind];
Chris@82: }
Chris@82: 
Chris@82: static void print(const problem *ego_, printer *p)
Chris@82: {
Chris@82:      const problem_rdft *ego = (const problem_rdft *) ego_;
Chris@82:      int i;
Chris@82:      p->print(p, "(rdft %d %D %T %T", 
Chris@82: 	      X(ialignment_of)(ego->I),
Chris@82: 	      (INT)(ego->O - ego->I), 
Chris@82: 	      ego->sz,
Chris@82: 	      ego->vecsz);
Chris@82:      for (i = 0; i < ego->sz->rnk; ++i)
Chris@82: 	  p->print(p, " %d", (int)ego->kind[i]);
Chris@82:      p->print(p, ")");
Chris@82: }
Chris@82: 
Chris@82: static void zero(const problem *ego_)
Chris@82: {
Chris@82:      const problem_rdft *ego = (const problem_rdft *) ego_;
Chris@82:      tensor *sz = X(tensor_append)(ego->vecsz, ego->sz);
Chris@82:      X(rdft_zerotens)(sz, UNTAINT(ego->I));
Chris@82:      X(tensor_destroy)(sz);
Chris@82: }
Chris@82: 
Chris@82: static const problem_adt padt =
Chris@82: {
Chris@82:      PROBLEM_RDFT,
Chris@82:      hash,
Chris@82:      zero,
Chris@82:      print,
Chris@82:      destroy
Chris@82: };
Chris@82: 
Chris@82: /* Dimensions of size 1 that are not REDFT/RODFT are no-ops and can be
Chris@82:    eliminated.  REDFT/RODFT unit dimensions often have factors of 2.0
Chris@82:    and suchlike from normalization and phases, although in principle
Chris@82:    these constant factors from different dimensions could be combined. */
Chris@82: static int nontrivial(const iodim *d, rdft_kind kind)
Chris@82: {
Chris@82:      return (d->n > 1 || kind == R2HC11 || kind == HC2R11
Chris@82: 	     || (REODFT_KINDP(kind) && kind != REDFT01 && kind != RODFT01));
Chris@82: }
Chris@82: 
Chris@82: problem *X(mkproblem_rdft)(const tensor *sz, const tensor *vecsz,
Chris@82: 			   R *I, R *O, const rdft_kind *kind)
Chris@82: {
Chris@82:      problem_rdft *ego;
Chris@82:      int rnk = sz->rnk;
Chris@82:      int i;
Chris@82: 
Chris@82:      A(X(tensor_kosherp)(sz));
Chris@82:      A(X(tensor_kosherp)(vecsz));
Chris@82:      A(FINITE_RNK(sz->rnk));
Chris@82: 
Chris@82:      if (UNTAINT(I) == UNTAINT(O))
Chris@82: 	  I = O = JOIN_TAINT(I, O);
Chris@82: 
Chris@82:      if (I == O && !X(tensor_inplace_locations)(sz, vecsz))
Chris@82: 	  return X(mkproblem_unsolvable)();
Chris@82: 
Chris@82:      for (i = rnk = 0; i < sz->rnk; ++i) {
Chris@82:           A(sz->dims[i].n > 0);
Chris@82:           if (nontrivial(sz->dims + i, kind[i]))
Chris@82:                ++rnk;
Chris@82:      }
Chris@82: 
Chris@82: #if defined(STRUCT_HACK_KR)
Chris@82:      ego = (problem_rdft *) X(mkproblem)(sizeof(problem_rdft)
Chris@82: 					 + sizeof(rdft_kind)
Chris@82: 					 * (rnk > 0 ? rnk - 1u : 0u), &padt);
Chris@82: #elif defined(STRUCT_HACK_C99)
Chris@82:      ego = (problem_rdft *) X(mkproblem)(sizeof(problem_rdft)
Chris@82: 					 + sizeof(rdft_kind) * (unsigned)rnk, &padt);
Chris@82: #else
Chris@82:      ego = (problem_rdft *) X(mkproblem)(sizeof(problem_rdft), &padt);
Chris@82:      ego->kind = (rdft_kind *) MALLOC(sizeof(rdft_kind) * (unsigned)rnk, PROBLEMS);
Chris@82: #endif
Chris@82: 
Chris@82:      /* do compression and sorting as in X(tensor_compress), but take
Chris@82: 	transform kind into account (sigh) */
Chris@82:      ego->sz = X(mktensor)(rnk);
Chris@82:      for (i = rnk = 0; i < sz->rnk; ++i) {
Chris@82:           if (nontrivial(sz->dims + i, kind[i])) {
Chris@82: 	       ego->kind[rnk] = kind[i];
Chris@82:                ego->sz->dims[rnk++] = sz->dims[i];
Chris@82: 	  }
Chris@82:      }
Chris@82:      for (i = 0; i + 1 < rnk; ++i) {
Chris@82: 	  int j;
Chris@82: 	  for (j = i + 1; j < rnk; ++j)
Chris@82: 	       if (X(dimcmp)(ego->sz->dims + i, ego->sz->dims + j) > 0) {
Chris@82: 		    iodim dswap;
Chris@82: 		    rdft_kind kswap;
Chris@82: 		    dswap = ego->sz->dims[i];
Chris@82: 		    ego->sz->dims[i] = ego->sz->dims[j];
Chris@82: 		    ego->sz->dims[j] = dswap;
Chris@82: 		    kswap = ego->kind[i];
Chris@82: 		    ego->kind[i] = ego->kind[j];
Chris@82: 		    ego->kind[j] = kswap;
Chris@82: 	       }
Chris@82:      }
Chris@82: 
Chris@82:      for (i = 0; i < rnk; ++i)
Chris@82: 	  if (ego->sz->dims[i].n == 2 && (ego->kind[i] == REDFT00
Chris@82: 					  || ego->kind[i] == DHT
Chris@82: 					  || ego->kind[i] == HC2R))
Chris@82: 	       ego->kind[i] = R2HC; /* size-2 transforms are equivalent */
Chris@82: 
Chris@82:      ego->vecsz = X(tensor_compress_contiguous)(vecsz);
Chris@82:      ego->I = I;
Chris@82:      ego->O = O;
Chris@82: 
Chris@82:      A(FINITE_RNK(ego->sz->rnk));
Chris@82: 
Chris@82:      return &(ego->super);
Chris@82: }
Chris@82: 
Chris@82: /* Same as X(mkproblem_rdft), but also destroy input tensors. */
Chris@82: problem *X(mkproblem_rdft_d)(tensor *sz, tensor *vecsz,
Chris@82: 			     R *I, R *O, const rdft_kind *kind)
Chris@82: {
Chris@82:      problem *p = X(mkproblem_rdft)(sz, vecsz, I, O, kind);
Chris@82:      X(tensor_destroy2)(vecsz, sz);
Chris@82:      return p;
Chris@82: }
Chris@82: 
Chris@82: /* As above, but for rnk <= 1 only and takes a scalar kind parameter */
Chris@82: problem *X(mkproblem_rdft_1)(const tensor *sz, const tensor *vecsz,
Chris@82: 			     R *I, R *O, rdft_kind kind)
Chris@82: {
Chris@82:      A(sz->rnk <= 1);
Chris@82:      return X(mkproblem_rdft)(sz, vecsz, I, O, &kind);
Chris@82: }
Chris@82: 
Chris@82: problem *X(mkproblem_rdft_1_d)(tensor *sz, tensor *vecsz,
Chris@82: 			       R *I, R *O, rdft_kind kind)
Chris@82: {
Chris@82:      A(sz->rnk <= 1);
Chris@82:      return X(mkproblem_rdft_d)(sz, vecsz, I, O, &kind);
Chris@82: }
Chris@82: 
Chris@82: /* create a zero-dimensional problem */
Chris@82: problem *X(mkproblem_rdft_0_d)(tensor *vecsz, R *I, R *O)
Chris@82: {
Chris@82:      return X(mkproblem_rdft_d)(X(mktensor_0d)(), vecsz, I, O, 
Chris@82: 				(const rdft_kind *)0);
Chris@82: }