cannam@167: /*
cannam@167:  * Copyright (c) 2003, 2007-14 Matteo Frigo
cannam@167:  * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
cannam@167:  *
cannam@167:  * This program is free software; you can redistribute it and/or modify
cannam@167:  * it under the terms of the GNU General Public License as published by
cannam@167:  * the Free Software Foundation; either version 2 of the License, or
cannam@167:  * (at your option) any later version.
cannam@167:  *
cannam@167:  * This program is distributed in the hope that it will be useful,
cannam@167:  * but WITHOUT ANY WARRANTY; without even the implied warranty of
cannam@167:  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
cannam@167:  * GNU General Public License for more details.
cannam@167:  *
cannam@167:  * You should have received a copy of the GNU General Public License
cannam@167:  * along with this program; if not, write to the Free Software
cannam@167:  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
cannam@167:  *
cannam@167:  */
cannam@167: 
cannam@167: #include "kernel/ifftw.h"
cannam@167: 
cannam@167: /* in place square transposition, iterative */
cannam@167: void X(transpose)(R *I, INT n, INT s0, INT s1, INT vl)
cannam@167: {
cannam@167:      INT i0, i1, v;
cannam@167: 
cannam@167:      switch (vl) {
cannam@167: 	 case 1:
cannam@167: 	      for (i1 = 1; i1 < n; ++i1) {
cannam@167: 		   for (i0 = 0; i0 < i1; ++i0) {
cannam@167: 			R x0 = I[i1 * s0 + i0 * s1];
cannam@167: 			R y0 = I[i1 * s1 + i0 * s0];
cannam@167: 			I[i1 * s1 + i0 * s0] = x0;
cannam@167: 			I[i1 * s0 + i0 * s1] = y0;
cannam@167: 		   }
cannam@167: 	      }
cannam@167: 	      break;
cannam@167: 	 case 2:
cannam@167: 	      for (i1 = 1; i1 < n; ++i1) {
cannam@167: 		   for (i0 = 0; i0 < i1; ++i0) {
cannam@167: 			R x0 = I[i1 * s0 + i0 * s1];
cannam@167: 			R x1 = I[i1 * s0 + i0 * s1 + 1];
cannam@167: 			R y0 = I[i1 * s1 + i0 * s0];
cannam@167: 			R y1 = I[i1 * s1 + i0 * s0 + 1];
cannam@167: 			I[i1 * s1 + i0 * s0] = x0;
cannam@167: 			I[i1 * s1 + i0 * s0 + 1] = x1;
cannam@167: 			I[i1 * s0 + i0 * s1] = y0;
cannam@167: 			I[i1 * s0 + i0 * s1 + 1] = y1;
cannam@167: 		   }
cannam@167: 	      }
cannam@167: 	      break;
cannam@167: 	 default:
cannam@167: 	      for (i1 = 1; i1 < n; ++i1) {
cannam@167: 		   for (i0 = 0; i0 < i1; ++i0) {
cannam@167: 			for (v = 0; v < vl; ++v) {
cannam@167: 			     R x0 = I[i1 * s0 + i0 * s1 + v];
cannam@167: 			     R y0 = I[i1 * s1 + i0 * s0 + v];
cannam@167: 			     I[i1 * s1 + i0 * s0 + v] = x0;
cannam@167: 			     I[i1 * s0 + i0 * s1 + v] = y0;
cannam@167: 			}
cannam@167: 		   }
cannam@167: 	      }
cannam@167: 	      break;
cannam@167:      }
cannam@167: }
cannam@167: 
cannam@167: struct transpose_closure {
cannam@167:      R *I;
cannam@167:      INT s0, s1, vl, tilesz;
cannam@167:      R *buf0, *buf1; 
cannam@167: };
cannam@167: 
cannam@167: static void dotile(INT n0l, INT n0u, INT n1l, INT n1u, void *args)
cannam@167: {
cannam@167:      struct transpose_closure *k = (struct transpose_closure *)args;
cannam@167:      R *I = k->I;
cannam@167:      INT s0 = k->s0, s1 = k->s1, vl = k->vl;
cannam@167:      INT i0, i1, v;
cannam@167: 
cannam@167:      switch (vl) {
cannam@167: 	 case 1:
cannam@167: 	      for (i1 = n1l; i1 < n1u; ++i1) {
cannam@167: 		   for (i0 = n0l; i0 < n0u; ++i0) {
cannam@167: 			R x0 = I[i1 * s0 + i0 * s1];
cannam@167: 			R y0 = I[i1 * s1 + i0 * s0];
cannam@167: 			I[i1 * s1 + i0 * s0] = x0;
cannam@167: 			I[i1 * s0 + i0 * s1] = y0;
cannam@167: 		   }
cannam@167: 	      }
cannam@167: 	      break;
cannam@167: 	 case 2:
cannam@167: 	      for (i1 = n1l; i1 < n1u; ++i1) {
cannam@167: 		   for (i0 = n0l; i0 < n0u; ++i0) {
cannam@167: 			R x0 = I[i1 * s0 + i0 * s1];
cannam@167: 			R x1 = I[i1 * s0 + i0 * s1 + 1];
cannam@167: 			R y0 = I[i1 * s1 + i0 * s0];
cannam@167: 			R y1 = I[i1 * s1 + i0 * s0 + 1];
cannam@167: 			I[i1 * s1 + i0 * s0] = x0;
cannam@167: 			I[i1 * s1 + i0 * s0 + 1] = x1;
cannam@167: 			I[i1 * s0 + i0 * s1] = y0;
cannam@167: 			I[i1 * s0 + i0 * s1 + 1] = y1;
cannam@167: 		   }
cannam@167: 	      }
cannam@167: 	      break;
cannam@167: 	 default:
cannam@167: 	      for (i1 = n1l; i1 < n1u; ++i1) {
cannam@167: 		   for (i0 = n0l; i0 < n0u; ++i0) {
cannam@167: 			for (v = 0; v < vl; ++v) {
cannam@167: 			     R x0 = I[i1 * s0 + i0 * s1 + v];
cannam@167: 			     R y0 = I[i1 * s1 + i0 * s0 + v];
cannam@167: 			     I[i1 * s1 + i0 * s0 + v] = x0;
cannam@167: 			     I[i1 * s0 + i0 * s1 + v] = y0;
cannam@167: 			}
cannam@167: 		   }
cannam@167: 	      }
cannam@167:      }
cannam@167: }
cannam@167: 
cannam@167: static void dotile_buf(INT n0l, INT n0u, INT n1l, INT n1u, void *args)
cannam@167: {
cannam@167:      struct transpose_closure *k = (struct transpose_closure *)args;
cannam@167:      X(cpy2d_ci)(k->I + n0l * k->s0 + n1l * k->s1,
cannam@167: 		 k->buf0,
cannam@167: 		 n0u - n0l, k->s0, k->vl,
cannam@167: 		 n1u - n1l, k->s1, k->vl * (n0u - n0l),
cannam@167: 		 k->vl);
cannam@167:      X(cpy2d_ci)(k->I + n0l * k->s1 + n1l * k->s0,
cannam@167: 		 k->buf1,
cannam@167: 		 n0u - n0l, k->s1, k->vl,
cannam@167: 		 n1u - n1l, k->s0, k->vl * (n0u - n0l),
cannam@167: 		 k->vl);
cannam@167:      X(cpy2d_co)(k->buf1,
cannam@167: 		 k->I + n0l * k->s0 + n1l * k->s1,
cannam@167: 		 n0u - n0l, k->vl, k->s0,
cannam@167: 		 n1u - n1l, k->vl * (n0u - n0l), k->s1,
cannam@167: 		 k->vl);
cannam@167:      X(cpy2d_co)(k->buf0,
cannam@167: 		 k->I + n0l * k->s1 + n1l * k->s0,
cannam@167: 		 n0u - n0l, k->vl, k->s1,
cannam@167: 		 n1u - n1l, k->vl * (n0u - n0l), k->s0,
cannam@167: 		 k->vl);
cannam@167: }
cannam@167: 
cannam@167: static void transpose_rec(R *I, INT n,
cannam@167: 			  void (*f)(INT n0l, INT n0u, INT n1l, INT n1u,
cannam@167: 				    void *args),
cannam@167: 			  struct transpose_closure *k)
cannam@167: {
cannam@167:    tail:
cannam@167:      if (n > 1) {
cannam@167: 	  INT n2 = n / 2;
cannam@167: 	  k->I = I;
cannam@167: 	  X(tile2d)(0, n2, n2, n, k->tilesz, f, k);
cannam@167: 	  transpose_rec(I, n2, f, k);
cannam@167: 	  I += n2 * (k->s0 + k->s1); n -= n2; goto tail;
cannam@167:      }
cannam@167: }
cannam@167: 
cannam@167: void X(transpose_tiled)(R *I, INT n, INT s0, INT s1, INT vl) 
cannam@167: {
cannam@167:      struct transpose_closure k;
cannam@167:      k.s0 = s0;
cannam@167:      k.s1 = s1;
cannam@167:      k.vl = vl;
cannam@167:      /* two blocks must be in cache, to be swapped */
cannam@167:      k.tilesz = X(compute_tilesz)(vl, 2);
cannam@167:      k.buf0 = k.buf1 = 0; /* unused */
cannam@167:      transpose_rec(I, n, dotile, &k);
cannam@167: }
cannam@167: 
cannam@167: void X(transpose_tiledbuf)(R *I, INT n, INT s0, INT s1, INT vl) 
cannam@167: {
cannam@167:      struct transpose_closure k;
cannam@167:      /* Assume that the the rows of I conflict into the same cache
cannam@167:         lines, and therefore we don't need to reserve cache space for
cannam@167:         the input.  If the rows don't conflict, there is no reason
cannam@167: 	to use tiledbuf at all.*/
cannam@167:      R buf0[CACHESIZE / (2 * sizeof(R))];
cannam@167:      R buf1[CACHESIZE / (2 * sizeof(R))];
cannam@167:      k.s0 = s0;
cannam@167:      k.s1 = s1;
cannam@167:      k.vl = vl;
cannam@167:      k.tilesz = X(compute_tilesz)(vl, 2);
cannam@167:      k.buf0 = buf0;
cannam@167:      k.buf1 = buf1;
cannam@167:      A(k.tilesz * k.tilesz * vl * sizeof(R) <= sizeof(buf0));
cannam@167:      A(k.tilesz * k.tilesz * vl * sizeof(R) <= sizeof(buf1));
cannam@167:      transpose_rec(I, n, dotile_buf, &k);
cannam@167: }
cannam@167: