Mercurial > hg > sv-dependency-builds

annotate src/fftw-3.3.8/kernel/transpose.c @ 82:d0c2a83c1364

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add FFTW 3.3.8 source, and a Linux build
author Chris Cannam
date Tue, 19 Nov 2019 14:52:55 +0000
parents
children
rev   line source
Chris@82 1 /*
Chris@82 2 * Copyright (c) 2003, 2007-14 Matteo Frigo
Chris@82 3 * Copyright (c) 2003, 2007-14 Massachusetts Institute of Technology
Chris@82 4 *
Chris@82 5 * This program is free software; you can redistribute it and/or modify
Chris@82 6 * it under the terms of the GNU General Public License as published by
Chris@82 7 * the Free Software Foundation; either version 2 of the License, or
Chris@82 8 * (at your option) any later version.
Chris@82 9 *
Chris@82 10 * This program is distributed in the hope that it will be useful,
Chris@82 11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
Chris@82 12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
Chris@82 13 * GNU General Public License for more details.
Chris@82 14 *
Chris@82 15 * You should have received a copy of the GNU General Public License
Chris@82 16 * along with this program; if not, write to the Free Software
Chris@82 17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
Chris@82 18 *
Chris@82 19 */
Chris@82 20
Chris@82 21 #include "kernel/ifftw.h"
Chris@82 22
Chris@82 23 /* in place square transposition, iterative */
Chris@82 24 void X(transpose)(R *I, INT n, INT s0, INT s1, INT vl)
Chris@82 25 {
Chris@82 26 INT i0, i1, v;
Chris@82 27
Chris@82 28 switch (vl) {
Chris@82 29 case 1:
Chris@82 30 for (i1 = 1; i1 < n; ++i1) {
Chris@82 31 for (i0 = 0; i0 < i1; ++i0) {
Chris@82 32 R x0 = I[i1 * s0 + i0 * s1];
Chris@82 33 R y0 = I[i1 * s1 + i0 * s0];
Chris@82 34 I[i1 * s1 + i0 * s0] = x0;
Chris@82 35 I[i1 * s0 + i0 * s1] = y0;
Chris@82 36 }
Chris@82 37 }
Chris@82 38 break;
Chris@82 39 case 2:
Chris@82 40 for (i1 = 1; i1 < n; ++i1) {
Chris@82 41 for (i0 = 0; i0 < i1; ++i0) {
Chris@82 42 R x0 = I[i1 * s0 + i0 * s1];
Chris@82 43 R x1 = I[i1 * s0 + i0 * s1 + 1];
Chris@82 44 R y0 = I[i1 * s1 + i0 * s0];
Chris@82 45 R y1 = I[i1 * s1 + i0 * s0 + 1];
Chris@82 46 I[i1 * s1 + i0 * s0] = x0;
Chris@82 47 I[i1 * s1 + i0 * s0 + 1] = x1;
Chris@82 48 I[i1 * s0 + i0 * s1] = y0;
Chris@82 49 I[i1 * s0 + i0 * s1 + 1] = y1;
Chris@82 50 }
Chris@82 51 }
Chris@82 52 break;
Chris@82 53 default:
Chris@82 54 for (i1 = 1; i1 < n; ++i1) {
Chris@82 55 for (i0 = 0; i0 < i1; ++i0) {
Chris@82 56 for (v = 0; v < vl; ++v) {
Chris@82 57 R x0 = I[i1 * s0 + i0 * s1 + v];
Chris@82 58 R y0 = I[i1 * s1 + i0 * s0 + v];
Chris@82 59 I[i1 * s1 + i0 * s0 + v] = x0;
Chris@82 60 I[i1 * s0 + i0 * s1 + v] = y0;
Chris@82 61 }
Chris@82 62 }
Chris@82 63 }
Chris@82 64 break;
Chris@82 65 }
Chris@82 66 }
Chris@82 67
Chris@82 68 struct transpose_closure {
Chris@82 69 R *I;
Chris@82 70 INT s0, s1, vl, tilesz;
Chris@82 71 R *buf0, *buf1;
Chris@82 72 };
Chris@82 73
Chris@82 74 static void dotile(INT n0l, INT n0u, INT n1l, INT n1u, void *args)
Chris@82 75 {
Chris@82 76 struct transpose_closure *k = (struct transpose_closure *)args;
Chris@82 77 R *I = k->I;
Chris@82 78 INT s0 = k->s0, s1 = k->s1, vl = k->vl;
Chris@82 79 INT i0, i1, v;
Chris@82 80
Chris@82 81 switch (vl) {
Chris@82 82 case 1:
Chris@82 83 for (i1 = n1l; i1 < n1u; ++i1) {
Chris@82 84 for (i0 = n0l; i0 < n0u; ++i0) {
Chris@82 85 R x0 = I[i1 * s0 + i0 * s1];
Chris@82 86 R y0 = I[i1 * s1 + i0 * s0];
Chris@82 87 I[i1 * s1 + i0 * s0] = x0;
Chris@82 88 I[i1 * s0 + i0 * s1] = y0;
Chris@82 89 }
Chris@82 90 }
Chris@82 91 break;
Chris@82 92 case 2:
Chris@82 93 for (i1 = n1l; i1 < n1u; ++i1) {
Chris@82 94 for (i0 = n0l; i0 < n0u; ++i0) {
Chris@82 95 R x0 = I[i1 * s0 + i0 * s1];
Chris@82 96 R x1 = I[i1 * s0 + i0 * s1 + 1];
Chris@82 97 R y0 = I[i1 * s1 + i0 * s0];
Chris@82 98 R y1 = I[i1 * s1 + i0 * s0 + 1];
Chris@82 99 I[i1 * s1 + i0 * s0] = x0;
Chris@82 100 I[i1 * s1 + i0 * s0 + 1] = x1;
Chris@82 101 I[i1 * s0 + i0 * s1] = y0;
Chris@82 102 I[i1 * s0 + i0 * s1 + 1] = y1;
Chris@82 103 }
Chris@82 104 }
Chris@82 105 break;
Chris@82 106 default:
Chris@82 107 for (i1 = n1l; i1 < n1u; ++i1) {
Chris@82 108 for (i0 = n0l; i0 < n0u; ++i0) {
Chris@82 109 for (v = 0; v < vl; ++v) {
Chris@82 110 R x0 = I[i1 * s0 + i0 * s1 + v];
Chris@82 111 R y0 = I[i1 * s1 + i0 * s0 + v];
Chris@82 112 I[i1 * s1 + i0 * s0 + v] = x0;
Chris@82 113 I[i1 * s0 + i0 * s1 + v] = y0;
Chris@82 114 }
Chris@82 115 }
Chris@82 116 }
Chris@82 117 }
Chris@82 118 }
Chris@82 119
Chris@82 120 static void dotile_buf(INT n0l, INT n0u, INT n1l, INT n1u, void *args)
Chris@82 121 {
Chris@82 122 struct transpose_closure *k = (struct transpose_closure *)args;
Chris@82 123 X(cpy2d_ci)(k->I + n0l * k->s0 + n1l * k->s1,
Chris@82 124 k->buf0,
Chris@82 125 n0u - n0l, k->s0, k->vl,
Chris@82 126 n1u - n1l, k->s1, k->vl * (n0u - n0l),
Chris@82 127 k->vl);
Chris@82 128 X(cpy2d_ci)(k->I + n0l * k->s1 + n1l * k->s0,
Chris@82 129 k->buf1,
Chris@82 130 n0u - n0l, k->s1, k->vl,
Chris@82 131 n1u - n1l, k->s0, k->vl * (n0u - n0l),
Chris@82 132 k->vl);
Chris@82 133 X(cpy2d_co)(k->buf1,
Chris@82 134 k->I + n0l * k->s0 + n1l * k->s1,
Chris@82 135 n0u - n0l, k->vl, k->s0,
Chris@82 136 n1u - n1l, k->vl * (n0u - n0l), k->s1,
Chris@82 137 k->vl);
Chris@82 138 X(cpy2d_co)(k->buf0,
Chris@82 139 k->I + n0l * k->s1 + n1l * k->s0,
Chris@82 140 n0u - n0l, k->vl, k->s1,
Chris@82 141 n1u - n1l, k->vl * (n0u - n0l), k->s0,
Chris@82 142 k->vl);
Chris@82 143 }
Chris@82 144
Chris@82 145 static void transpose_rec(R *I, INT n,
Chris@82 146 void (*f)(INT n0l, INT n0u, INT n1l, INT n1u,
Chris@82 147 void *args),
Chris@82 148 struct transpose_closure *k)
Chris@82 149 {
Chris@82 150 tail:
Chris@82 151 if (n > 1) {
Chris@82 152 INT n2 = n / 2;
Chris@82 153 k->I = I;
Chris@82 154 X(tile2d)(0, n2, n2, n, k->tilesz, f, k);
Chris@82 155 transpose_rec(I, n2, f, k);
Chris@82 156 I += n2 * (k->s0 + k->s1); n -= n2; goto tail;
Chris@82 157 }
Chris@82 158 }
Chris@82 159
Chris@82 160 void X(transpose_tiled)(R *I, INT n, INT s0, INT s1, INT vl)
Chris@82 161 {
Chris@82 162 struct transpose_closure k;
Chris@82 163 k.s0 = s0;
Chris@82 164 k.s1 = s1;
Chris@82 165 k.vl = vl;
Chris@82 166 /* two blocks must be in cache, to be swapped */
Chris@82 167 k.tilesz = X(compute_tilesz)(vl, 2);
Chris@82 168 k.buf0 = k.buf1 = 0; /* unused */
Chris@82 169 transpose_rec(I, n, dotile, &k);
Chris@82 170 }
Chris@82 171
Chris@82 172 void X(transpose_tiledbuf)(R *I, INT n, INT s0, INT s1, INT vl)
Chris@82 173 {
Chris@82 174 struct transpose_closure k;
Chris@82 175 /* Assume that the the rows of I conflict into the same cache
Chris@82 176 lines, and therefore we don't need to reserve cache space for
Chris@82 177 the input. If the rows don't conflict, there is no reason
Chris@82 178 to use tiledbuf at all.*/
Chris@82 179 R buf0[CACHESIZE / (2 * sizeof(R))];
Chris@82 180 R buf1[CACHESIZE / (2 * sizeof(R))];
Chris@82 181 k.s0 = s0;
Chris@82 182 k.s1 = s1;
Chris@82 183 k.vl = vl;
Chris@82 184 k.tilesz = X(compute_tilesz)(vl, 2);
Chris@82 185 k.buf0 = buf0;
Chris@82 186 k.buf1 = buf1;
Chris@82 187 A(k.tilesz * k.tilesz * vl * sizeof(R) <= sizeof(buf0));
Chris@82 188 A(k.tilesz * k.tilesz * vl * sizeof(R) <= sizeof(buf1));
Chris@82 189 transpose_rec(I, n, dotile_buf, &k);
Chris@82 190 }
Chris@82 191