From 29d7dbac9269132e5d183d530b7717646c49708e Mon Sep 17 00:00:00 2001
From: mcodev31 <mcodev31@gmail.com>
Date: Mon, 9 Nov 2015 03:22:15 +0100
Subject: [PATCH] threshold scaling and fix possible null pointer deref

threshold scaling based on found vectors
fix possible null pointer dereference in context
minor fixes
---
 examples/example.c    | 22 ++++++++++--------
 src/context.c         | 53 ++++++++++++++++++++++---------------------
 src/context.h         |  2 +-
 src/equivalence_set.c | 10 +++++++-
 src/linalg.c          |  3 +++
 src/point_group.c     |  2 +-
 src/subspace.c        | 39 ++++++++++++++++---------------
 src/symmetry.c        |  3 +--
 8 files changed, 75 insertions(+), 59 deletions(-)

diff --git a/examples/example.c b/examples/example.c
index 7568a35..89c1aca 100644
--- a/examples/example.c
+++ b/examples/example.c
@@ -70,13 +70,11 @@ int example(const char* in_file, msym_thresholds_t *thresholds){
     int length = read_xyz(in_file, &elements);
     if(length <= 0) return -1;
     
-    //orbitalsl = sizeof(orbitals)/sizeof(char*);
-    //bfsl = orbitalsl*length;
-    //bfs = calloc(bfsl, sizeof(msym_basis_function_t));
-    double (*psalcs)[bfsl] = NULL; //calloc(bfsl, sizeof(*salcs)); // SALCs in matrix form, and input for symmetrization
-    double *pcmem = NULL; // calloc(bfsl, sizeof(double)); // Some temporary memory
-    int *pspecies = NULL; // calloc(bfsl, sizeof(*species));
-    msym_partner_function_t *ppf = NULL; //calloc(bfsl, sizeof(*pf));
+    
+    double (*psalcs)[bfsl] = NULL; // SALCs in matrix form, and input for symmetrization
+    double *pcmem = NULL; // Some temporary memory
+    int *pspecies = NULL;
+    msym_partner_function_t *ppf = NULL;
     
     /* Create a context */
     msym_context ctx = msymCreateContext();
@@ -306,9 +304,13 @@ int example(const char* in_file, msym_thresholds_t *thresholds){
             do{
                 
                 int sel_ss = 0, sel_salc = 0;
-                for(int i = 0; i < msrsl;i++) printf("\t [%d] %s (%d SALCs with %d partner functions each)\n",i, mct->s[msrs[i].s].name, msrs[i].salcl, mct->s[msrs[i].s].d);
-                
-                do {printf("\nChoose subspace [0-%d]:",msrsl-1);} while(scanf(" %d", &sel_ss) <= 0 || sel_ss < 0 || sel_ss >= msrsl);
+                for(int i = 0; i < msrsl;i++){
+                    if(msrs[i].salcl > 0)
+                        printf("\t [%d] %s (%d SALCs with %d partner functions each)\n",i, mct->s[msrs[i].s].name, msrs[i].salcl, mct->s[msrs[i].s].d);
+                    else
+                        printf("\t [-] %s (no SALCs of this symmetry species)\n",mct->s[msrs[i].s].name);
+                }
+                do {printf("\nChoose subspace [0-%d]:",msrsl-1);} while(scanf(" %d", &sel_ss) <= 0 || sel_ss < 0 || sel_ss >= msrsl || msrs[sel_ss].salcl <= 0);
                 
                 int salcl = msrs[sel_ss].salcl;
                 
diff --git a/src/context.c b/src/context.c
index 04ce205..dbb377d 100644
--- a/src/context.c
+++ b/src/context.c
@@ -117,7 +117,7 @@ msym_error_t msymSetThresholds(msym_context ctx, const msym_thresholds_t *thresh
 
 msym_error_t msymGetThresholds(msym_context ctx, const msym_thresholds_t **thresholds){
     msym_error_t ret = MSYM_SUCCESS;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     if(ctx->thresholds == NULL) ret = MSYM_INVALID_THRESHOLD;
     *thresholds = ctx->thresholds;
 err:
@@ -126,7 +126,7 @@ msym_error_t msymGetThresholds(msym_context ctx, const msym_thresholds_t **thres
 
 msym_error_t msymSetElements(msym_context ctx, int length, msym_element_t *elements){
     msym_error_t ret = MSYM_SUCCESS;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     /* Allow manual setting of point group before elements */
     if(NULL != ctx->es) ctxDestroyPointGroup(ctx);
     return ctxSetElements(ctx, length, elements);
@@ -137,7 +137,7 @@ msym_error_t msymSetElements(msym_context ctx, int length, msym_element_t *eleme
 msym_error_t msymGetElements(msym_context ctx, int *length, msym_element_t **elements){
     msym_error_t ret = MSYM_SUCCESS;
     msym_element_t *relements = NULL;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     if(ctx->elements == NULL || ctx->ext.elements == NULL) {ret = MSYM_INVALID_ELEMENTS;goto err;}
     
     *elements = ctx->ext.elements;
@@ -164,7 +164,7 @@ msym_error_t msymGetEquivalenceSets(msym_context ctx, int *length, const msym_eq
 
 msym_error_t msymGetBasisFunctions(msym_context ctx, int *length, msym_basis_function_t **basis){
     msym_error_t ret = MSYM_SUCCESS;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     if(ctx->basis == NULL) {
         msymSetErrorDetails("Found no basis functions");
         ret = MSYM_INVALID_BASIS_FUNCTIONS;
@@ -179,7 +179,7 @@ msym_error_t msymGetBasisFunctions(msym_context ctx, int *length, msym_basis_fun
 
 msym_error_t msymSetBasisFunctions(msym_context ctx, int length, msym_basis_function_t *basis){
     msym_error_t ret = MSYM_SUCCESS;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     if(ctx->elements == NULL) {ret = MSYM_INVALID_ELEMENTS;goto err;}
     ctxDestroyBasisFunctions(ctx);
     ctx->basis = malloc(sizeof(msym_basis_function_t[length]));
@@ -229,7 +229,7 @@ msym_error_t msymSetBasisFunctions(msym_context ctx, int length, msym_basis_func
 
 msym_error_t msymGetPointGroupName(msym_context ctx, int l, char buf[l]){
     msym_error_t ret = MSYM_SUCCESS;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     if(ctx->pg == NULL) {ret = MSYM_INVALID_POINT_GROUP;goto err;}
     snprintf(buf, l, "%s",ctx->pg->name);
 err:
@@ -238,7 +238,7 @@ msym_error_t msymGetPointGroupName(msym_context ctx, int l, char buf[l]){
 
 msym_error_t msymGetPointGroupType(msym_context ctx, msym_point_group_type_t *t, int *n){
     msym_error_t ret = MSYM_SUCCESS;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     if(ctx->pg == NULL) {ret = MSYM_INVALID_POINT_GROUP;goto err;}
     
     *t = ctx->pg->type;
@@ -254,7 +254,7 @@ msym_error_t msymGetSubgroups(msym_context ctx, int *sgl, const msym_subgroup_t
     msym_subgroup_t *gsg = NULL;
     int gsgl = 0;
     
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     if(ctx->pg == NULL) {ret = MSYM_INVALID_POINT_GROUP;goto err;}
     if(ctx->pg->perm == NULL && !(isLinearPointGroup(ctx->pg) && !isLinearSubgroup(ctx->pg))) {
         ret = MSYM_INVALID_PERMUTATION;
@@ -337,7 +337,7 @@ msym_error_t msymGetCharacterTable(msym_context ctx, const msym_character_table_
 
 msym_error_t msymGetCenterOfMass(msym_context ctx, double v[3]){
     msym_error_t ret = MSYM_SUCCESS;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     if(ctx->elements == NULL) {ret = MSYM_INVALID_ELEMENTS;goto err;}
     vcopy(ctx->cm, v);
 err:
@@ -355,7 +355,7 @@ msym_error_t msymSetCenterOfMass(msym_context ctx, double cm[3]){
 
 msym_error_t msymGetGeometry(msym_context ctx, msym_geometry_t *geometry){
     msym_error_t ret = MSYM_SUCCESS;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     if(ctx->elements == NULL) {ret = MSYM_INVALID_ELEMENTS;goto err;}
     if(ctx->geometry == MSYM_GEOMETRY_UNKNOWN) {ret = MSYM_INVALID_GEOMETRY;goto err;}
     *geometry = ctx->geometry;
@@ -365,7 +365,7 @@ msym_error_t msymGetGeometry(msym_context ctx, msym_geometry_t *geometry){
 
 msym_error_t msymGetPrincipalMoments(msym_context ctx, double eigval[3]){
     msym_error_t ret = MSYM_SUCCESS;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     if(ctx->elements == NULL) {ret = MSYM_INVALID_ELEMENTS;goto err;}
     vcopy(ctx->eigval, eigval);
 err:
@@ -373,7 +373,7 @@ msym_error_t msymGetPrincipalMoments(msym_context ctx, double eigval[3]){
 }
 msym_error_t msymGetPrincipalAxes(msym_context ctx, double eigvec[3][3]){
     msym_error_t ret = MSYM_SUCCESS;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     if(ctx->elements == NULL) {ret = MSYM_INVALID_ELEMENTS;goto err;}
     mcopy(ctx->eigvec, eigvec);
 err:
@@ -382,9 +382,9 @@ msym_error_t msymGetPrincipalAxes(msym_context ctx, double eigvec[3][3]){
 
 msym_error_t msymGetRadius(msym_context ctx, double *radius){
     msym_error_t ret = MSYM_SUCCESS;
-    double r = 0.0;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     if(ctx->elements == NULL) {ret = MSYM_INVALID_ELEMENTS;goto err;}
+    double r = 0.0;
     for(int i = 0;i < ctx->elementsl;i++){
         double abs = vabs(ctx->elements[i].v);
         r = r > abs ? r : abs;
@@ -398,7 +398,7 @@ msym_error_t msymGetRadius(msym_context ctx, double *radius){
 msym_error_t msymGetSymmetryOperations(msym_context ctx, int *sopsl, const msym_symmetry_operation_t **sops){
     msym_error_t ret = MSYM_SUCCESS;
     msym_symmetry_operation_t *rsops = NULL;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     if(ctx->pg == NULL || ctx->pg->sops == NULL) {ret = MSYM_INVALID_POINT_GROUP;goto err;}
     
     *sops = ctx->pg->sops;
@@ -413,12 +413,12 @@ msym_error_t msymGetSymmetryOperations(msym_context ctx, int *sopsl, const msym_
 
 msym_error_t msymReleaseContext(msym_context ctx){
     msym_error_t ret = MSYM_SUCCESS;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     free(ctx->thresholds);
     ctxDestroyElements(ctx);
     ctxDestroyPointGroup(ctx);
     free(ctx);
-err:
+//err:
     return ret;
 }
 
@@ -433,7 +433,7 @@ msym_error_t msymReleaseContext(msym_context ctx){
 msym_error_t ctxSetElements(msym_context ctx, int length, msym_element_t elements[length]){
     msym_error_t ret = MSYM_SUCCESS;
     msym_thresholds_t *thresholds = NULL;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     
     ctxDestroyElements(ctx);
     
@@ -466,6 +466,7 @@ msym_error_t ctxSetElements(msym_context ctx, int length, msym_element_t element
     
     return ret;
 err:
+    
     free(ctx->elements);
     free(ctx->pelements);
     free(ctx->ext.elements);
@@ -480,8 +481,8 @@ msym_error_t ctxSetElements(msym_context ctx, int length, msym_element_t element
 
 msym_error_t ctxGetThresholds(msym_context ctx, msym_thresholds_t **thresholds){
     msym_error_t ret = MSYM_SUCCESS;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
-    if(ctx->thresholds == NULL) ret = MSYM_INVALID_THRESHOLD;
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
+    if(ctx->thresholds == NULL) {ret = MSYM_INVALID_THRESHOLD; goto err;}
     msym_thresholds_t *t = ctx->thresholds;
     if(t->angle < 1.0 && !signbit(t->angle) &&
        t->equivalence < 1.0 && !signbit(t->equivalence) &&
@@ -501,7 +502,7 @@ msym_error_t ctxGetThresholds(msym_context ctx, msym_thresholds_t **thresholds){
 
 msym_error_t ctxGetElements(msym_context ctx, int *l, msym_element_t **elements){
     msym_error_t ret = MSYM_SUCCESS;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     if(ctx->elements == NULL) {ret = MSYM_INVALID_ELEMENTS; goto err;}
     *elements = (msym_element_t *) ctx->elements;
     *l = ctx->elementsl;
@@ -511,7 +512,7 @@ msym_error_t ctxGetElements(msym_context ctx, int *l, msym_element_t **elements)
 
 msym_error_t ctxGetExternalElements(msym_context ctx, int *l, msym_element_t **elements){
     msym_error_t ret = MSYM_SUCCESS;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     if(ctx->ext.elements == NULL) {ret = MSYM_INVALID_ELEMENTS; goto err;}
     *elements = (msym_element_t *) ctx->ext.elements;
     *l = ctx->elementsl;
@@ -537,7 +538,7 @@ msym_error_t ctxUpdateExternalElementCoordinates(msym_context ctx){
 
 msym_error_t ctxGetInternalElement(msym_context ctx, msym_element_t *ext, msym_element_t **element){
     msym_error_t ret = MSYM_SUCCESS;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     if(ctx->ext.elements == NULL) {ret = MSYM_INVALID_ELEMENTS;goto err;}
     if(ext < ctx->ext.elements || ext >= ctx->ext.elements+ctx->elementsl){
         msymSetErrorDetails("Element pointer (%p) outside memory block (%p -> %p)", ext, ctx->ext.elements, ctx->ext.elements + ctx->elementsl);
@@ -551,7 +552,7 @@ msym_error_t ctxGetInternalElement(msym_context ctx, msym_element_t *ext, msym_e
 
 msym_error_t ctxGetSubgroups(msym_context ctx, int *sgl, msym_subgroup_t **sg){
     msym_error_t ret = MSYM_SUCCESS;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     if(ctx->sg == NULL) {ret = MSYM_INVALID_SUBGROUPS;goto err;}
     *sg = ctx->sg;
     *sgl = ctx->sgl;
@@ -561,7 +562,7 @@ msym_error_t ctxGetSubgroups(msym_context ctx, int *sgl, msym_subgroup_t **sg){
 
 msym_error_t ctxSetSubgroups(msym_context ctx, int sgl, msym_subgroup_t *sg){
     msym_error_t ret = MSYM_SUCCESS;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     ctxDestroySubgroups(ctx);
     ctx->sg = sg;
     ctx->sgl = sgl;
@@ -572,7 +573,7 @@ msym_error_t ctxSetSubgroups(msym_context ctx, int sgl, msym_subgroup_t *sg){
 
 msym_error_t ctxGetElementPtrs(msym_context ctx, int *l, msym_element_t ***pelements){
     msym_error_t ret = MSYM_SUCCESS;
-    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT;goto err;}
+    if(ctx == NULL) {ret = MSYM_INVALID_CONTEXT; return ret;}
     if(ctx->pelements == NULL) {ret = MSYM_INVALID_ELEMENTS; goto err;}
     *pelements = (msym_element_t **) ctx->pelements;
     *l = ctx->elementsl;
diff --git a/src/context.h b/src/context.h
index 6531644..d967883 100644
--- a/src/context.h
+++ b/src/context.h
@@ -22,7 +22,7 @@
 #define DEFAULT_EQUIVALENCE_THRESHOLD 5.0e-4
 #define DEFAULT_EIGFACT_THRESHOLD 1.0e-3
 #define DEFAULT_PERMUTATION_THRESHOLD 5.0e-3
-#define DEFAULT_ORTHOGONALIZATION_THRESHOLD 0.01
+#define DEFAULT_ORTHOGONALIZATION_THRESHOLD 1.0e-2
 
 
 msym_error_t ctxGetThresholds(msym_context ctx, msym_thresholds_t **thresholds);
diff --git a/src/equivalence_set.c b/src/equivalence_set.c
index 95b8ab9..293ed30 100644
--- a/src/equivalence_set.c
+++ b/src/equivalence_set.c
@@ -48,6 +48,13 @@ msym_error_t generateEquivalenceSet(msym_point_group_t *pg, int length, msym_ele
     msym_equivalence_set_t *ges = calloc(length,sizeof(msym_equivalence_set_t));
     int gel = 0;
     int gesl = 0;
+    
+    if(pg->order <= 0){
+        ret = MSYM_INVALID_POINT_GROUP;
+        msymSetErrorDetails("Point group of zero order when determining equivalence set");
+        goto err;
+    }
+
     for(int i = 0;i < length;i++){
         double ev[3];
         vsub(elements[i].v, cm, ev);
@@ -83,7 +90,7 @@ msym_error_t generateEquivalenceSet(msym_point_group_t *pg, int length, msym_ele
             }
         }
         
-        if(pg->order % aes->length != 0){
+        if(!aes->length || (pg->order % aes->length != 0)){
             msymSetErrorDetails("Equivalence set length (%d) not a divisor of point group order (%d)",pg->order);
             ret = MSYM_INVALID_EQUIVALENCE_SET;
             goto err;
@@ -185,6 +192,7 @@ msym_error_t findPointGroupEquivalenceSets(msym_point_group_t *pg, int length, m
     free(pelements);
     return ret;
 err:
+    free(ges);
     free(pelements);
     return ret;
 
diff --git a/src/linalg.c b/src/linalg.c
index cb83cfc..2b53449 100644
--- a/src/linalg.c
+++ b/src/linalg.c
@@ -535,6 +535,7 @@ int mgs(int l, const double m[l][l], double o[l][l], int n, double t){
 int mgs2(int l, int lm, const double m[l][l], double o[l][l], int n, double t){
     
     int nm = n + lm + MGS_ADD;
+    double ts = l/(1.0 + l);
     for(int i = 0; i < l && n < nm;i++){
         if(vlabs(l, m[i]) < t){
             continue;
@@ -550,8 +551,10 @@ int mgs2(int l, int lm, const double m[l][l], double o[l][l], int n, double t){
                 for(int k = 0;k < l;k++) o[n][k] -= c*o[j][k];
             }
             if(vlabs(l, o[n]) >= t){
+                t *= ts;
                 vlnorm(l, o[n]);
                 n++;
+                
             }
         }
     }
diff --git a/src/point_group.c b/src/point_group.c
index 3728f7f..4f2baea 100644
--- a/src/point_group.c
+++ b/src/point_group.c
@@ -793,8 +793,8 @@ msym_error_t transformAxes(msym_point_group_type_t type, int n, msym_symmetry_op
         case (MSYM_POINT_GROUP_TYPE_Ih)  : {
             msym_symmetry_operation_t *dprimary = NULL;
             msym_symmetry_operation_t *sop;
+            double z = -2.0;
             for(sop = sops; sop < (sops + sopsl); sop++){
-                double z = -2.0;
                 if(sop->type == PROPER_ROTATION && sop->order == 2){
                     double v[3];
                     vnorm2(sop->v,v);
diff --git a/src/subspace.c b/src/subspace.c
index 583b04a..88794a2 100644
--- a/src/subspace.c
+++ b/src/subspace.c
@@ -140,7 +140,7 @@ msym_error_t generatePermutationSubspaces(msym_point_group_t *pg, msym_permutati
             }
         }
         
-        nirl = mgs2(dim, vspan,proj, ss, oirl, thresholds->orthogonalization/dim);
+        nirl = mgs2(dim, vspan,proj, ss, oirl, thresholds->orthogonalization);
         if(nirl - oirl != vspan){
             debug_printTransform(dim, dim, ss);
             ret = MSYM_SUBSPACE_ERROR;
@@ -183,7 +183,7 @@ msym_error_t generateSubspaces(msym_point_group_t *pg, msym_permutation_t perm[p
         if(MSYM_SUCCESS != (ret = generateProjectionOperator(irrepd,sopsl,cmem,perm,ld,lrsops,proj))) goto err;
         
         if(irrepd == 1){
-            nirl = mgs2(dim, pgvspan, proj, ss, oirl, thresholds->orthogonalization/dim);
+            nirl = mgs2(dim, pgvspan, proj, ss, oirl, thresholds->orthogonalization);
             if(nirl - oirl != pgvspan){
                 debug_printTransform(dim, dim, ss);
                 ret = MSYM_SUBSPACE_ERROR;
@@ -193,12 +193,12 @@ msym_error_t generateSubspaces(msym_point_group_t *pg, msym_permutation_t perm[p
             pss[k][0] = &ss[oirl];
         } else if(!(icosahedral && irrepd == 5)){
             
-            int pgnirl = mgs2(dim, pgvspan, proj, sspg, 0, thresholds->orthogonalization/dim);
+            int pgnirl = mgs2(dim, pgvspan, proj, sspg, 0, thresholds->orthogonalization);
             for(int d = 0; d < irrepd;d++,oirl = nirl){
                 
                 if(MSYM_SUCCESS != (ret = generateProjectionOperator(1,sopsl,sgc[k][d],perm,ld,lrsops,proj))) goto err;
                 
-                int sgnirl = mgs2(dim, sgd[k][d], proj, sssg, 0, thresholds->orthogonalization/dim);
+                int sgnirl = mgs2(dim, sgd[k][d], proj, sssg, 0, thresholds->orthogonalization);
                 
                 if(MSYM_SUCCESS != (ret = projectLinearlyIndependent(dim, pgnirl, sspg, sgnirl, sssg, thresholds, cmem, mem, ss, &nirl))) goto err;
                 
@@ -214,11 +214,11 @@ msym_error_t generateSubspaces(msym_point_group_t *pg, msym_permutation_t perm[p
             }
         } else {
             int idim[] = {1,2,2}, sdim[] = {3,4}, ssd = 0;
-            int pgnirl = mgs2(dim, pgvspan, proj, sspg, 0, thresholds->orthogonalization/dim);
+            int pgnirl = mgs2(dim, pgvspan, proj, sspg, 0, thresholds->orthogonalization);
             for(int d = 0; d < 3;d++,oirl = nirl){
                 if(MSYM_SUCCESS != (ret = generateProjectionOperator(idim[d],sopsl,sgc[k][d],perm,ld,lrsops,proj))) goto err;
                 
-                int sgnirl = mgs2(dim, sgd[k][d], proj, sssg, 0, thresholds->orthogonalization/dim);
+                int sgnirl = mgs2(dim, sgd[k][d], proj, sssg, 0, thresholds->orthogonalization);
                 
                 int id = idim[d];
                 if(id > 1){
@@ -230,7 +230,7 @@ msym_error_t generateSubspaces(msym_point_group_t *pg, msym_permutation_t perm[p
                         if(MSYM_SUCCESS != (ret = generateProjectionOperator(1,sopsl,sgc[k][sid],perm,ld,lrsops,proj))) goto err;
                         
                         // sspg, sssg and proj are taken, use mem, and take proj as mem after
-                        int ignirl = mgs2(dim, sgd[k][sid], proj, mem, 0, thresholds->orthogonalization/dim);
+                        int ignirl = mgs2(dim, sgd[k][sid], proj, mem, 0, thresholds->orthogonalization);
                         
                         if(MSYM_SUCCESS != (ret = projectLinearlyIndependent(dim, sgnirl, sssg, ignirl, mem, thresholds, cmem, proj, ss, &nirl))) goto err;
                         
@@ -291,7 +291,7 @@ msym_error_t generateSubspacesMatrix(msym_point_group_t *pg, msym_permutation_t
         if(MSYM_SUCCESS != (ret = generateProjectionOperator(irrepd,sopsl,cmem,perm,ld,lrsops,projpg))) goto err;
         
         if(irrepd == 1){
-            nirl = mgs2(dim, pgvspan, projpg, ss, oirl, thresholds->orthogonalization/dim);
+            nirl = mgs2(dim, pgvspan, projpg, ss, oirl, thresholds->orthogonalization/sqrt(dim));
             if(nirl - oirl != pgvspan){
                 debug_printTransform(dim, dim, ss);
                 ret = MSYM_SUBSPACE_ERROR;
@@ -310,7 +310,7 @@ msym_error_t generateSubspacesMatrix(msym_point_group_t *pg, msym_permutation_t
                 trace = mltrace(dim, mem);
                 mlscale(span[k]/trace, dim, mem, mem);
                 
-                nirl = mgs2(dim, span[k], mem, ss, oirl, thresholds->orthogonalization/dim);
+                nirl = mgs2(dim, span[k], mem, ss, oirl, thresholds->orthogonalization/sqrt(dim));
                 if(nirl - oirl != span[k]){
                     debug_printTransform(dim, dim, ss);
                     ret = MSYM_SUBSPACE_ERROR;
@@ -340,7 +340,7 @@ msym_error_t generateSubspacesMatrix(msym_point_group_t *pg, msym_permutation_t
                         trace = mltrace(dim, mem);
                         mlscale(span[k]/trace, dim, mem, mem); // We might have small components in these subspaces
                         
-                        nirl = mgs2(dim, span[k], mem, ss, oirl, thresholds->orthogonalization/dim);
+                        nirl = mgs2(dim, span[k], mem, ss, oirl, thresholds->orthogonalization/sqrt(dim));
                         if(nirl - oirl != span[k]){
                             debug_printTransform(dim, dim, mem);
                             ret = MSYM_SUBSPACE_ERROR;
@@ -352,7 +352,7 @@ msym_error_t generateSubspacesMatrix(msym_point_group_t *pg, msym_permutation_t
                     }
                 } else {
                     
-                    nirl = mgs2(dim, span[k], projig, ss, oirl, thresholds->orthogonalization/dim);
+                    nirl = mgs2(dim, span[k], projig, ss, oirl, thresholds->orthogonalization/sqrt(dim));
                     if(nirl - oirl != span[k]){
                         debug_printTransform(dim, dim, ss);
                         ret = MSYM_SUBSPACE_ERROR;
@@ -593,7 +593,6 @@ msym_error_t getSplittingFieldCharacters(msym_point_group_t *pg, const msym_subg
             for(int i = 0;i < sg->order;i++){
                 if(&pg->sops[s] != sg->sops[i]) continue;
                 if(pg->sops[s].type == IDENTITY){
-                    e = 1;
                     c[0][s] = c[3][s] = c[4][s] = 1;
                     c[1][s] = c[2][s] = 2;
                     found++;
@@ -638,7 +637,7 @@ msym_error_t projectLinearlyIndependent(int dim, int vdim, double v[vdim][dim],
     
     int mdim = vdim > udim ? udim : vdim;
 
-    int nirl = mgs2(dim, mdim, mem, o, *oirl, thresholds->orthogonalization/dim);
+    int nirl = mgs2(dim, mdim, mem, o, *oirl, thresholds->orthogonalization/sqrt(dim));
     
     for(int i = *oirl; i < nirl;i++) vlnorm(dim, o[i]);
     
@@ -956,7 +955,6 @@ msym_error_t generateSubrepresentationSpaces(msym_point_group_t *pg, int sgl, co
     const msym_subgroup_t **rsg = calloc(ct->d, sizeof(*rsg));
     
     int (*sgd)[5] = calloc(ct->d,sizeof(*sgd));
-    int (*sgdASDASDASD)[5] = calloc(ct->d,sizeof(*sgdASDASDASD));
     int *ispan = calloc(ct->d, sizeof(*ispan));                               // decoposed total span of symmetrized basis (int)
     int (*iespan)[lmax+1][ct->d] = calloc(esl, sizeof(*iespan));
     int (*ipspan)[ct->d] = calloc(esl, sizeof(*ipspan));               // span of permutation operators
@@ -991,7 +989,7 @@ msym_error_t generateSubrepresentationSpaces(msym_point_group_t *pg, int sgl, co
     for(int k = 0; k < ct->d;k++){
         if(ct->s[k].d > 1){
             if(MSYM_SUCCESS != (ret = findSplittingFieldSubgroup(pg, k, sgl, sg, thresholds, &rsg[k]))) goto err;
-            if(MSYM_SUCCESS != (ret = getSplittingFieldCharacters(pg, rsg[k], sgc[k], sgdASDASDASD[k]))) goto err; //TODO: remove sgd
+            if(MSYM_SUCCESS != (ret = getSplittingFieldCharacters(pg, rsg[k], sgc[k], sgd[k]))) goto err; //TODO: remove sgd
         }
     }
     
@@ -1098,7 +1096,7 @@ msym_error_t generateSubrepresentationSpaces(msym_point_group_t *pg, int sgl, co
                 }
             }
             
-            nirl = mgs2(d, vspan,lproj, st[pg->order], oirl, thresholds->orthogonalization/basisl);
+            nirl = mgs2(d, vspan,lproj, st[pg->order], oirl, thresholds->orthogonalization);
             
             if(nirl - oirl != vspan){
                 debug_printTransform(d, d, st[pg->order]);
@@ -1226,7 +1224,7 @@ msym_error_t generateSubrepresentationSpaces(msym_point_group_t *pg, int sgl, co
                             if(sdvl != sspan){
                                 debug_printTransform(vspan, dim, dss);
                                 ret = MSYM_SUBSPACE_ERROR;
-                                msymSetErrorDetails("Linear projection subspace of dimension (%d) inconsistent with span (%d) in %s component %d",sdvl,sd,ct->s[dk].name,d);
+                                msymSetErrorDetails("Linear projection subspace of dimension (%d) inconsistent with span (%d) in %s component %d",sdvl,sspan,ct->s[dk].name,d);
                                 goto err;
                             }
                             
@@ -1284,12 +1282,12 @@ msym_error_t generateSubrepresentationSpaces(msym_point_group_t *pg, int sgl, co
     free(iespan);
     free(ipspan);
     free(ibspan);
+    free(sspmem);
     free(psspmem);
     free(lssp);
     free(rsg);
     free(sgc);
     free(sgd);
-    free(sgdASDASDASD);
     free(isalc);
     free(esnmax);
     free(esbfmap);
@@ -1314,10 +1312,12 @@ msym_error_t generateSubrepresentationSpaces(msym_point_group_t *pg, int sgl, co
     free(iespan);
     free(ipspan);
     free(ibspan);
+    free(sspmem);
     free(psspmem);
     free(lssp);
     free(rsg);
     free(sgc);
+    free(sgd);
     free(ispan);
     free(isalc);
     free(esnmax);
@@ -1824,6 +1824,7 @@ msym_error_t generateSubrepresentationSpacesLowMem(msym_point_group_t *pg, int s
     free(rsg);
     free(sgc);
     free(sgd);
+    free(mpih);
     free(mdcomp);
     free(mdproj);
     free(mdfound);
@@ -1854,6 +1855,8 @@ msym_error_t generateSubrepresentationSpacesLowMem(msym_point_group_t *pg, int s
     free(mdec);
     free(rsg);
     free(sgc);
+    free(sgd);
+    free(mpih);
     free(mdcomp);
     free(mdproj);
     free(mdfound);
diff --git a/src/symmetry.c b/src/symmetry.c
index 9619a2c..d5c29e8 100644
--- a/src/symmetry.c
+++ b/src/symmetry.c
@@ -673,7 +673,7 @@ msym_error_t findSymmetryCubic(msym_equivalence_set_t *es, double cm[3], double
         double d;
     } *pairs = malloc(sizeof(struct _pair[150]));*/
     double c2d = 0, c4d = 0, sigmad = 0;
-    int found = 0, nsigma = 0, esigma = 0, inversion = 0, nc[6] = {0,0,0,0,0,0}, ec[6] = {0,0,0,0,0,0}, *ncb = &nsigma, *ecb = &esigma,*nc3b = &nsigma;
+    int found = 0, nsigma = 0, esigma = 0, inversion = 0, nc[6] = {0,0,0,0,0,0}, ec[6] = {0,0,0,0,0,0}, *ncb = &nsigma, *nc3b = &nsigma;
     msym_symmetry_operation_t **(ac[6]);
     double thetac[6] = {0.0,M_PI,M_PI/2,M_PI/3,M_PI/4,M_PI/5};
     
@@ -901,7 +901,6 @@ msym_error_t findSymmetryCubic(msym_equivalence_set_t *es, double cm[3], double
                 case 15 :
                     ec[5] = 6; ec[3] = 10; ec[2] = 15;
                     ncb = &(nc[2]);
-                    ecb = &(ec[2]);
                     cb = ac[2];
                     break;
                 default :