From 526426388ef470b69d02d858b89ca10ec86d11ae Mon Sep 17 00:00:00 2001 From: Martin Reinecke Date: Thu, 17 Jan 2013 11:45:43 +0100 Subject: [PATCH] output a line for easy use with gnuplot --- libsharp/sharp_testsuite.c | 101 +++++++++++++++++++++++-------------- 1 file changed, 62 insertions(+), 39 deletions(-) diff --git a/libsharp/sharp_testsuite.c b/libsharp/sharp_testsuite.c index 866a6b6..84864da 100644 --- a/libsharp/sharp_testsuite.c +++ b/libsharp/sharp_testsuite.c @@ -34,6 +34,9 @@ #include "mpi.h" #include "sharp_mpi.h" #endif +#ifdef _OPENMP +#include +#endif #include "sharp.h" #include "sharp_geomhelpers.h" #include "sharp_almhelpers.h" @@ -41,6 +44,7 @@ #include "sharp_announce.h" #include "sharp_core.h" #include "memusage.h" +#include "sharp_vecsupport.h" typedef complex double dcmplx; @@ -252,65 +256,63 @@ static int sharp_get_mlim (int lmax, int spin, double sth, double cth, return (int)(res+0.5); } -static void get_infos (const char *gname, int lmax, int mmax, int gpar1, - int gpar2, sharp_geom_info **ginfo, sharp_alm_info **ainfo) +static void get_infos (const char *gname, int lmax, int *mmax, int *gpar1, + int *gpar2, sharp_geom_info **ginfo, sharp_alm_info **ainfo) { UTIL_ASSERT(lmax>=0,"lmax must not be negative"); - if (mmax<0) mmax=lmax; - UTIL_ASSERT(mmax<=lmax,"mmax larger than lmax"); + if (*mmax<0) *mmax=lmax; + UTIL_ASSERT(*mmax<=lmax,"mmax larger than lmax"); - if (mytask==0) printf ("lmax: %d, mmax: %d\n",lmax,mmax); + if (mytask==0) printf ("lmax: %d, mmax: %d\n",lmax,*mmax); - sharp_make_triangular_alm_info(lmax,mmax,1,ainfo); + sharp_make_triangular_alm_info(lmax,*mmax,1,ainfo); #ifdef USE_MPI reduce_alm_info(*ainfo); #endif if (strcmp(gname,"healpix")==0) { - int nside=gpar1; - if (nside<1) nside=lmax/2; - if (nside==0) ++nside; - sharp_make_healpix_geom_info (nside, 1, ginfo); - if (mytask==0) printf ("HEALPix grid, nside=%d\n",nside); + if (*gpar1<1) *gpar1=lmax/2; + if (*gpar1==0) ++(*gpar1); + sharp_make_healpix_geom_info (*gpar1, 1, ginfo); + if (mytask==0) printf ("HEALPix grid, nside=%d\n",*gpar1); } else if (strcmp(gname,"gauss")==0) { - int nlat=gpar1, nlon=gpar2; - if (nlat<1) nlat=lmax+1; - if (nlon<1) nlon=2*mmax+1; - sharp_make_gauss_geom_info (nlat, nlon, 0., 1, nlon, ginfo); - if (mytask==0) printf ("Gauss-Legendre grid, nlat=%d, nlon=%d\n",nlat,nlon); + if (*gpar1<1) *gpar1=lmax+1; + if (*gpar2<1) *gpar2=2*(*mmax)+1; + sharp_make_gauss_geom_info (*gpar1, *gpar2, 0., 1, *gpar2, ginfo); + if (mytask==0) + printf ("Gauss-Legendre grid, nlat=%d, nlon=%d\n",*gpar1,*gpar2); } else if (strcmp(gname,"fejer1")==0) { - int nlat=gpar1, nlon=gpar2; - if (nlat<1) nlat=2*lmax+1; - if (nlon<1) nlon=2*mmax+1; - sharp_make_fejer1_geom_info (nlat, nlon, 0., 1, nlon, ginfo); - if (mytask==0) printf ("Fejer1 grid, nlat=%d, nlon=%d\n",nlat,nlon); + if (*gpar1<1) *gpar1=2*lmax+1; + if (*gpar2<1) *gpar2=2*(*mmax)+1; + sharp_make_fejer1_geom_info (*gpar1, *gpar2, 0., 1, *gpar2, ginfo); + if (mytask==0) printf ("Fejer1 grid, nlat=%d, nlon=%d\n",*gpar1,*gpar2); } else if (strcmp(gname,"fejer2")==0) { - int nlat=gpar1, nlon=gpar2; - if (nlat<1) nlat=2*lmax+1; - if (nlon<1) nlon=2*mmax+1; - sharp_make_fejer2_geom_info (nlat, nlon, 0., 1, nlon, ginfo); - if (mytask==0) printf ("Fejer2 grid, nlat=%d, nlon=%d\n",nlat,nlon); + if (*gpar1<1) *gpar1=2*lmax+1; + if (*gpar2<1) *gpar2=2*(*mmax)+1; + sharp_make_fejer2_geom_info (*gpar1, *gpar2, 0., 1, *gpar2, ginfo); + if (mytask==0) printf ("Fejer2 grid, nlat=%d, nlon=%d\n",*gpar1,*gpar2); } else if (strcmp(gname,"cc")==0) { - int nlat=gpar1, nlon=gpar2; - if (nlat<1) nlat=2*lmax+1; - if (nlon<1) nlon=2*mmax+1; - sharp_make_fejer1_geom_info (nlat, nlon, 0., 1, nlon, ginfo); - if (mytask==0) printf("Clenshaw-Curtis grid, nlat=%d, nlon=%d\n",nlat,nlon); + if (*gpar1<1) *gpar1=2*lmax+1; + if (*gpar2<1) *gpar2=2*(*mmax)+1; + sharp_make_cc_geom_info (*gpar1, *gpar2, 0., 1, *gpar2, ginfo); + if (mytask==0) + printf("Clenshaw-Curtis grid, nlat=%d, nlon=%d\n",*gpar1,*gpar2); } else if (strcmp(gname,"smallgauss")==0) { - int nlat=gpar1, nlon=gpar2; + int nlat=*gpar1, nlon=*gpar2; if (nlat<1) nlat=lmax+1; - if (nlon<1) nlon=2*mmax+1; + if (nlon<1) nlon=2*(*mmax)+1; + *gpar1=nlat; *gpar2=nlon; sharp_make_gauss_geom_info (nlat, nlon, 0., 1, nlon, ginfo); ptrdiff_t npix_o=get_npix(*ginfo); size_t ofs=0; @@ -526,7 +528,8 @@ static void sharp_acctest(void) sharp_geom_info *ginfo; sharp_alm_info *ainfo; - get_infos ("gauss", 127, 127, 128, 256, &ginfo, &ainfo); + int lmax=127, mmax=127, nlat=128, nlon=256; + get_infos ("gauss", lmax, &mmax, &nlat, &nlon, &ginfo, &ainfo); for (int nv=1; nv<=6; ++nv) for (int ntrans=1; ntrans<=6; ++ntrans) { @@ -557,7 +560,7 @@ static void sharp_test (int argc, const char **argv) sharp_geom_info *ginfo; sharp_alm_info *ainfo; - get_infos (argv[2], lmax, mmax, gpar1, gpar2, &ginfo, &ainfo); + get_infos (argv[2], lmax, &mmax, &gpar1, &gpar2, &ginfo, &ainfo); int ncomp = ntrans*((spin==0) ? 1 : 2); double t_a2m, t_m2a; @@ -575,9 +578,6 @@ static void sharp_test (int argc, const char **argv) for (int i=0; i