From 549d1c35e1777f7a9ba2d9380a18e4e1b0815b3e Mon Sep 17 00:00:00 2001
From: Martin Reinecke <martin@mpa-garching.mpg.de>
Date: Mon, 5 Nov 2012 13:13:20 +0100
Subject: [PATCH] tweak alm <-> alm_tmp

---
 libsharp/sharp.c | 75 ++++++++++++++++++++++++++++++++++++------------
 1 file changed, 56 insertions(+), 19 deletions(-)

diff --git a/libsharp/sharp.c b/libsharp/sharp.c
index 69a4bfa..d833516 100644
--- a/libsharp/sharp.c
+++ b/libsharp/sharp.c
@@ -415,18 +415,36 @@ static void dealloc_almtmp (sharp_job *job)
 static void alm2almtmp (sharp_job *job, int lmax, int mi)
   {
   if (job->type!=SHARP_MAP2ALM)
-    for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)
+    {
+    ptrdiff_t ofs=job->ainfo->mvstart[mi];
+    int stride=job->ainfo->stride;
+    if (job->spin==0)
       {
-      ptrdiff_t aidx = sharp_alm_index(job->ainfo,l,mi);
-      double fct = job->norm_l[l];
-      for (int i=0; i<job->ntrans*job->nalm; ++i)
-        if (job->fde==DOUBLE)
-          job->almtmp[job->ntrans*job->nalm*l+i]
-            = ((dcmplx *)job->alm[i])[aidx]*fct;
-        else
-          job->almtmp[job->ntrans*job->nalm*l+i]
-            = ((fcmplx *)job->alm[i])[aidx]*fct;
+      if (job->fde==DOUBLE)
+        for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)
+          for (int i=0; i<job->ntrans*job->nalm; ++i)
+            job->almtmp[job->ntrans*job->nalm*l+i]
+              = ((dcmplx *)job->alm[i])[ofs+l*stride];
+      else
+        for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)
+          for (int i=0; i<job->ntrans*job->nalm; ++i)
+            job->almtmp[job->ntrans*job->nalm*l+i]
+              = ((fcmplx *)job->alm[i])[ofs+l*stride];
       }
+    else
+      {
+      if (job->fde==DOUBLE)
+        for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)
+          for (int i=0; i<job->ntrans*job->nalm; ++i)
+            job->almtmp[job->ntrans*job->nalm*l+i]
+              = ((dcmplx *)job->alm[i])[ofs+l*stride]*job->norm_l[l];
+      else
+        for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)
+          for (int i=0; i<job->ntrans*job->nalm; ++i)
+            job->almtmp[job->ntrans*job->nalm*l+i]
+              = ((fcmplx *)job->alm[i])[ofs+l*stride]*job->norm_l[l];
+      }
+    }
   else
     SET_ARRAY(job->almtmp,job->ntrans*job->nalm*job->ainfo->mval[mi],
               job->ntrans*job->nalm*(lmax+1),0.);
@@ -435,16 +453,33 @@ static void alm2almtmp (sharp_job *job, int lmax, int mi)
 static void almtmp2alm (sharp_job *job, int lmax, int mi)
   {
   if (job->type != SHARP_MAP2ALM) return;
-  for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)
+  ptrdiff_t ofs=job->ainfo->mvstart[mi];
+  int stride=job->ainfo->stride;
+  if (job->spin==0)
     {
-    ptrdiff_t aidx = sharp_alm_index(job->ainfo,l,mi);
-    for (int i=0;i<job->ntrans*job->nalm;++i)
-      if (job->fde==DOUBLE)
-        ((dcmplx *)job->alm[i])[aidx] +=
-          job->almtmp[job->ntrans*job->nalm*l+i]*job->norm_l[l];
-      else
-        ((fcmplx *)job->alm[i])[aidx] +=
-          (fcmplx)(job->almtmp[job->ntrans*job->nalm*l+i]*job->norm_l[l]);
+    if (job->fde==DOUBLE)
+      for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)
+        for (int i=0;i<job->ntrans*job->nalm;++i)
+          ((dcmplx *)job->alm[i])[ofs+l*stride] +=
+            job->almtmp[job->ntrans*job->nalm*l+i];
+    else
+      for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)
+        for (int i=0;i<job->ntrans*job->nalm;++i)
+          ((fcmplx *)job->alm[i])[ofs+l*stride] +=
+            (fcmplx)(job->almtmp[job->ntrans*job->nalm*l+i]);
+    }
+  else
+    {
+    if (job->fde==DOUBLE)
+      for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)
+        for (int i=0;i<job->ntrans*job->nalm;++i)
+          ((dcmplx *)job->alm[i])[ofs+l*stride] +=
+            job->almtmp[job->ntrans*job->nalm*l+i]*job->norm_l[l];
+    else
+      for (int l=job->ainfo->mval[mi]; l<=lmax; ++l)
+        for (int i=0;i<job->ntrans*job->nalm;++i)
+          ((fcmplx *)job->alm[i])[ofs+l*stride] +=
+            (fcmplx)(job->almtmp[job->ntrans*job->nalm*l+i]*job->norm_l[l]);
     }
   }
 
@@ -600,6 +635,8 @@ static int sharp_oracle (sharp_jobtype type, int spin, int ntrans)
   int nrings=(lmax+1)/4;
   int ppring=1;
 
+  spin = (spin!=0) ? 2 : 0;
+
   ptrdiff_t npix=(ptrdiff_t)nrings*ppring;
   sharp_geom_info *tinfo;
   sharp_make_gauss_geom_info (nrings, ppring, 1, ppring, &tinfo);