Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Martin Reinecke
ducc
Commits
aaea79a7
Commit
aaea79a7
authored
Jan 15, 2020
by
Martin Reinecke
Browse files
introduce abstract alm_info interface
parent
ce3cb52f
Changes
7
Hide whitespace changes
Inline
Side-by-side
libsharp2/sharp.cc
View file @
aaea79a7
...
...
@@ -197,29 +197,29 @@ struct ringhelper
}
};
sharp_alm_info
::
sharp_alm_info
(
size_t
lmax_
,
size_t
nm_
,
ptrdiff_t
stride_
,
const
size_t
*
mval_
,
const
ptrdiff_t
*
mstart
)
:
lmax
(
lmax_
),
nm
(
nm
_
),
mval
(
nm_
),
mvstart
(
nm
),
stride
(
stride_
)
sharp_
standard_
alm_info
::
sharp_
standard_
alm_info
(
size_t
lmax_
_
,
size_t
nm_
,
ptrdiff_t
stride_
,
const
size_t
*
mval_
_
,
const
ptrdiff_t
*
mstart
)
:
lmax
_
(
lmax__
),
mval
_
(
nm_
),
mvstart
(
nm
_
),
stride
(
stride_
)
{
for
(
size_t
mi
=
0
;
mi
<
nm
;
++
mi
)
for
(
size_t
mi
=
0
;
mi
<
nm
_
;
++
mi
)
{
mval
[
mi
]
=
mval_
[
mi
];
mval
_
[
mi
]
=
mval_
_
[
mi
];
mvstart
[
mi
]
=
mstart
[
mi
];
}
}
sharp_alm_info
::
sharp_alm_info
(
size_t
lmax_
,
size_t
mmax
,
ptrdiff_t
stride_
,
sharp_
standard_
alm_info
::
sharp_
standard_
alm_info
(
size_t
lmax_
_
,
size_t
mmax
_
,
ptrdiff_t
stride_
,
const
ptrdiff_t
*
mstart
)
:
lmax
(
lmax_
),
nm
(
mmax
+
1
),
mval
(
mmax
+
1
),
mvstart
(
mmax
+
1
),
stride
(
stride_
)
:
lmax
_
(
lmax_
_
),
mval
_
(
mmax
_
+
1
),
mvstart
(
mmax
_
+
1
),
stride
(
stride_
)
{
for
(
size_t
i
=
0
;
i
<=
mmax
;
++
i
)
for
(
size_t
i
=
0
;
i
<=
mmax
_
;
++
i
)
{
mval
[
i
]
=
i
;
mval
_
[
i
]
=
i
;
mvstart
[
i
]
=
mstart
[
i
];
}
}
ptrdiff_t
sharp_alm_info
::
index
(
int
l
,
int
mi
)
ptrdiff_t
sharp_
standard_
alm_info
::
index
(
int
l
,
int
mi
)
{
return
mvstart
[
mi
]
+
stride
*
l
;
}
...
...
@@ -278,18 +278,18 @@ sharp_standard_geom_info::sharp_standard_geom_info(size_t nrings, const size_t *
/* This currently requires all m values from 0 to nm-1 to be present.
It might be worthwhile to relax this criterion such that holes in the m
distribution are permissible. */
static
size_t
sharp_
get_
mmax
(
const
vector
<
size_t
>
&
mval
)
size_t
sharp_
standard_alm_info
::
mmax
(
)
const
{
//FIXME: if gaps are allowed, we have to search the maximum m in the array
auto
nm
=
mval
.
size
();
vector
<
bool
>
mcheck
(
nm
,
false
);
for
(
auto
m_cur
:
mval
)
auto
nm
_
=
mval
_
.
size
();
vector
<
bool
>
mcheck
(
nm
_
,
false
);
for
(
auto
m_cur
:
mval
_
)
{
MR_assert
(
m_cur
<
nm
,
"not all m values are present"
);
MR_assert
(
m_cur
<
nm
_
,
"not all m values are present"
);
MR_assert
(
mcheck
[
m_cur
]
==
false
,
"duplicate m value"
);
mcheck
[
m_cur
]
=
true
;
}
return
nm
-
1
;
return
nm
_
-
1
;
}
MRUTIL_NOINLINE
void
sharp_standard_geom_info
::
clear_map
(
double
*
map
)
const
...
...
@@ -315,21 +315,17 @@ MRUTIL_NOINLINE void sharp_standard_geom_info::clear_map (float *map) const
}
}
MRUTIL_NOINLINE
static
void
clear_alm
(
const
sharp_alm_info
*
ainfo
,
void
*
alm
,
int
flags
)
void
sharp_standard_alm_info
::
clear_alm
(
dcmplx
*
alm
)
const
{
for
(
size_t
mi
=
0
;
mi
<
ainfo
->
nm
;
++
mi
)
{
auto
m
=
ainfo
->
mval
[
mi
];
ptrdiff_t
mvstart
=
ainfo
->
mvstart
[
mi
];
ptrdiff_t
stride
=
ainfo
->
stride
;
if
(
flags
&
SHARP_DP
)
for
(
size_t
l
=
m
;
l
<=
ainfo
->
lmax
;
++
l
)
reinterpret_cast
<
dcmplx
*>
(
alm
)[
mvstart
+
l
*
stride
]
=
0.
;
else
for
(
size_t
l
=
m
;
l
<=
ainfo
->
lmax
;
++
l
)
reinterpret_cast
<
fcmplx
*>
(
alm
)[
mvstart
+
l
*
stride
]
=
0.
;
}
for
(
size_t
mi
=
0
;
mi
<
mval_
.
size
();
++
mi
)
for
(
size_t
l
=
mval_
[
mi
];
l
<=
lmax_
;
++
l
)
reinterpret_cast
<
dcmplx
*>
(
alm
)[
mvstart
[
mi
]
+
l
*
stride
]
=
0.
;
}
void
sharp_standard_alm_info
::
clear_alm
(
fcmplx
*
alm
)
const
{
for
(
size_t
mi
=
0
;
mi
<
mval_
.
size
();
++
mi
)
for
(
size_t
l
=
mval_
[
mi
];
l
<=
lmax_
;
++
l
)
reinterpret_cast
<
fcmplx
*>
(
alm
)[
mvstart
[
mi
]
+
l
*
stride
]
=
0.
;
}
MRUTIL_NOINLINE
void
sharp_job
::
init_output
()
...
...
@@ -337,7 +333,8 @@ MRUTIL_NOINLINE void sharp_job::init_output()
if
(
flags
&
SHARP_ADD
)
return
;
if
(
type
==
SHARP_MAP2ALM
)
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
clear_alm
(
ainfo
,
alm
[
i
],
flags
);
(
flags
&
SHARP_DP
)
?
ainfo
->
clear_alm
(
reinterpret_cast
<
dcmplx
*>
(
alm
[
i
]))
:
ainfo
->
clear_alm
(
reinterpret_cast
<
fcmplx
*>
(
alm
[
i
]));
else
for
(
size_t
i
=
0
;
i
<
nmaps
;
++
i
)
(
flags
&
SHARP_DP
)
?
ginfo
->
clear_map
(
reinterpret_cast
<
double
*>
(
map
[
i
]))
...
...
@@ -368,99 +365,78 @@ void sharp_job::alloc_almtmp (size_t lmax, vector<dcmplx> &data)
almtmp
=
data
.
data
();
}
void
sharp_standard_alm_info
::
get_alm
(
size_t
mi
,
const
dcmplx
*
alm
,
dcmplx
*
almtmp
,
size_t
nalm
)
const
{
for
(
auto
l
=
mval_
[
mi
];
l
<=
lmax_
;
++
l
)
almtmp
[
nalm
*
l
]
=
alm
[
mvstart
[
mi
]
+
l
*
stride
];
}
void
sharp_standard_alm_info
::
get_alm
(
size_t
mi
,
const
fcmplx
*
alm
,
dcmplx
*
almtmp
,
size_t
nalm
)
const
{
for
(
auto
l
=
mval_
[
mi
];
l
<=
lmax_
;
++
l
)
almtmp
[
nalm
*
l
]
=
alm
[
mvstart
[
mi
]
+
l
*
stride
];
}
void
sharp_standard_alm_info
::
add_alm
(
size_t
mi
,
const
dcmplx
*
almtmp
,
dcmplx
*
alm
,
size_t
nalm
)
const
{
for
(
auto
l
=
mval_
[
mi
];
l
<=
lmax_
;
++
l
)
alm
[
mvstart
[
mi
]
+
l
*
stride
]
+=
almtmp
[
nalm
*
l
];
}
void
sharp_standard_alm_info
::
add_alm
(
size_t
mi
,
const
dcmplx
*
almtmp
,
fcmplx
*
alm
,
size_t
nalm
)
const
{
for
(
auto
l
=
mval_
[
mi
];
l
<=
lmax_
;
++
l
)
alm
[
mvstart
[
mi
]
+
l
*
stride
]
+=
fcmplx
(
almtmp
[
nalm
*
l
]);
}
MRUTIL_NOINLINE
void
sharp_job
::
alm2almtmp
(
size_t
lmax
,
size_t
mi
)
{
if
(
type
!=
SHARP_MAP2ALM
)
{
auto
ofs
=
ainfo
->
mvstart
[
mi
];
auto
stride
=
ainfo
->
stride
;
auto
m
=
ainfo
->
mval
[
mi
];
auto
m
=
ainfo
->
mval
(
mi
);
auto
lmin
=
(
m
<
spin
)
?
spin
:
m
;
if
(
spin
==
0
)
if
(
flags
&
SHARP_DP
)
{
if
(
flags
&
SHARP_DP
)
{
for
(
auto
l
=
m
;
l
<
lmin
;
++
l
)
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
almtmp
[
nalm
*
l
+
i
]
=
0
;
for
(
auto
l
=
lmin
;
l
<=
lmax
;
++
l
)
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
almtmp
[
nalm
*
l
+
i
]
=
reinterpret_cast
<
dcmplx
**>
(
alm
)[
i
][
ofs
+
l
*
stride
];
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
almtmp
[
nalm
*
(
lmax
+
1
)
+
i
]
=
0
;
}
else
{
for
(
auto
l
=
m
;
l
<
lmin
;
++
l
)
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
almtmp
[
nalm
*
l
+
i
]
=
0
;
for
(
auto
l
=
lmin
;
l
<=
lmax
;
++
l
)
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
almtmp
[
nalm
*
l
+
i
]
=
reinterpret_cast
<
fcmplx
**>
(
alm
)[
i
][
ofs
+
l
*
stride
];
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
ainfo
->
get_alm
(
mi
,
reinterpret_cast
<
dcmplx
**>
(
alm
)[
i
],
almtmp
+
i
,
nalm
);
for
(
auto
l
=
m
;
l
<
lmin
;
++
l
)
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
almtmp
[
nalm
*
(
lmax
+
1
)
+
i
]
=
0
;
}
almtmp
[
nalm
*
l
+
i
]
=
0
;
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
almtmp
[
nalm
*
(
lmax
+
1
)
+
i
]
=
0
;
}
else
{
if
(
flags
&
SHARP_DP
)
{
for
(
auto
l
=
m
;
l
<
lmin
;
++
l
)
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
almtmp
[
nalm
*
l
+
i
]
=
0
;
for
(
auto
l
=
lmin
;
l
<=
lmax
;
++
l
)
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
almtmp
[
nalm
*
l
+
i
]
=
reinterpret_cast
<
dcmplx
**>
(
alm
)[
i
][
ofs
+
l
*
stride
]
*
norm_l
[
l
];
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
almtmp
[
nalm
*
(
lmax
+
1
)
+
i
]
=
0
;
}
else
{
for
(
auto
l
=
m
;
l
<
lmin
;
++
l
)
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
almtmp
[
nalm
*
l
+
i
]
=
0
;
for
(
auto
l
=
lmin
;
l
<=
lmax
;
++
l
)
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
almtmp
[
nalm
*
l
+
i
]
=
dcmplx
(
reinterpret_cast
<
fcmplx
**>
(
alm
)[
i
][
ofs
+
l
*
stride
])
*
norm_l
[
l
];
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
ainfo
->
get_alm
(
mi
,
reinterpret_cast
<
fcmplx
**>
(
alm
)[
i
],
almtmp
+
i
,
nalm
);
for
(
auto
l
=
m
;
l
<
lmin
;
++
l
)
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
almtmp
[
nalm
*
(
lmax
+
1
)
+
i
]
=
0
;
}
almtmp
[
nalm
*
l
+
i
]
=
0
;
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
almtmp
[
nalm
*
(
lmax
+
1
)
+
i
]
=
0
;
}
if
(
spin
>
0
)
for
(
auto
l
=
lmin
;
l
<=
lmax
;
++
l
)
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
almtmp
[
nalm
*
l
+
i
]
*=
norm_l
[
l
];
}
else
for
(
size_t
i
=
nalm
*
ainfo
->
mval
[
mi
]
;
i
<
nalm
*
(
lmax
+
2
);
++
i
)
for
(
size_t
i
=
nalm
*
ainfo
->
mval
(
mi
)
;
i
<
nalm
*
(
lmax
+
2
);
++
i
)
almtmp
[
i
]
=
0
;
}
MRUTIL_NOINLINE
void
sharp_job
::
almtmp2alm
(
size_t
lmax
,
size_t
mi
)
{
if
(
type
!=
SHARP_MAP2ALM
)
return
;
auto
ofs
=
ainfo
->
mvstart
[
mi
];
auto
stride
=
ainfo
->
stride
;
auto
m
=
ainfo
->
mval
[
mi
];
auto
m
=
ainfo
->
mval
(
mi
);
auto
lmin
=
(
m
<
spin
)
?
spin
:
m
;
if
(
spin
==
0
)
{
if
(
flags
&
SHARP_DP
)
for
(
auto
l
=
lmin
;
l
<=
lmax
;
++
l
)
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
((
dcmplx
**
)
alm
)[
i
][
ofs
+
l
*
stride
]
+=
almtmp
[
nalm
*
l
+
i
];
else
for
(
auto
l
=
lmin
;
l
<=
lmax
;
++
l
)
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
((
fcmplx
**
)
alm
)[
i
][
ofs
+
l
*
stride
]
+=
fcmplx
(
almtmp
[
nalm
*
l
+
i
]);
}
if
(
spin
>
0
)
for
(
auto
l
=
lmin
;
l
<=
lmax
;
++
l
)
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
almtmp
[
nalm
*
l
+
i
]
*=
norm_l
[
l
];
if
(
flags
&
SHARP_DP
)
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
ainfo
->
add_alm
(
mi
,
almtmp
+
i
,
reinterpret_cast
<
dcmplx
**>
(
alm
)[
i
],
nalm
);
else
{
if
(
flags
&
SHARP_DP
)
for
(
auto
l
=
lmin
;
l
<=
lmax
;
++
l
)
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
((
dcmplx
**
)
alm
)[
i
][
ofs
+
l
*
stride
]
+=
almtmp
[
nalm
*
l
+
i
]
*
norm_l
[
l
];
else
for
(
auto
l
=
lmin
;
l
<=
lmax
;
++
l
)
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
((
fcmplx
**
)
alm
)[
i
][
ofs
+
l
*
stride
]
+=
fcmplx
(
almtmp
[
nalm
*
l
+
i
]
*
norm_l
[
l
]);
}
for
(
size_t
i
=
0
;
i
<
nalm
;
++
i
)
ainfo
->
add_alm
(
mi
,
almtmp
+
i
,
reinterpret_cast
<
fcmplx
**>
(
alm
)[
i
],
nalm
);
}
//virtual
...
...
@@ -496,26 +472,26 @@ void sharp_standard_geom_info::get_ring(bool weighted, size_t iring, const float
ringtmp
[
m
]
=
p1
[
m
*
stride
]
*
wgt
;
}
MRUTIL_NOINLINE
void
sharp_job
::
ringtmp2ring
(
const
sharp_geom_info
&
ginfo
,
size_t
iring
,
MRUTIL_NOINLINE
void
sharp_job
::
ringtmp2ring
(
size_t
iring
,
const
vector
<
double
>
&
ringtmp
,
ptrdiff_t
rstride
)
{
if
(
flags
&
SHARP_DP
)
for
(
size_t
i
=
0
;
i
<
nmaps
;
++
i
)
ginfo
.
add_ring
(
flags
&
SHARP_USE_WEIGHTS
,
iring
,
&
ringtmp
[
i
*
rstride
+
1
],
((
double
**
)
map
)[
i
]);
ginfo
->
add_ring
(
flags
&
SHARP_USE_WEIGHTS
,
iring
,
&
ringtmp
[
i
*
rstride
+
1
],
((
double
**
)
map
)[
i
]);
else
for
(
size_t
i
=
0
;
i
<
nmaps
;
++
i
)
ginfo
.
add_ring
(
flags
&
SHARP_USE_WEIGHTS
,
iring
,
&
ringtmp
[
i
*
rstride
+
1
],
((
float
**
)
map
)[
i
]);
ginfo
->
add_ring
(
flags
&
SHARP_USE_WEIGHTS
,
iring
,
&
ringtmp
[
i
*
rstride
+
1
],
((
float
**
)
map
)[
i
]);
}
MRUTIL_NOINLINE
void
sharp_job
::
ring2ringtmp
(
const
sharp_geom_info
&
ginfo
,
size_t
iring
,
MRUTIL_NOINLINE
void
sharp_job
::
ring2ringtmp
(
size_t
iring
,
vector
<
double
>
&
ringtmp
,
ptrdiff_t
rstride
)
{
if
(
flags
&
SHARP_DP
)
for
(
size_t
i
=
0
;
i
<
nmaps
;
++
i
)
ginfo
.
get_ring
(
flags
&
SHARP_USE_WEIGHTS
,
iring
,
((
double
**
)
map
)[
i
],
&
ringtmp
[
i
*
rstride
+
1
]);
ginfo
->
get_ring
(
flags
&
SHARP_USE_WEIGHTS
,
iring
,
((
double
**
)
map
)[
i
],
&
ringtmp
[
i
*
rstride
+
1
]);
else
for
(
size_t
i
=
0
;
i
<
nmaps
;
++
i
)
ginfo
.
get_ring
(
flags
&
SHARP_USE_WEIGHTS
,
iring
,
((
float
**
)
map
)[
i
],
&
ringtmp
[
i
*
rstride
+
1
]);
ginfo
->
get_ring
(
flags
&
SHARP_USE_WEIGHTS
,
iring
,
((
float
**
)
map
)[
i
],
&
ringtmp
[
i
*
rstride
+
1
]);
}
//FIXME: set phase to zero if not SHARP_MAP2ALM?
...
...
@@ -532,13 +508,13 @@ MRUTIL_NOINLINE void sharp_job::map2phase (size_t mmax, size_t llim, size_t ulim
while
(
auto
rng
=
sched
.
getNext
())
for
(
auto
ith
=
rng
.
lo
+
llim
;
ith
<
rng
.
hi
+
llim
;
++
ith
)
{
int
dim2
=
s_th
*
(
ith
-
llim
);
ring2ringtmp
(
*
ginfo
,
ginfo
->
pair
(
ith
).
r1
,
ringtmp
,
rstride
);
ring2ringtmp
(
ginfo
->
pair
(
ith
).
r1
,
ringtmp
,
rstride
);
for
(
size_t
i
=
0
;
i
<
nmaps
;
++
i
)
helper
.
ring2phase
(
*
ginfo
,
ginfo
->
pair
(
ith
).
r1
,
&
ringtmp
[
i
*
rstride
],
mmax
,
&
phase
[
dim2
+
2
*
i
],
pstride
);
if
(
ginfo
->
pair
(
ith
).
r2
!=~
size_t
(
0
))
{
ring2ringtmp
(
*
ginfo
,
ginfo
->
pair
(
ith
).
r2
,
ringtmp
,
rstride
);
ring2ringtmp
(
ginfo
->
pair
(
ith
).
r2
,
ringtmp
,
rstride
);
for
(
size_t
i
=
0
;
i
<
nmaps
;
++
i
)
helper
.
ring2phase
(
*
ginfo
,
ginfo
->
pair
(
ith
).
r2
,
&
ringtmp
[
i
*
rstride
],
mmax
,
&
phase
[
dim2
+
2
*
i
+
1
],
pstride
);
...
...
@@ -563,13 +539,13 @@ MRUTIL_NOINLINE void sharp_job::phase2map (size_t mmax, size_t llim, size_t ulim
for
(
size_t
i
=
0
;
i
<
nmaps
;
++
i
)
helper
.
phase2ring
(
*
ginfo
,
ginfo
->
pair
(
ith
).
r1
,
&
ringtmp
[
i
*
rstride
],
mmax
,
&
phase
[
dim2
+
2
*
i
],
pstride
);
ringtmp2ring
(
*
ginfo
,
ginfo
->
pair
(
ith
).
r1
,
ringtmp
,
rstride
);
ringtmp2ring
(
ginfo
->
pair
(
ith
).
r1
,
ringtmp
,
rstride
);
if
(
ginfo
->
pair
(
ith
).
r2
!=~
size_t
(
0
))
{
for
(
size_t
i
=
0
;
i
<
nmaps
;
++
i
)
helper
.
phase2ring
(
*
ginfo
,
ginfo
->
pair
(
ith
).
r2
,
&
ringtmp
[
i
*
rstride
],
mmax
,
&
phase
[
dim2
+
2
*
i
+
1
],
pstride
);
ringtmp2ring
(
*
ginfo
,
ginfo
->
pair
(
ith
).
r2
,
ringtmp
,
rstride
);
ringtmp2ring
(
ginfo
->
pair
(
ith
).
r2
,
ringtmp
,
rstride
);
}
}
});
/* end of parallel region */
...
...
@@ -579,8 +555,8 @@ MRUTIL_NOINLINE void sharp_job::execute()
{
mr
::
SimpleTimer
timer
;
opcnt
=
0
;
size_t
lmax
=
ainfo
->
lmax
,
mmax
=
sharp_get_mmax
(
ainfo
->
m
val
);
size_t
lmax
=
ainfo
->
lmax
()
,
mmax
=
ainfo
->
m
max
(
);
norm_l
=
(
type
==
SHARP_ALM2MAP_DERIV1
)
?
sharp_Ylmgen
::
get_d1norm
(
lmax
)
:
...
...
@@ -615,7 +591,7 @@ MRUTIL_NOINLINE void sharp_job::execute()
/* map->phase where necessary */
map2phase
(
mmax
,
llim
,
ulim
);
mr
::
execDynamic
(
ainfo
->
nm
,
0
,
1
,
[
&
](
mr
::
Scheduler
&
sched
)
mr
::
execDynamic
(
ainfo
->
nm
()
,
0
,
1
,
[
&
](
mr
::
Scheduler
&
sched
)
{
sharp_job
ljob
=
*
this
;
ljob
.
opcnt
=
0
;
...
...
@@ -654,7 +630,7 @@ void sharp_job::build_common (sharp_jobtype type,
if
(
type
==
SHARP_Yt
)
type
=
SHARP_MAP2ALM
;
if
(
type
==
SHARP_WY
)
{
type
=
SHARP_ALM2MAP
;
flags
|=
SHARP_USE_WEIGHTS
;
}
MR_assert
(
spin
<=
alm_info
.
lmax
,
"bad spin"
);
MR_assert
(
spin
<=
alm_info
.
lmax
()
,
"bad spin"
);
this
->
type
=
type
;
this
->
spin
=
spin
;
nmaps
=
(
type
==
SHARP_ALM2MAP_DERIV1
)
?
2
:
((
spin
>
0
)
?
2
:
1
);
...
...
libsharp2/sharp.h
View file @
aaea79a7
...
...
@@ -28,6 +28,7 @@
#ifndef SHARP_SHARP_H
#define SHARP_SHARP_H
#include
<complex>
#include
<cstddef>
#include
<vector>
#include
<memory>
...
...
@@ -107,48 +108,74 @@ class sharp_standard_geom_info: public sharp_geom_info
/*! \defgroup almgroup Helpers for dealing with a_lm */
/*! \{ */
class
sharp_alm_info
{
public:
~
sharp_alm_info
()
{}
virtual
size_t
lmax
()
const
=
0
;
virtual
size_t
mmax
()
const
=
0
;
virtual
size_t
nm
()
const
=
0
;
virtual
size_t
mval
(
size_t
i
)
const
=
0
;
virtual
void
clear_alm
(
std
::
complex
<
double
>
*
alm
)
const
=
0
;
virtual
void
clear_alm
(
std
::
complex
<
float
>
*
alm
)
const
=
0
;
virtual
void
get_alm
(
size_t
mi
,
const
std
::
complex
<
double
>
*
alm
,
std
::
complex
<
double
>
*
almtmp
,
size_t
nalm
)
const
=
0
;
virtual
void
get_alm
(
size_t
mi
,
const
std
::
complex
<
float
>
*
alm
,
std
::
complex
<
double
>
*
almtmp
,
size_t
nalm
)
const
=
0
;
virtual
void
add_alm
(
size_t
mi
,
const
std
::
complex
<
double
>
*
almtmp
,
std
::
complex
<
double
>
*
alm
,
size_t
nalm
)
const
=
0
;
virtual
void
add_alm
(
size_t
mi
,
const
std
::
complex
<
double
>
*
almtmp
,
std
::
complex
<
float
>
*
alm
,
size_t
nalm
)
const
=
0
;
};
/*! \internal
Helper type for index calculation in a_lm arrays. */
struct
sharp_alm_info
class
sharp_standard_alm_info
:
public
sharp_alm_info
{
/*! Maximum \a l index of the array */
size_t
lmax
;
/*! Number of different \a m values in this object */
size_t
nm
;
/*! Array with \a nm entries containing the individual m values */
std
::
vector
<
size_t
>
mval
;
/*! Array with \a nm entries containing the (hypothetical) indices of
the coefficients with quantum numbers 0,\a mval[i] */
std
::
vector
<
ptrdiff_t
>
mvstart
;
/*! Stride between a_lm and a_(l+1),m */
ptrdiff_t
stride
;
/*! Creates an a_lm data structure from the following parameters:
\param lmax maximum \a l quantum number (>=0)
\param mmax maximum \a m quantum number (0<= \a mmax <= \a lmax)
\param stride the stride between entries with identical \a m, and \a l
differing by 1.
\param mstart the index of the (hypothetical) coefficient with the
quantum numbers 0,\a m. Must have \a mmax+1 entries.
*/
sharp_alm_info
(
size_t
lmax_
,
size_t
mmax
,
ptrdiff_t
stride_
,
const
ptrdiff_t
*
mstart
);
/*! Creates an a_lm data structure which from the following parameters:
\param lmax maximum \a l quantum number (\a >=0)
\param nm number of different \a m (\a 0<=nm<=lmax+1)
\param stride the stride between entries with identical \a m, and \a l
differing by 1.
\param mval array with \a nm entries containing the individual m values
\param mvstart array with \a nm entries containing the (hypothetical)
indices of the coefficients with the quantum numbers 0,\a mval[i]
*/
sharp_alm_info
(
size_t
lmax_
,
size_t
nm_
,
ptrdiff_t
stride_
,
const
size_t
*
mval_
,
const
ptrdiff_t
*
mvstart
);
/*! Returns the index of the coefficient with quantum numbers \a l,
\a mval_[mi].
\note for a \a sharp_alm_info generated by sharp_make_alm_info() this is
the index for the coefficient with the quantum numbers \a l, \a mi. */
ptrdiff_t
index
(
int
l
,
int
mi
);
private:
/*! Maximum \a l index of the array */
size_t
lmax_
;
/*! Array with \a nm entries containing the individual m values */
std
::
vector
<
size_t
>
mval_
;
/*! Array with \a nm entries containing the (hypothetical) indices of
the coefficients with quantum numbers 0,\a mval[i] */
std
::
vector
<
ptrdiff_t
>
mvstart
;
/*! Stride between a_lm and a_(l+1),m */
ptrdiff_t
stride
;
public:
/*! Creates an a_lm data structure from the following parameters:
\param lmax maximum \a l quantum number (>=0)
\param mmax maximum \a m quantum number (0<= \a mmax <= \a lmax)
\param stride the stride between entries with identical \a m, and \a l
differing by 1.
\param mstart the index of the (hypothetical) coefficient with the
quantum numbers 0,\a m. Must have \a mmax+1 entries.
*/
sharp_standard_alm_info
(
size_t
lmax__
,
size_t
mmax_
,
ptrdiff_t
stride_
,
const
ptrdiff_t
*
mstart
);
/*! Creates an a_lm data structure which from the following parameters:
\param lmax maximum \a l quantum number (\a >=0)
\param nm number of different \a m (\a 0<=nm<=lmax+1)
\param stride the stride between entries with identical \a m, and \a l
differing by 1.
\param mval array with \a nm entries containing the individual m values
\param mvstart array with \a nm entries containing the (hypothetical)
indices of the coefficients with the quantum numbers 0,\a mval[i]
*/
sharp_standard_alm_info
(
size_t
lmax__
,
size_t
nm__
,
ptrdiff_t
stride_
,
const
size_t
*
mval__
,
const
ptrdiff_t
*
mvstart_
);
/*! Returns the index of the coefficient with quantum numbers \a l,
\a mval_[mi].
\note for a \a sharp_alm_info generated by sharp_make_alm_info() this is
the index for the coefficient with the quantum numbers \a l, \a mi. */
ptrdiff_t
index
(
int
l
,
int
mi
);
virtual
size_t
lmax
()
const
{
return
lmax_
;
}
virtual
size_t
mmax
()
const
;
virtual
size_t
nm
()
const
{
return
mval_
.
size
();
}
virtual
size_t
mval
(
size_t
i
)
const
{
return
mval_
[
i
];
}
virtual
void
clear_alm
(
std
::
complex
<
double
>
*
alm
)
const
;
virtual
void
clear_alm
(
std
::
complex
<
float
>
*
alm
)
const
;
virtual
void
get_alm
(
size_t
mi
,
const
std
::
complex
<
double
>
*
alm
,
std
::
complex
<
double
>
*
almtmp
,
size_t
nalm
)
const
;
virtual
void
get_alm
(
size_t
mi
,
const
std
::
complex
<
float
>
*
alm
,
std
::
complex
<
double
>
*
almtmp
,
size_t
nalm
)
const
;
virtual
void
add_alm
(
size_t
mi
,
const
std
::
complex
<
double
>
*
almtmp
,
std
::
complex
<
double
>
*
alm
,
size_t
nalm
)
const
;
virtual
void
add_alm
(
size_t
mi
,
const
std
::
complex
<
double
>
*
almtmp
,
std
::
complex
<
float
>
*
alm
,
size_t
nalm
)
const
;
};
/*! \} */
...
...
libsharp2/sharp_almhelpers.cc
View file @
aaea79a7
...
...
@@ -29,19 +29,19 @@
using
namespace
std
;
unique_ptr
<
sharp_alm_info
>
sharp_make_triangular_alm_info
(
int
lmax
,
int
mmax
,
int
stride
)
unique_ptr
<
sharp_
standard_
alm_info
>
sharp_make_triangular_alm_info
(
int
lmax
,
int
mmax
,
int
stride
)
{
vector
<
ptrdiff_t
>
mvstart
(
mmax
+
1
);
ptrdiff_t
tval
=
2
*
lmax
+
1
;
for
(
ptrdiff_t
m
=
0
;
m
<=
mmax
;
++
m
)
mvstart
[
m
]
=
stride
*
((
m
*
(
tval
-
m
))
>>
1
);
return
make_unique
<
sharp_alm_info
>
(
lmax
,
mmax
,
stride
,
mvstart
.
data
());
return
make_unique
<
sharp_
standard_
alm_info
>
(
lmax
,
mmax
,
stride
,
mvstart
.
data
());
}
unique_ptr
<
sharp_alm_info
>
sharp_make_rectangular_alm_info
(
int
lmax
,
int
mmax
,
int
stride
)
unique_ptr
<
sharp_
standard_
alm_info
>
sharp_make_rectangular_alm_info
(
int
lmax
,
int
mmax
,
int
stride
)
{
vector
<
ptrdiff_t
>
mvstart
(
mmax
+
1
);
for
(
ptrdiff_t
m
=
0
;
m
<=
mmax
;
++
m
)
mvstart
[
m
]
=
stride
*
m
*
(
lmax
+
1
);
return
make_unique
<
sharp_alm_info
>
(
lmax
,
mmax
,
stride
,
mvstart
.
data
());
return
make_unique
<
sharp_
standard_
alm_info
>
(
lmax
,
mmax
,
stride
,
mvstart
.
data
());
}
libsharp2/sharp_almhelpers.h
View file @
aaea79a7
...
...
@@ -35,11 +35,11 @@
/*! Initialises an a_lm data structure according to the scheme used by
Healpix_cxx.
\ingroup almgroup */
std
::
unique_ptr
<
sharp_alm_info
>
sharp_make_triangular_alm_info
(
int
lmax
,
int
mmax
,
int
stride
);
std
::
unique_ptr
<
sharp_
standard_
alm_info
>
sharp_make_triangular_alm_info
(
int
lmax
,
int
mmax
,
int
stride
);
/*! Initialises an a_lm data structure according to the scheme used by
Fortran Healpix
\ingroup almgroup */
std
::
unique_ptr
<
sharp_alm_info
>
sharp_make_rectangular_alm_info
(
int
lmax
,
int
mmax
,
int
stride
);
std
::
unique_ptr
<
sharp_
standard_
alm_info
>
sharp_make_rectangular_alm_info
(
int
lmax
,
int
mmax
,
int
stride
);
#endif
libsharp2/sharp_core_inc.cc
View file @
aaea79a7
...
...
@@ -943,7 +943,7 @@ MRUTIL_NOINLINE static void inner_loop_a2m(sharp_job &job, const int *ispair,
const
double
*
cth_
,
const
double
*
sth_
,
size_t
llim
,
size_t
ulim
,
sharp_Ylmgen
&
gen
,
size_t
mi
,
const
size_t
*
mlim
)
{
const
size_t
m
=
job
.
ainfo
->
mval
[
mi
]
;
const
size_t
m
=
job
.
ainfo
->
mval
(
mi
)
;
gen
.
prepare
(
m
);
switch
(
job
.
type
)
...
...
@@ -1102,7 +1102,7 @@ MRUTIL_NOINLINE static void inner_loop_m2a(sharp_job &job, const int *ispair,
const
double
*
cth_
,
const
double
*
sth_
,
size_t
llim
,
size_t
ulim
,
sharp_Ylmgen
&
gen
,
size_t
mi
,
const
size_t
*
mlim
)
{
const
size_t
m
=
job
.
ainfo
->
mval
[
mi
]
;
const
size_t
m
=
job
.
ainfo
->
mval
(
mi
)
;
gen
.
prepare
(
m
);
switch
(
job
.
type
)
...
...
libsharp2/sharp_internal.h
View file @
aaea79a7
...
...
@@ -58,9 +58,9 @@ struct sharp_job
void
init_output
();
void
alm2almtmp
(
size_t
lmax
,
size_t
mi
);
void
almtmp2alm
(
size_t
lmax
,
size_t
mi
);
void
ring2ringtmp
(
const
sharp_geom_info
&
ginfo
,
size_t
iring
,
std
::
vector
<
double
>
&
ringtmp
,
void
ring2ringtmp
(
size_t
iring
,
std
::
vector
<
double
>
&
ringtmp
,
ptrdiff_t
rstride
);
void
ringtmp2ring
(
const
sharp_geom_info
&
ginfo
,
size_t
iring
,
const
std
::
vector
<
double
>
&
ringtmp
,
ptrdiff_t
rstride
);
void
ringtmp2ring
(
size_t
iring
,
const
std
::
vector
<
double
>
&
ringtmp
,
ptrdiff_t
rstride
);
void
map2phase
(
size_t
mmax
,
size_t
llim
,
size_t
ulim
);
void
phase2map
(
size_t
mmax
,
size_t
llim
,
size_t
ulim
);
void
execute
();
...
...
test/sharp2_testsuite.cc
View file @
aaea79a7
...
...
@@ -86,14 +86,14 @@ static double drand (double min, double max, unsigned *state)
return
min
+
(
max
-
min
)
*
(
*
state
)
/
(
0x7fffffff
+
1.0
);
}
static
void
random_alm
(
dcmplx
*
alm
,
sharp_alm_info
&
helper
,
size_t
spin
,
size_t
cnt
)
static
void
random_alm
(
dcmplx
*
alm
,
sharp_
standard_
alm_info
&
helper
,
size_t
spin
,
size_t
cnt
)
{
// FIXME: re-introduce multi-threading?
for
(
size_t
mi
=
0
;
mi
<
helper
.
nm
;
++
mi
)
for
(
size_t
mi
=
0
;
mi
<
helper
.
nm
()
;
++
mi
)
{
auto
m
=
helper
.
mval
[
mi
]
;
auto
m
=
helper
.
mval
(
mi
)
;
unsigned
state
=
1234567u
*
unsigned
(
cnt
)
+
8912u
*
unsigned
(
m
);
// random seed
for
(
auto
l
=
m
;
l
<=
helper
.
lmax
;
++
l
)
for
(
auto
l
=
m
;
l
<=
helper
.
lmax
()
;
++
l
)
{
if
((
l
<
spin
)
&&
(
m
<
spin
))
alm
[
helper
.
index
(
l
,
mi
)]
=
0.
;
...
...
@@ -130,8 +130,8 @@ static double totalMem()
static
ptrdiff_t
get_nalms
(
const
sharp_alm_info
&
ainfo
)
{
ptrdiff_t
res
=
0
;
for
(
size_t
i
=
0
;
i
<
ainfo
.
nm
;
++
i
)
res
+=
ainfo
.
lmax
-
ainfo
.
mval
[
i
]
+
1
;
for
(
size_t
i
=
0
;
i
<
ainfo
.
nm
()
;
++
i
)
res
+=
ainfo
.
lmax
()
-
ainfo
.
mval
(
i
)
+
1
;
return
res
;
}
...
...
@@ -200,7 +200,7 @@ static int good_fft_size(int n)
}
static
void
get_infos
(
const
string
&
gname
,
int
lmax
,
int
&
mmax
,
int
&
gpar1
,
int
&
gpar2
,
unique_ptr
<
sharp_geom_info
>
&
ginfo
,
unique_ptr
<
sharp_alm_info
>
&
ainfo
,
int
verbose
)
int
&
gpar2
,
unique_ptr
<
sharp_geom_info
>
&
ginfo
,
unique_ptr
<
sharp_
standard_
alm_info
>
&
ainfo
,
int
verbose
)
{
MR_assert
(
lmax
>=
0
,
"lmax must not be negative"
);
if
(
mmax
<
0
)
mmax
=
lmax
;
...
...
@@ -365,7 +365,7 @@ static void check_sign_scale(void)
"error"
);
}
static
void
do_sht
(
sharp_geom_info
&
ginfo
,
sharp_alm_info
&
ainfo
,
static
void
do_sht
(
sharp_geom_info
&
ginfo
,
sharp_
standard_
alm_info
&
ainfo
,