From 9cf2f7dffd17a21888f94cd22e1428dc39ca0a3c Mon Sep 17 00:00:00 2001
From: ims <ims@strw.leidenuniv.nl>
Date: Thu, 20 Apr 2017 10:27:24 +0200
Subject: [PATCH 1/2] Slightly rearranged stateq.

In stateq.c I moved the function getmatrix() to before stateq(), which allows its template to be removed from lime.h.
---
 src/lime.h   |   2 +-
 src/stateq.c | 170 ++++++++++++++++++++++++++-------------------------
 2 files changed, 87 insertions(+), 85 deletions(-)

diff --git a/src/lime.h b/src/lime.h
index b8973a1..d6e3005 100644
--- a/src/lime.h
+++ b/src/lime.h
@@ -337,7 +337,7 @@ void	getclosest(double, double, double, long*, long*, double*, double*, double*)
 double	geterf(const double, const double);
 void	getjbar(int, molData*, struct grid*, const int, configInfo*, struct blendInfo, int, gridPointData*, double*);
 void	getMass(configInfo*, struct grid*, const gsl_rng*);
-void	getmatrix(int, gsl_matrix*, molData*, struct grid*, int, gridPointData*);
+//void	getmatrix(int, gsl_matrix*, molData*, struct grid*, int, gridPointData*);
 void	getVelocities(configInfo *, struct grid *);
 void	getVelocities_pregrid(configInfo *, struct grid *);
 void	gridPopsInit(configInfo*, molData*, struct grid*);
diff --git a/src/stateq.c b/src/stateq.c
index bd64125..a669086 100644
--- a/src/stateq.c
+++ b/src/stateq.c
@@ -13,6 +13,92 @@
 #include <gsl/gsl_errno.h>
 
 
+/*....................................................................*/
+void
+getmatrix(int id, gsl_matrix *matrix, molData *md, struct grid *gp, int ispec, gridPointData *mp){
+  int ti,k,l,li,ipart,di;
+  struct getmatrix {
+    double *ctot;
+    gsl_matrix * colli;
+  } *partner;
+
+  partner = malloc(sizeof(struct getmatrix)*md[ispec].npart);
+
+  /* Initialize matrix with zeros */
+  for(ipart=0;ipart<md[ispec].npart;ipart++){
+    partner[ipart].colli = gsl_matrix_alloc(md[ispec].nlev+1,md[ispec].nlev+1);
+    if(md[ispec].nlev>0) partner[ipart].ctot  = malloc(sizeof(double)*md[ispec].nlev);
+    else {
+      if(!silent)bail_out("Matrix initialization error in stateq");
+      exit(0);
+    }
+    for(k=0;k<md[ispec].nlev+1;k++){
+      for(l=0;l<md[ispec].nlev+1;l++){
+        gsl_matrix_set(matrix, k, l, 0.);
+        gsl_matrix_set(partner[ipart].colli, k, l, 0.);
+      }
+    }
+  }
+
+  /* Populate matrix with radiative transitions */
+  for(li=0;li<md[ispec].nline;li++){
+    k=md[ispec].lau[li];
+    l=md[ispec].lal[li];
+    gsl_matrix_set(matrix, k, k, gsl_matrix_get(matrix, k, k)+md[ispec].beinstu[li]*mp[ispec].jbar[li]+md[ispec].aeinst[li]);
+    gsl_matrix_set(matrix, l, l, gsl_matrix_get(matrix, l, l)+md[ispec].beinstl[li]*mp[ispec].jbar[li]);
+    gsl_matrix_set(matrix, k, l, gsl_matrix_get(matrix, k, l)-md[ispec].beinstl[li]*mp[ispec].jbar[li]);
+    gsl_matrix_set(matrix, l, k, gsl_matrix_get(matrix, l, k)-md[ispec].beinstu[li]*mp[ispec].jbar[li]-md[ispec].aeinst[li]);
+  }
+
+  /* Populate matrix with collisional transitions */
+  for(ipart=0;ipart<md[ispec].npart;ipart++){
+    struct cpData part = md[ispec].part[ipart];
+    double *downrates = part.down;
+    for(ti=0;ti<part.ntrans;ti++){
+      int coeff_index = ti*part.ntemp + gp[id].mol[ispec].partner[ipart].t_binlow;
+      double down = downrates[coeff_index]\
+                  + gp[id].mol[ispec].partner[ipart].interp_coeff*(downrates[coeff_index+1]\
+                  - downrates[coeff_index]);
+      double up = down*md[ispec].gstat[part.lcu[ti]]/md[ispec].gstat[part.lcl[ti]]\
+                *exp(-HCKB*(md[ispec].eterm[part.lcu[ti]]-md[ispec].eterm[part.lcl[ti]])/gp[id].t[0]);
+      gsl_matrix_set(partner[ipart].colli, part.lcu[ti], part.lcl[ti], down);
+      gsl_matrix_set(partner[ipart].colli, part.lcl[ti], part.lcu[ti], up);
+    }
+
+    for(k=0;k<md[ispec].nlev;k++){
+      partner[ipart].ctot[k]=0.;
+      for(l=0;l<md[ispec].nlev;l++)
+        partner[ipart].ctot[k] += gsl_matrix_get(partner[ipart].colli,k,l);
+    }
+  }
+
+  for(k=0;k<md[ispec].nlev;k++){
+    for(ipart=0;ipart<md[ispec].npart;ipart++){
+      di = md[ispec].part[ipart].densityIndex;
+      if(di>=0)
+        gsl_matrix_set(matrix,k,k,gsl_matrix_get(matrix,k,k)+gp[id].dens[di]*partner[ipart].ctot[k]);
+    }
+    for(l=0;l<md[ispec].nlev;l++){
+      if(k!=l){
+        for(ipart=0;ipart<md[ispec].npart;ipart++){
+          di = md[ispec].part[ipart].densityIndex;
+          if(di>=0)
+            gsl_matrix_set(matrix,k,l,gsl_matrix_get(matrix,k,l)-gp[id].dens[di]*gsl_matrix_get(partner[ipart].colli,l,k));
+        }
+      }
+    }
+    gsl_matrix_set(matrix, md[ispec].nlev, k, 1.);
+    gsl_matrix_set(matrix, k, md[ispec].nlev, 0.);
+  }
+
+  for(ipart=0;ipart<md[ispec].npart;ipart++){
+    gsl_matrix_free(partner[ipart].colli);
+    free(partner[ipart].ctot);
+  }
+  free(partner);
+}
+
+/*....................................................................*/
 void
 stateq(int id, struct grid *gp, molData *md, const int ispec, configInfo *par\
   , struct blendInfo blends, int nextMolWithBlend, gridPointData *mp\
@@ -104,88 +190,4 @@ stateq(int id, struct grid *gp, molData *md, const int ispec, configInfo *par\
   free(oopop);
 }
 
-void
-getmatrix(int id, gsl_matrix *matrix, molData *md, struct grid *gp, int ispec, gridPointData *mp){
-  int ti,k,l,li,ipart,di;
-  struct getmatrix {
-    double *ctot;
-    gsl_matrix * colli;
-  } *partner;
-
-  partner = malloc(sizeof(struct getmatrix)*md[ispec].npart);
-
-  /* Initialize matrix with zeros */
-  for(ipart=0;ipart<md[ispec].npart;ipart++){
-    partner[ipart].colli = gsl_matrix_alloc(md[ispec].nlev+1,md[ispec].nlev+1);
-    if(md[ispec].nlev>0) partner[ipart].ctot  = malloc(sizeof(double)*md[ispec].nlev);
-    else {
-      if(!silent)bail_out("Matrix initialization error in stateq");
-      exit(0);
-    }
-    for(k=0;k<md[ispec].nlev+1;k++){
-      for(l=0;l<md[ispec].nlev+1;l++){
-        gsl_matrix_set(matrix, k, l, 0.);
-        gsl_matrix_set(partner[ipart].colli, k, l, 0.);
-      }
-    }
-  }
-
-  /* Populate matrix with radiative transitions */
-  for(li=0;li<md[ispec].nline;li++){
-    k=md[ispec].lau[li];
-    l=md[ispec].lal[li];
-    gsl_matrix_set(matrix, k, k, gsl_matrix_get(matrix, k, k)+md[ispec].beinstu[li]*mp[ispec].jbar[li]+md[ispec].aeinst[li]);
-    gsl_matrix_set(matrix, l, l, gsl_matrix_get(matrix, l, l)+md[ispec].beinstl[li]*mp[ispec].jbar[li]);
-    gsl_matrix_set(matrix, k, l, gsl_matrix_get(matrix, k, l)-md[ispec].beinstl[li]*mp[ispec].jbar[li]);
-    gsl_matrix_set(matrix, l, k, gsl_matrix_get(matrix, l, k)-md[ispec].beinstu[li]*mp[ispec].jbar[li]-md[ispec].aeinst[li]);
-  }
-
-  /* Populate matrix with collisional transitions */
-  for(ipart=0;ipart<md[ispec].npart;ipart++){
-    struct cpData part = md[ispec].part[ipart];
-    double *downrates = part.down;
-    for(ti=0;ti<part.ntrans;ti++){
-      int coeff_index = ti*part.ntemp + gp[id].mol[ispec].partner[ipart].t_binlow;
-      double down = downrates[coeff_index]\
-                  + gp[id].mol[ispec].partner[ipart].interp_coeff*(downrates[coeff_index+1]\
-                  - downrates[coeff_index]);
-      double up = down*md[ispec].gstat[part.lcu[ti]]/md[ispec].gstat[part.lcl[ti]]\
-                *exp(-HCKB*(md[ispec].eterm[part.lcu[ti]]-md[ispec].eterm[part.lcl[ti]])/gp[id].t[0]);
-      gsl_matrix_set(partner[ipart].colli, part.lcu[ti], part.lcl[ti], down);
-      gsl_matrix_set(partner[ipart].colli, part.lcl[ti], part.lcu[ti], up);
-    }
-
-    for(k=0;k<md[ispec].nlev;k++){
-      partner[ipart].ctot[k]=0.;
-      for(l=0;l<md[ispec].nlev;l++)
-        partner[ipart].ctot[k] += gsl_matrix_get(partner[ipart].colli,k,l);
-    }
-  }
-
-  for(k=0;k<md[ispec].nlev;k++){
-    for(ipart=0;ipart<md[ispec].npart;ipart++){
-      di = md[ispec].part[ipart].densityIndex;
-      if(di>=0)
-        gsl_matrix_set(matrix,k,k,gsl_matrix_get(matrix,k,k)+gp[id].dens[di]*partner[ipart].ctot[k]);
-    }
-    for(l=0;l<md[ispec].nlev;l++){
-      if(k!=l){
-        for(ipart=0;ipart<md[ispec].npart;ipart++){
-          di = md[ispec].part[ipart].densityIndex;
-          if(di>=0)
-            gsl_matrix_set(matrix,k,l,gsl_matrix_get(matrix,k,l)-gp[id].dens[di]*gsl_matrix_get(partner[ipart].colli,l,k));
-        }
-      }
-    }
-    gsl_matrix_set(matrix, md[ispec].nlev, k, 1.);
-    gsl_matrix_set(matrix, k, md[ispec].nlev, 0.);
-  }
-
-  for(ipart=0;ipart<md[ispec].npart;ipart++){
-    gsl_matrix_free(partner[ipart].colli);
-    free(partner[ipart].ctot);
-  }
-  free(partner);
-}
-
 

From 3a7817e725c68e71377686c7731b3e4e9ac3c346 Mon Sep 17 00:00:00 2001
From: ims <ims@strw.leidenuniv.nl>
Date: Thu, 20 Apr 2017 12:17:27 +0200
Subject: [PATCH 2/2] Moved unchanging calculation out of loop.

In stateq.c the part of the matrix involving collision coefficients does not depend on the results of getjbar() and thus has been moved outside the loop now starting at line 141 (old line 45) in stateq.c. This results in a slight speedup.
---
 src/lime.h   |  1 -
 src/stateq.c | 93 +++++++++++++++++++++++++++++++---------------------
 2 files changed, 55 insertions(+), 39 deletions(-)

diff --git a/src/lime.h b/src/lime.h
index d6e3005..4f021ea 100644
--- a/src/lime.h
+++ b/src/lime.h
@@ -337,7 +337,6 @@ void	getclosest(double, double, double, long*, long*, double*, double*, double*)
 double	geterf(const double, const double);
 void	getjbar(int, molData*, struct grid*, const int, configInfo*, struct blendInfo, int, gridPointData*, double*);
 void	getMass(configInfo*, struct grid*, const gsl_rng*);
-//void	getmatrix(int, gsl_matrix*, molData*, struct grid*, int, gridPointData*);
 void	getVelocities(configInfo *, struct grid *);
 void	getVelocities_pregrid(configInfo *, struct grid *);
 void	gridPopsInit(configInfo*, molData*, struct grid*);
diff --git a/src/stateq.c b/src/stateq.c
index a669086..87033b0 100644
--- a/src/stateq.c
+++ b/src/stateq.c
@@ -12,44 +12,34 @@
 #include <gsl/gsl_permutation.h>
 #include <gsl/gsl_errno.h>
 
+struct collPartMatrixType {
+  double *ctot;
+  gsl_matrix *colli;
+};
 
 /*....................................................................*/
 void
-getmatrix(int id, gsl_matrix *matrix, molData *md, struct grid *gp, int ispec, gridPointData *mp){
-  int ti,k,l,li,ipart,di;
-  struct getmatrix {
-    double *ctot;
-    gsl_matrix * colli;
-  } *partner;
+getFixedMatrix(molData *md, int ispec, struct grid *gp, int id, struct collPartMatrixType **collPartMat){
+  int ipart,k,l,ti;
 
-  partner = malloc(sizeof(struct getmatrix)*md[ispec].npart);
+  (*collPartMat) = malloc(sizeof(**collPartMat)*md[ispec].npart);
 
   /* Initialize matrix with zeros */
   for(ipart=0;ipart<md[ispec].npart;ipart++){
-    partner[ipart].colli = gsl_matrix_alloc(md[ispec].nlev+1,md[ispec].nlev+1);
-    if(md[ispec].nlev>0) partner[ipart].ctot  = malloc(sizeof(double)*md[ispec].nlev);
-    else {
-      if(!silent)bail_out("Matrix initialization error in stateq");
-      exit(0);
+    if(md[ispec].nlev<=0){
+      if(!silent) bail_out("Matrix initialization error in stateq");
+      exit(1);
     }
+
+    (*collPartMat)[ipart].colli = gsl_matrix_alloc(md[ispec].nlev+1,md[ispec].nlev+1);
+    (*collPartMat)[ipart].ctot  = malloc(sizeof(double)*md[ispec].nlev);
     for(k=0;k<md[ispec].nlev+1;k++){
       for(l=0;l<md[ispec].nlev+1;l++){
-        gsl_matrix_set(matrix, k, l, 0.);
-        gsl_matrix_set(partner[ipart].colli, k, l, 0.);
+        gsl_matrix_set((*collPartMat)[ipart].colli, k, l, 0.0);
       }
     }
   }
 
-  /* Populate matrix with radiative transitions */
-  for(li=0;li<md[ispec].nline;li++){
-    k=md[ispec].lau[li];
-    l=md[ispec].lal[li];
-    gsl_matrix_set(matrix, k, k, gsl_matrix_get(matrix, k, k)+md[ispec].beinstu[li]*mp[ispec].jbar[li]+md[ispec].aeinst[li]);
-    gsl_matrix_set(matrix, l, l, gsl_matrix_get(matrix, l, l)+md[ispec].beinstl[li]*mp[ispec].jbar[li]);
-    gsl_matrix_set(matrix, k, l, gsl_matrix_get(matrix, k, l)-md[ispec].beinstl[li]*mp[ispec].jbar[li]);
-    gsl_matrix_set(matrix, l, k, gsl_matrix_get(matrix, l, k)-md[ispec].beinstu[li]*mp[ispec].jbar[li]-md[ispec].aeinst[li]);
-  }
-
   /* Populate matrix with collisional transitions */
   for(ipart=0;ipart<md[ispec].npart;ipart++){
     struct cpData part = md[ispec].part[ipart];
@@ -61,41 +51,58 @@ getmatrix(int id, gsl_matrix *matrix, molData *md, struct grid *gp, int ispec, g
                   - downrates[coeff_index]);
       double up = down*md[ispec].gstat[part.lcu[ti]]/md[ispec].gstat[part.lcl[ti]]\
                 *exp(-HCKB*(md[ispec].eterm[part.lcu[ti]]-md[ispec].eterm[part.lcl[ti]])/gp[id].t[0]);
-      gsl_matrix_set(partner[ipart].colli, part.lcu[ti], part.lcl[ti], down);
-      gsl_matrix_set(partner[ipart].colli, part.lcl[ti], part.lcu[ti], up);
+      gsl_matrix_set((*collPartMat)[ipart].colli, part.lcu[ti], part.lcl[ti], down);
+      gsl_matrix_set((*collPartMat)[ipart].colli, part.lcl[ti], part.lcu[ti], up);
     }
 
     for(k=0;k<md[ispec].nlev;k++){
-      partner[ipart].ctot[k]=0.;
+      (*collPartMat)[ipart].ctot[k]=0.0;
       for(l=0;l<md[ispec].nlev;l++)
-        partner[ipart].ctot[k] += gsl_matrix_get(partner[ipart].colli,k,l);
+        (*collPartMat)[ipart].ctot[k] += gsl_matrix_get((*collPartMat)[ipart].colli,k,l);
+    }
+  }
+}
+
+/*....................................................................*/
+void
+getMatrix(int id, gsl_matrix *matrix, molData *md, struct grid *gp, int ispec, gridPointData *mp, struct collPartMatrixType *collPartMat){
+  int k,l,li,ipart,di;
+
+  /* Initialize matrix with zeros */
+  for(k=0;k<md[ispec].nlev+1;k++){
+    for(l=0;l<md[ispec].nlev+1;l++){
+      gsl_matrix_set(matrix, k, l, 0.);
     }
   }
 
+  /* Populate matrix with radiative transitions */
+  for(li=0;li<md[ispec].nline;li++){
+    k=md[ispec].lau[li];
+    l=md[ispec].lal[li];
+    gsl_matrix_set(matrix, k, k, gsl_matrix_get(matrix, k, k)+md[ispec].beinstu[li]*mp[ispec].jbar[li]+md[ispec].aeinst[li]);
+    gsl_matrix_set(matrix, l, l, gsl_matrix_get(matrix, l, l)+md[ispec].beinstl[li]*mp[ispec].jbar[li]);
+    gsl_matrix_set(matrix, k, l, gsl_matrix_get(matrix, k, l)-md[ispec].beinstl[li]*mp[ispec].jbar[li]);
+    gsl_matrix_set(matrix, l, k, gsl_matrix_get(matrix, l, k)-md[ispec].beinstu[li]*mp[ispec].jbar[li]-md[ispec].aeinst[li]);
+  }
+
   for(k=0;k<md[ispec].nlev;k++){
     for(ipart=0;ipart<md[ispec].npart;ipart++){
       di = md[ispec].part[ipart].densityIndex;
       if(di>=0)
-        gsl_matrix_set(matrix,k,k,gsl_matrix_get(matrix,k,k)+gp[id].dens[di]*partner[ipart].ctot[k]);
+        gsl_matrix_set(matrix,k,k,gsl_matrix_get(matrix,k,k)+gp[id].dens[di]*collPartMat[ipart].ctot[k]);
     }
     for(l=0;l<md[ispec].nlev;l++){
       if(k!=l){
         for(ipart=0;ipart<md[ispec].npart;ipart++){
           di = md[ispec].part[ipart].densityIndex;
           if(di>=0)
-            gsl_matrix_set(matrix,k,l,gsl_matrix_get(matrix,k,l)-gp[id].dens[di]*gsl_matrix_get(partner[ipart].colli,l,k));
+            gsl_matrix_set(matrix,k,l,gsl_matrix_get(matrix,k,l)-gp[id].dens[di]*gsl_matrix_get(collPartMat[ipart].colli,l,k));
         }
       }
     }
     gsl_matrix_set(matrix, md[ispec].nlev, k, 1.);
     gsl_matrix_set(matrix, k, md[ispec].nlev, 0.);
   }
-
-  for(ipart=0;ipart<md[ispec].npart;ipart++){
-    gsl_matrix_free(partner[ipart].colli);
-    free(partner[ipart].ctot);
-  }
-  free(partner);
 }
 
 /*....................................................................*/
@@ -104,10 +111,11 @@ stateq(int id, struct grid *gp, molData *md, const int ispec, configInfo *par\
   , struct blendInfo blends, int nextMolWithBlend, gridPointData *mp\
   , double *halfFirstDs, _Bool *luWarningGiven){
 
-  int t,s,iter,status;
+  int t,s,iter,status,ipart;
   double *opop,*oopop,*tempNewPop=NULL;
   double diff;
   char errStr[80];
+  struct collPartMatrixType *collPartMat=NULL;
 
   gsl_matrix *matrix = gsl_matrix_alloc(md[ispec].nlev+1, md[ispec].nlev+1);
   gsl_matrix *reduc  = gsl_matrix_alloc(md[ispec].nlev, md[ispec].nlev);
@@ -128,10 +136,12 @@ stateq(int id, struct grid *gp, molData *md, const int ispec, configInfo *par\
   diff=1;
   iter=0;
 
+  getFixedMatrix(md,ispec,gp,id,&collPartMat);
+
   while((diff>TOL && iter<MAXITER) || iter<5){
     getjbar(id,md,gp,ispec,par,blends,nextMolWithBlend,mp,halfFirstDs);
 
-    getmatrix(id,matrix,md,gp,ispec,mp);
+    getMatrix(id,matrix,md,gp,ispec,mp,collPartMat);
     for(s=0;s<md[ispec].nlev;s++){
       for(t=0;t<md[ispec].nlev-1;t++){
         gsl_matrix_set(reduc,t,s,gsl_matrix_get(matrix,t,s));
@@ -180,6 +190,12 @@ stateq(int id, struct grid *gp, molData *md, const int ispec, configInfo *par\
     iter++;
   }
 
+  for(ipart=0;ipart<md[ispec].npart;ipart++){
+    gsl_matrix_free(collPartMat[ipart].colli);
+    free(collPartMat[ipart].ctot);
+  }
+  free(collPartMat);
+
   gsl_matrix_free(matrix);
   gsl_matrix_free(reduc);
   gsl_vector_free(rhVec);
@@ -191,3 +207,4 @@ stateq(int id, struct grid *gp, molData *md, const int ispec, configInfo *par\
 }
 
 
+