1 /*
2  *  This file is part of libsharp.
3  *
4  *  libsharp is free software; you can redistribute it and/or modify
5  *  it under the terms of the GNU General Public License as published by
6  *  the Free Software Foundation; either version 2 of the License, or
7  *  (at your option) any later version.
8  *
9  *  libsharp is distributed in the hope that it will be useful,
10  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
11  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  *  GNU General Public License for more details.
13  *
14  *  You should have received a copy of the GNU General Public License
15  *  along with libsharp; if not, write to the Free Software
16  *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
17  */
18 
19 /*
20  *  libsharp is being developed at the Max-Planck-Institut fuer Astrophysik
21  *  and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
22  *  (DLR).
23  */
24 
25 /*! \file sharp_core.c
26  *  Computational core
27  *
28  *  Copyright (C) 2012 Max-Planck-Society
29  *  \author Martin Reinecke
30  */
31 
32 #include <complex.h>
33 #include <math.h>
34 #include <string.h>
35 #include "sharp_vecsupport.h"
36 #include "sharp_complex_hacks.h"
37 #include "sharp_ylmgen_c.h"
38 #include "sharp.h"
39 #include "sharp_core.h"
40 #include "c_utils.h"
41 
42 typedef complex double dcmplx;
43 
44 // must be in the range [0;6]
45 #define MAXJOB_SPECIAL 2
46 
47 #define XCONCAT2(a,b) a##_##b
48 #define CONCAT2(a,b) XCONCAT2(a,b)
49 #define XCONCAT3(a,b,c) a##_##b##_##c
50 #define CONCAT3(a,b,c) XCONCAT3(a,b,c)
51 
52 #define nvec 1
53 #include "sharp_inchelper1.inc.c"
54 #undef nvec
55 
56 #define nvec 2
57 #include "sharp_inchelper1.inc.c"
58 #undef nvec
59 
60 #define nvec 3
61 #include "sharp_inchelper1.inc.c"
62 #undef nvec
63 
64 #define nvec 4
65 #include "sharp_inchelper1.inc.c"
66 #undef nvec
67 
68 #define nvec 5
69 #include "sharp_inchelper1.inc.c"
70 #undef nvec
71 
72 #define nvec 6
73 #include "sharp_inchelper1.inc.c"
74 #undef nvec
75 
inner_loop(sharp_job * job,const int * ispair,const double * cth,const double * sth,int llim,int ulim,sharp_Ylmgen_C * gen,int mi,const int * idx)76 void inner_loop (sharp_job *job, const int *ispair,const double *cth,
77   const double *sth, int llim, int ulim, sharp_Ylmgen_C *gen, int mi,
78   const int *idx)
79   {
80   int njobs=job->ntrans, nv=job->flags&SHARP_NVMAX;
81   if (njobs<=MAXJOB_SPECIAL)
82     {
83     switch (njobs*16+nv)
84       {
85 #if ((MAXJOB_SPECIAL>=1)&&(SHARP_MAXTRANS>=1))
86       case 0x11:
87         CONCAT3(inner_loop,1,1) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
88         return;
89       case 0x12:
90         CONCAT3(inner_loop,2,1) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
91         return;
92       case 0x13:
93         CONCAT3(inner_loop,3,1) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
94         return;
95       case 0x14:
96         CONCAT3(inner_loop,4,1) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
97         return;
98       case 0x15:
99         CONCAT3(inner_loop,5,1) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
100         return;
101       case 0x16:
102         CONCAT3(inner_loop,6,1) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
103         return;
104 #endif
105 #if ((MAXJOB_SPECIAL>=2)&&(SHARP_MAXTRANS>=2))
106       case 0x21:
107         CONCAT3(inner_loop,1,2) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
108         return;
109       case 0x22:
110         CONCAT3(inner_loop,2,2) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
111         return;
112       case 0x23:
113         CONCAT3(inner_loop,3,2) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
114         return;
115       case 0x24:
116         CONCAT3(inner_loop,4,2) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
117         return;
118       case 0x25:
119         CONCAT3(inner_loop,5,2) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
120         return;
121       case 0x26:
122         CONCAT3(inner_loop,6,2) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
123         return;
124 #endif
125 #if ((MAXJOB_SPECIAL>=3)&&(SHARP_MAXTRANS>=3))
126       case 0x31:
127         CONCAT3(inner_loop,1,3) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
128         return;
129       case 0x32:
130         CONCAT3(inner_loop,2,3) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
131         return;
132       case 0x33:
133         CONCAT3(inner_loop,3,3) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
134         return;
135       case 0x34:
136         CONCAT3(inner_loop,4,3) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
137         return;
138       case 0x35:
139         CONCAT3(inner_loop,5,3) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
140         return;
141       case 0x36:
142         CONCAT3(inner_loop,6,3) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
143         return;
144 #endif
145 #if ((MAXJOB_SPECIAL>=4)&&(SHARP_MAXTRANS>=4))
146       case 0x41:
147         CONCAT3(inner_loop,1,4) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
148         return;
149       case 0x42:
150         CONCAT3(inner_loop,2,4) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
151         return;
152       case 0x43:
153         CONCAT3(inner_loop,3,4) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
154         return;
155       case 0x44:
156         CONCAT3(inner_loop,4,4) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
157         return;
158       case 0x45:
159         CONCAT3(inner_loop,5,4) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
160         return;
161       case 0x46:
162         CONCAT3(inner_loop,6,4) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
163         return;
164 #endif
165 #if ((MAXJOB_SPECIAL>=5)&&(SHARP_MAXTRANS>=5))
166       case 0x51:
167         CONCAT3(inner_loop,1,5) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
168         return;
169       case 0x52:
170         CONCAT3(inner_loop,2,5) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
171         return;
172       case 0x53:
173         CONCAT3(inner_loop,3,5) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
174         return;
175       case 0x54:
176         CONCAT3(inner_loop,4,5) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
177         return;
178       case 0x55:
179         CONCAT3(inner_loop,5,5) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
180         return;
181       case 0x56:
182         CONCAT3(inner_loop,6,5) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
183         return;
184 #endif
185 #if ((MAXJOB_SPECIAL>=6)&&(SHARP_MAXTRANS>=6))
186       case 0x61:
187         CONCAT3(inner_loop,1,6) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
188         return;
189       case 0x62:
190         CONCAT3(inner_loop,2,6) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
191         return;
192       case 0x63:
193         CONCAT3(inner_loop,3,6) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
194         return;
195       case 0x64:
196         CONCAT3(inner_loop,4,6) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
197         return;
198       case 0x65:
199         CONCAT3(inner_loop,5,6) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
200         return;
201       case 0x66:
202         CONCAT3(inner_loop,6,6) (job, ispair,cth,sth,llim,ulim,gen,mi,idx);
203         return;
204 #endif
205       }
206     }
207 #if (SHARP_MAXTRANS>MAXJOB_SPECIAL)
208   else
209     {
210     switch (nv)
211       {
212       case 1:
213         CONCAT2(inner_loop,1)
214           (job, ispair,cth,sth,llim,ulim,gen,mi,idx,job->ntrans);
215         return;
216       case 2:
217         CONCAT2(inner_loop,2)
218           (job, ispair,cth,sth,llim,ulim,gen,mi,idx,job->ntrans);
219         return;
220       case 3:
221         CONCAT2(inner_loop,3)
222           (job, ispair,cth,sth,llim,ulim,gen,mi,idx,job->ntrans);
223         return;
224       case 4:
225         CONCAT2(inner_loop,4)
226           (job, ispair,cth,sth,llim,ulim,gen,mi,idx,job->ntrans);
227         return;
228       case 5:
229         CONCAT2(inner_loop,5)
230           (job, ispair,cth,sth,llim,ulim,gen,mi,idx,job->ntrans);
231         return;
232       case 6:
233         CONCAT2(inner_loop,6)
234           (job, ispair,cth,sth,llim,ulim,gen,mi,idx,job->ntrans);
235         return;
236       }
237     }
238 #endif
239   UTIL_FAIL("Incorrect vector parameters");
240   }
241