cannam@95: /*
cannam@95:  * Copyright (c) 2003, 2007-11 Matteo Frigo
cannam@95:  * Copyright (c) 2003, 2007-11 Massachusetts Institute of Technology
cannam@95:  *
cannam@95:  * This program is free software; you can redistribute it and/or modify
cannam@95:  * it under the terms of the GNU General Public License as published by
cannam@95:  * the Free Software Foundation; either version 2 of the License, or
cannam@95:  * (at your option) any later version.
cannam@95:  *
cannam@95:  * This program is distributed in the hope that it will be useful,
cannam@95:  * but WITHOUT ANY WARRANTY; without even the implied warranty of
cannam@95:  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
cannam@95:  * GNU General Public License for more details.
cannam@95:  *
cannam@95:  * You should have received a copy of the GNU General Public License
cannam@95:  * along with this program; if not, write to the Free Software
cannam@95:  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
cannam@95:  *
cannam@95:  */
cannam@95: 
cannam@95: 
cannam@95: #include "ifftw.h"
cannam@95: 
cannam@95: static int signof(INT x)
cannam@95: {
cannam@95:      if (x < 0) return -1;
cannam@95:      if (x == 0) return 0;
cannam@95:      /* if (x > 0) */ return 1;
cannam@95: }
cannam@95: 
cannam@95: /* total order among iodim's */
cannam@95: int X(dimcmp)(const iodim *a, const iodim *b)
cannam@95: {
cannam@95:      INT sai = X(iabs)(a->is), sbi = X(iabs)(b->is);
cannam@95:      INT sao = X(iabs)(a->os), sbo = X(iabs)(b->os);
cannam@95:      INT sam = X(imin)(sai, sao), sbm = X(imin)(sbi, sbo);
cannam@95: 
cannam@95:      /* in descending order of min{istride, ostride} */
cannam@95:      if (sam != sbm)
cannam@95: 	  return signof(sbm - sam);
cannam@95: 
cannam@95:      /* in case of a tie, in descending order of istride */
cannam@95:      if (sbi != sai)
cannam@95:           return signof(sbi - sai);
cannam@95: 
cannam@95:      /* in case of a tie, in descending order of ostride */
cannam@95:      if (sbo != sao)
cannam@95:           return signof(sbo - sao);
cannam@95: 
cannam@95:      /* in case of a tie, in ascending order of n */
cannam@95:      return signof(a->n - b->n);
cannam@95: }
cannam@95: 
cannam@95: static void canonicalize(tensor *x)
cannam@95: {
cannam@95:      if (x->rnk > 1) {
cannam@95: 	  qsort(x->dims, (size_t)x->rnk, sizeof(iodim),
cannam@95: 		(int (*)(const void *, const void *))X(dimcmp));
cannam@95:      }
cannam@95: }
cannam@95: 
cannam@95: static int compare_by_istride(const iodim *a, const iodim *b)
cannam@95: {
cannam@95:      INT sai = X(iabs)(a->is), sbi = X(iabs)(b->is);
cannam@95: 
cannam@95:      /* in descending order of istride */
cannam@95:      return signof(sbi - sai);
cannam@95: }
cannam@95: 
cannam@95: static tensor *really_compress(const tensor *sz)
cannam@95: {
cannam@95:      int i, rnk;
cannam@95:      tensor *x;
cannam@95: 
cannam@95:      A(FINITE_RNK(sz->rnk));
cannam@95:      for (i = rnk = 0; i < sz->rnk; ++i) {
cannam@95:           A(sz->dims[i].n > 0);
cannam@95:           if (sz->dims[i].n != 1)
cannam@95:                ++rnk;
cannam@95:      }
cannam@95: 
cannam@95:      x = X(mktensor)(rnk);
cannam@95:      for (i = rnk = 0; i < sz->rnk; ++i) {
cannam@95:           if (sz->dims[i].n != 1)
cannam@95:                x->dims[rnk++] = sz->dims[i];
cannam@95:      }
cannam@95:      return x;
cannam@95: }
cannam@95: 
cannam@95: /* Like tensor_copy, but eliminate n == 1 dimensions, which
cannam@95:    never affect any transform or transform vector.
cannam@95:  
cannam@95:    Also, we sort the tensor into a canonical order of decreasing
cannam@95:    strides (see X(dimcmp) for an exact definition).  In general,
cannam@95:    processing a loop/array in order of decreasing stride will improve
cannam@95:    locality.  Both forward and backwards traversal of the tensor are
cannam@95:    considered e.g. by vrank-geq1, so sorting in increasing
cannam@95:    vs. decreasing order is not really important. */
cannam@95: tensor *X(tensor_compress)(const tensor *sz)
cannam@95: {
cannam@95:      tensor *x = really_compress(sz);
cannam@95:      canonicalize(x);
cannam@95:      return x;
cannam@95: }
cannam@95: 
cannam@95: /* Return whether the strides of a and b are such that they form an
cannam@95:    effective contiguous 1d array.  Assumes that a.is >= b.is. */
cannam@95: static int strides_contig(iodim *a, iodim *b)
cannam@95: {
cannam@95:      return (a->is == b->is * b->n && a->os == b->os * b->n);
cannam@95: }
cannam@95: 
cannam@95: /* Like tensor_compress, but also compress into one dimension any
cannam@95:    group of dimensions that form a contiguous block of indices with
cannam@95:    some stride.  (This can safely be done for transform vector sizes.) */
cannam@95: tensor *X(tensor_compress_contiguous)(const tensor *sz)
cannam@95: {
cannam@95:      int i, rnk;
cannam@95:      tensor *sz2, *x;
cannam@95: 
cannam@95:      if (X(tensor_sz)(sz) == 0) 
cannam@95: 	  return X(mktensor)(RNK_MINFTY);
cannam@95: 
cannam@95:      sz2 = really_compress(sz);
cannam@95:      A(FINITE_RNK(sz2->rnk));
cannam@95: 
cannam@95:      if (sz2->rnk <= 1) { /* nothing to compress. */ 
cannam@95: 	  if (0) {
cannam@95: 	       /* this call is redundant, because "sz->rnk <= 1" implies
cannam@95: 		  that the tensor is already canonical, but I am writing
cannam@95: 		  it explicitly because "logically" we need to canonicalize
cannam@95: 		  the tensor before returning. */
cannam@95: 	       canonicalize(sz2);
cannam@95: 	  }
cannam@95:           return sz2;
cannam@95:      }
cannam@95: 
cannam@95:      /* sort in descending order of |istride|, so that compressible
cannam@95: 	dimensions appear contigously */
cannam@95:      qsort(sz2->dims, (size_t)sz2->rnk, sizeof(iodim),
cannam@95: 		(int (*)(const void *, const void *))compare_by_istride);
cannam@95: 
cannam@95:      /* compute what the rank will be after compression */
cannam@95:      for (i = rnk = 1; i < sz2->rnk; ++i)
cannam@95:           if (!strides_contig(sz2->dims + i - 1, sz2->dims + i))
cannam@95:                ++rnk;
cannam@95: 
cannam@95:      /* merge adjacent dimensions whenever possible */
cannam@95:      x = X(mktensor)(rnk);
cannam@95:      x->dims[0] = sz2->dims[0];
cannam@95:      for (i = rnk = 1; i < sz2->rnk; ++i) {
cannam@95:           if (strides_contig(sz2->dims + i - 1, sz2->dims + i)) {
cannam@95:                x->dims[rnk - 1].n *= sz2->dims[i].n;
cannam@95:                x->dims[rnk - 1].is = sz2->dims[i].is;
cannam@95:                x->dims[rnk - 1].os = sz2->dims[i].os;
cannam@95:           } else {
cannam@95:                A(rnk < x->rnk);
cannam@95:                x->dims[rnk++] = sz2->dims[i];
cannam@95:           }
cannam@95:      }
cannam@95: 
cannam@95:      X(tensor_destroy)(sz2);
cannam@95: 
cannam@95:      /* reduce to canonical form */
cannam@95:      canonicalize(x);
cannam@95:      return x;
cannam@95: }
cannam@95: 
cannam@95: /* The inverse of X(tensor_append): splits the sz tensor into
cannam@95:    tensor a followed by tensor b, where a's rank is arnk. */
cannam@95: void X(tensor_split)(const tensor *sz, tensor **a, int arnk, tensor **b)
cannam@95: {
cannam@95:      A(FINITE_RNK(sz->rnk) && FINITE_RNK(arnk));
cannam@95: 
cannam@95:      *a = X(tensor_copy_sub)(sz, 0, arnk);
cannam@95:      *b = X(tensor_copy_sub)(sz, arnk, sz->rnk - arnk);
cannam@95: }
cannam@95: 
cannam@95: /* TRUE if the two tensors are equal */
cannam@95: int X(tensor_equal)(const tensor *a, const tensor *b)
cannam@95: {
cannam@95:      if (a->rnk != b->rnk)
cannam@95: 	  return 0;
cannam@95: 
cannam@95:      if (FINITE_RNK(a->rnk)) {
cannam@95: 	  int i;
cannam@95: 	  for (i = 0; i < a->rnk; ++i) 
cannam@95: 	       if (0
cannam@95: 		   || a->dims[i].n != b->dims[i].n
cannam@95: 		   || a->dims[i].is != b->dims[i].is
cannam@95: 		   || a->dims[i].os != b->dims[i].os
cannam@95: 		    )
cannam@95: 		    return 0;
cannam@95:      }
cannam@95: 
cannam@95:      return 1;
cannam@95: }
cannam@95: 
cannam@95: /* TRUE if the sets of input and output locations described by
cannam@95:    (append sz vecsz) are the same */
cannam@95: int X(tensor_inplace_locations)(const tensor *sz, const tensor *vecsz)
cannam@95: {
cannam@95:      tensor *t = X(tensor_append)(sz, vecsz);
cannam@95:      tensor *ti = X(tensor_copy_inplace)(t, INPLACE_IS);
cannam@95:      tensor *to = X(tensor_copy_inplace)(t, INPLACE_OS);
cannam@95:      tensor *tic = X(tensor_compress_contiguous)(ti);
cannam@95:      tensor *toc = X(tensor_compress_contiguous)(to);
cannam@95: 
cannam@95:      int retval = X(tensor_equal)(tic, toc);
cannam@95: 
cannam@95:      X(tensor_destroy)(t);
cannam@95:      X(tensor_destroy4)(ti, to, tic, toc);
cannam@95: 
cannam@95:      return retval;
cannam@95: }