2
* Copyright (c) 2003, 2006 Matteo Frigo
3
* Copyright (c) 2003, 2006 Massachusetts Institute of Technology
5
* This program is free software; you can redistribute it and/or modify
6
* it under the terms of the GNU General Public License as published by
7
* the Free Software Foundation; either version 2 of the License, or
8
* (at your option) any later version.
10
* This program is distributed in the hope that it will be useful,
11
* but WITHOUT ANY WARRANTY; without even the implied warranty of
12
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13
* GNU General Public License for more details.
15
* You should have received a copy of the GNU General Public License
16
* along with this program; if not, write to the Free Software
17
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21
/* $Id: rank0-rdft2.c,v 1.13 2006-01-27 02:10:50 athena Exp $ */
23
/* plans for rank-0 RDFT2 (copy operations, plus setting 0 imag. parts) */
28
#include <string.h> /* for memcpy() */
42
static int applicable(const problem *p_)
44
const problem_rdft2 *p = (const problem_rdft2 *) p_;
48
|| (((p->r != p->rio && p->r != p->iio)
49
|| X(rdft2_inplace_strides)(p, RNK_MINFTY))
50
&& p->vecsz->rnk <= 1))
54
static void apply_r2hc(const plan *ego_, R *r, R *rio, R *iio)
56
const P *ego = (const P *) ego_;
58
INT ivs = ego->ivs, ovs = ego->ovs;
60
for (i = 4; i <= vl; i += 4) {
66
*rio = r0; rio += ovs;
67
*iio = K(0.0); iio += ovs;
68
*rio = r1; rio += ovs;
69
*iio = K(0.0); iio += ovs;
70
*rio = r2; rio += ovs;
71
*iio = K(0.0); iio += ovs;
72
*rio = r3; rio += ovs;
73
*iio = K(0.0); iio += ovs;
75
for (; i < vl + 4; ++i) {
78
*rio = r0; rio += ovs;
79
*iio = K(0.0); iio += ovs;
83
/* in-place r2hc rank-0: set imaginary parts of output to 0 */
84
static void apply_r2hc_inplace(const plan *ego_, R *r, R *rio, R *iio)
86
const P *ego = (const P *) ego_;
92
for (i = 4; i <= vl; i += 4) {
93
*iio = K(0.0); iio += ovs;
94
*iio = K(0.0); iio += ovs;
95
*iio = K(0.0); iio += ovs;
96
*iio = K(0.0); iio += ovs;
98
for (; i < vl + 4; ++i) {
99
*iio = K(0.0); iio += ovs;
103
/* a rank-0 HC2R rdft2 problem is just a copy from rio to r,
104
so we can use a rank-0 rdft plan */
105
static void apply_hc2r(const plan *ego_, R *r, R *rio, R *iio)
107
const P *ego = (const P *) ego_;
108
plan_rdft *cldcpy = (plan_rdft *) ego->cldcpy;
110
cldcpy->apply((plan *) cldcpy, rio, r);
113
static void awake(plan *ego_, enum wakefulness wakefulness)
117
X(plan_awake)(ego->cldcpy, wakefulness);
120
static void destroy(plan *ego_)
124
X(plan_destroy_internal)(ego->cldcpy);
127
static void print(const plan *ego_, printer *p)
129
const P *ego = (const P *) ego_;
131
p->print(p, "(rdft2-hc2r-rank0%(%p%))", ego->cldcpy);
133
p->print(p, "(rdft2-r2hc-rank0%v)", ego->vl);
136
static plan *mkplan(const solver *ego_, const problem *p_, planner *plnr)
138
const problem_rdft2 *p;
139
plan *cldcpy = (plan *) 0;
142
static const plan_adt padt = {
143
X(rdft2_solve), awake, print, destroy
151
p = (const problem_rdft2 *) p_;
153
if (p->kind == HC2R) {
154
cldcpy = X(mkplan_d)(plnr,
155
X(mkproblem_rdft_0_d)(
156
X(tensor_copy)(p->vecsz),
158
if (!cldcpy) return (plan *) 0;
161
pln = MKPLAN_RDFT2(P, &padt,
163
(p->r == p->rio ? apply_r2hc_inplace : apply_r2hc)
167
X(tensor_tornk1)(p->vecsz, &pln->vl, &pln->ivs, &pln->ovs);
168
pln->cldcpy = cldcpy;
170
if (p->kind == R2HC) {
171
/* vl loads, 2*vl stores */
172
X(ops_other)(3 * pln->vl, &pln->super.super.ops);
175
pln->super.super.ops = cldcpy->ops;
178
return &(pln->super.super);
181
static solver *mksolver(void)
183
static const solver_adt sadt = { PROBLEM_RDFT2, mkplan };
184
S *slv = MKSOLVER(S, &sadt);
185
return &(slv->super);
188
void X(rdft2_rank0_register)(planner *p)
190
REGISTER_SOLVER(p, mksolver());