diff --git a/src/DistributeMat.c b/src/DistributeMat.c
index ea4e889eefa213d91bd233ce8722cbd5a690a480..bb0c42868c878015561541d8e95d86ab4f5311fd 100644
--- a/src/DistributeMat.c
+++ b/src/DistributeMat.c
@@ -1,5 +1,6 @@
 #include "DistributeMat.h"
 #include "Permute.h"
+#include "FreeNonzeros.h"
 
 #ifdef USE_PATOH
 #include <patoh.h>
@@ -511,6 +512,43 @@ int DistributeMatrixMondriaan(struct sparsematrix *pT, int P, double eps, const
             fprintf(stderr, "DistributeMatrixMondriaan(): Unknown SplitMethod!\n");
             return FALSE;
         }
+        
+        /* Shift weight and procs */
+        for (j = k; j > i+1; j--) {
+            weight[j] = weight[j-1];
+            procs[j] = procs[j-1];
+        }
+
+        k++; /* new number of parts */
+
+        weight[i] = ComputeWeight(pT, pT->Pstart[i], pT->Pstart[i+1]-1, NULL, pOptions);
+        weight[i+1] = ComputeWeight(pT, pT->Pstart[i+1], pT->Pstart[i+2]-1, NULL, pOptions);
+        
+        if (weight[i] < 0 || weight[i + 1] < 0) {
+            fprintf(stderr, "DistributeMatrixMondriaan(): Unable to compute weights!\n");
+            return FALSE;
+        }
+        
+        procslo = procs[i]/2;
+        procshi = (procs[i]%2==0 ? procslo : procslo+1);
+        
+        if (weight[i] <= weight[i+1]) {
+            procs[i] = procslo;
+            procs[i+1] = procshi;
+        } else {
+            procs[i] = procshi;
+            procs[i+1] = procslo;
+        }
+        
+        /* Apply free nonzero search if enabled, but only for (symmetric) finegrain and mediumgrain strategies */
+        if(pOptions->ImproveFreeNonzeros == FreeNonzerosYes && (pOptions->SplitStrategy == FineGrain ||
+            pOptions->SplitStrategy == SFineGrain || pOptions->SplitStrategy == MediumGrain)) {
+            ImproveFreeNonzeros(pT, pOptions, procs, i, i+1);
+            
+            weight[i] = ComputeWeight(pT, pT->Pstart[i], pT->Pstart[i+1]-1, NULL, pOptions);
+            weight[i+1] = ComputeWeight(pT, pT->Pstart[i+1], pT->Pstart[i+2]-1, NULL, pOptions);
+        }
+        
 #ifdef INFO2
         printf("  Pstart[%d] = %ld ", i,  pT->Pstart[i]);
         printf("Pstart[%d] = %ld ", i+1,  pT->Pstart[i+1]);
@@ -592,34 +630,6 @@ int DistributeMatrixMondriaan(struct sparsematrix *pT, int P, double eps, const
                 }
             }
         }
-        
-        /* Shift weight and procs */
-        for (j = k; j > i+1; j--) {
-            weight[j] = weight[j-1];
-            procs[j] = procs[j-1];
-        }
-
-        k++; /* new number of parts */
-
-        weight[i] = ComputeWeight(pT, pT->Pstart[i], pT->Pstart[i+1]-1, NULL, pOptions);
-        weight[i+1] = ComputeWeight(pT, pT->Pstart[i+1], pT->Pstart[i+2]-1, NULL, pOptions);
-        
-        if (weight[i] < 0 || weight[i + 1] < 0) {
-            fprintf(stderr, "DistributeMatrixMondriaan(): Unable to compute weights!\n");
-            return FALSE;
-        }
-        
-        procslo = procs[i]/2;
-        procshi = (procs[i]%2==0 ? procslo : procslo+1);
-        
-        if (weight[i] <= weight[i+1]) {
-            procs[i] = procslo;
-            procs[i+1] = procshi;
-        } else { 
-            procs[i] = procshi;
-            procs[i+1] = procslo;
-        }
-
 
         /* Check if there is a part that is too large */
         done = TRUE;
diff --git a/src/FreeNonzeros.c b/src/FreeNonzeros.c
new file mode 100644
index 0000000000000000000000000000000000000000..544939d7c0d769803c2a3885eea7bd006003e9c5
--- /dev/null
+++ b/src/FreeNonzeros.c
@@ -0,0 +1,127 @@
+#include "FreeNonzeros.h"
+
+/**
+ * Swap nonzeros s and t of matrix pM
+ */
+void SwapNonzero(struct sparsematrix *pM, long s, long t) {
+	SwapLong(pM->i, s, t);
+	SwapLong(pM->j, s, t);
+	if(pM->MMTypeCode[2] != 'P')
+		SwapDouble(pM->ReValue, s, t);
+	if(pM->MMTypeCode[2] == 'C')
+		SwapDouble(pM->ImValue, s, t);
+} /* end SwapNonzero */
+
+/**
+ * Improve load balance by moving free nonzeros between two specified processors.
+ * This function works for general and symmetric matrices.
+ * It works when dummy nonzeros are present: dummies will not be moved as they do not contribute weight.
+ * It does not work when weights are based on column weights.
+ * 
+ * Input:
+ *   pOptions    Options struct
+ *   procs       abs(procs[i]) = Number of processors each part should still be distributed over
+ *   p1, p2      Processors to condider
+ * 
+ * Input/output:
+ *   pM          The matrix
+ * 
+ * Return value: FALSE if an error occurred, TRUE otherwise
+ */
+
+int ImproveFreeNonzeros(struct sparsematrix *pM, const struct opts *pOptions, const int *procs, int p1, int p2) {
+	
+	long nnz1 = ComputeWeight(pM, pM->Pstart[p1], pM->Pstart[p1+1]-1, NULL, pOptions);
+	long nnz2 = ComputeWeight(pM, pM->Pstart[p2], pM->Pstart[p2+1]-1, NULL, pOptions);
+	long i, j, t, nzWeight;
+	
+	int symmetric = pM->m == pM->n &&
+		(pM->MMTypeCode[3]=='S' || pM->MMTypeCode[3]=='K' || pM->MMTypeCode[3]=='H') &&
+		pOptions->SymmetricMatrix_UseSingleEntry == SingleEntYes;
+		
+	int dummies = pM->m == pM->n && pM->NrDummies > 0;
+	
+	if(nnz1/abs(procs[p1]) > nnz2/abs(procs[p2])) {
+		long tmp = p1;
+		p1 = p2;
+		p2 = tmp;
+		tmp = nnz1;
+		nnz1 = nnz2;
+		nnz2 = tmp;
+	}
+	
+	/* Now p2 is relatively larger than p1 */
+	if((nnz1+1)/(double)abs(procs[p1]) >= (nnz2-1)/(double)abs(procs[p2])) {
+		return TRUE;
+	}
+	
+	char *rowsP1 = (char *)malloc(pM->m * sizeof(char));
+	char *colsP1 = (char *)malloc(pM->n * sizeof(char));
+	
+	if(rowsP1 == NULL || colsP1 == NULL) {
+		fprintf(stderr, "ImproveFreeNonzeros(): Not enough memory!");
+		if(rowsP1 != NULL)
+			free(rowsP1);
+		if(colsP1 != NULL)
+			free(colsP1);
+		return FALSE;
+	}
+	
+	for(i=0; i<pM->m; ++i) {
+		rowsP1[i] = 0;
+	}
+	for(j=0; j<pM->n; ++j) {
+		colsP1[j] = 0;
+	}
+	
+	/* In what columns/rows does P(1) have nonzeros? */
+	for(t=pM->Pstart[p1]; t<pM->Pstart[p1+1]; ++t) {
+		rowsP1[pM->i[t]] = 1;
+		colsP1[pM->j[t]] = 1;
+		
+		if(symmetric) {
+			rowsP1[pM->j[t]] = 1;
+			colsP1[pM->i[t]] = 1;
+		}
+	}
+	
+	for(t=pM->Pstart[p2]; t<pM->Pstart[p2+1]; ++t) {
+		if(rowsP1[pM->i[t]] == 1 && colsP1[pM->j[t]] == 1) {
+			/* This nonzero is free. As p2 is relatively large, move it */
+			
+			if(dummies && pM->i[t] == pM->j[t] && pM->dummy[pM->i[t]]) {
+				continue;
+			}
+			
+			nzWeight = (symmetric && pM->i[t] != pM->j[t])?2:1;
+			
+			if((nnz1+nzWeight)/(double)abs(procs[p1]) > (nnz2-nzWeight)/(double)abs(procs[p2])) {
+				free(rowsP1);
+				free(colsP1);
+				return TRUE;
+			}
+			
+			/* Move the nonzero */
+			if(p2 > p1) {
+				SwapNonzero(pM, t, pM->Pstart[p2]);
+				++pM->Pstart[p2];
+			}
+			else {
+				SwapNonzero(pM, t, pM->Pstart[p2+1]-1);
+				--pM->Pstart[p2+1];
+			}
+			
+			/* Update bookkeeping */
+			nnz1 += nzWeight;
+			nnz2 -= nzWeight;
+			
+		}
+		
+	}
+	
+	/* Finish */
+	free(rowsP1);
+	free(colsP1);
+	
+	return TRUE;
+} /* end ImproveFreeNonzeros */
diff --git a/src/FreeNonzeros.h b/src/FreeNonzeros.h
new file mode 100644
index 0000000000000000000000000000000000000000..13939e44b936059a3cd67acc52e5e8258dc801f7
--- /dev/null
+++ b/src/FreeNonzeros.h
@@ -0,0 +1,17 @@
+#ifndef __FreeNonzeros_h__
+#define __FreeNonzeros_h__
+
+#include "Sort.h"
+#include "SparseMatrix.h"
+#include "DistributeMat.h"
+
+struct freeIndex {
+	long t;
+	long numProcs;
+};
+
+void SwapNonzero(struct sparsematrix *pM, long s, long t);
+
+int ImproveFreeNonzeros(struct sparsematrix *pM, const struct opts *pOptions, const int *procs, int p1, int p2);
+
+#endif /* __FreeNonzeros_h__ */
diff --git a/src/Makefile b/src/Makefile
index dcf6aadb69436fe4dfbc8e89c394eea9a7bce66d..f47acb9d86fe82624e7d19c85e810144c641f5ce 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -35,7 +35,8 @@ SUBSETSUM=SubsetSum ${HEAP}
 ZEROVOLUMESEARCH=ZeroVolumeSearch ${SUBSETSUM}
 DISTRIBUTE=${DISTRIBUTEMAT} ${DISTRIBUTEVEC} ${DISTRIBUTEVECORIGEQ} ${ZEROVOLUMESEARCH}
 CARTESIAN=Cartesian ${DISTRIBUTEVECLIB} ${SPARSEMATRIX}
-MONDRIAANLIBRARY=${DISTRIBUTE} ${CARTESIAN}
+FREENONZEROS=FreeNonzeros
+MONDRIAANLIBRARY=${DISTRIBUTE} ${CARTESIAN} ${FREENONZEROS}
 
 #targets
 lib/libMondriaan${MONDRIAANMAJORVERSION}.a: ${MONDRIAANLIBRARY:%=%.c} ${MONDRIAANLIBRARY:%=%.o} Mondriaan.h
diff --git a/src/Options.c b/src/Options.c
index ca5e46590052e542c4fa25b4e6b345f56f07d661..ad15c763761cabead52c8f79c2fea53ad633a2b0 100644
--- a/src/Options.c
+++ b/src/Options.c
@@ -131,6 +131,7 @@ char* GetDefaultOptionText() {
 "Metric                                         lambda1 \n"
 "Discard_Free_Nets                              yes \n"
 "Zero_Volume_Search                             no \n"
+"Improve_Free_Nonzeros                          yes \n"
 "SquareMatrix_DistributeVectorsEqual            no \n"
 "SquareMatrix_DistributeVectorsEqual_AddDummies yes \n"
 "SymmetricMatrix_UseSingleEntry                 no \n"
@@ -252,6 +253,12 @@ int ExportOptions(FILE *Out, const struct opts *Opts) {
     else return FALSE;
     fprintf(Out, "\n");
     
+    fprintf(Out, "Improve_Free_Nonzeros ");
+    if (Opts->ImproveFreeNonzeros == FreeNonzerosYes) fprintf(Out, "yes");
+    else if (Opts->ImproveFreeNonzeros == FreeNonzerosNo) fprintf(Out, "no");
+    else return FALSE;
+    fprintf(Out, "\n");
+    
     fprintf(Out, "SquareMatrix_DistributeVectorsEqual ");
     if (Opts->SquareMatrix_DistributeVectorsEqual == EqVecNo) fprintf(Out, "no");
     else if (Opts->SquareMatrix_DistributeVectorsEqual == EqVecYes) fprintf(Out, "yes");
@@ -476,6 +483,12 @@ int ExportOptionsToLaTeX(FILE *Out, const struct opts *Opts) {
     else fprintf(Out, "?");
     fprintf(Out, " \\\\\n");
     
+    fprintf(Out, "Improve-Free-Nonzeros & ");
+    if (Opts->ImproveFreeNonzeros == FreeNonzerosYes) fprintf(Out, "yes");
+    else if (Opts->ImproveFreeNonzeros == FreeNonzerosNo) fprintf(Out, "no");
+    else fprintf(Out, "?");
+    fprintf(Out, " \\\\\n");
+    
     fprintf(Out, "DistributeVectorsEqual & ");
     if (Opts->SquareMatrix_DistributeVectorsEqual == EqVecNo) fprintf(Out, "no");
     else if (Opts->SquareMatrix_DistributeVectorsEqual == EqVecYes) fprintf(Out, "yes");
@@ -870,6 +883,15 @@ int SetOption(struct opts *pOptions, const char *option, const char *value) {
             fprintf(stderr, "SetOptions(): unknown %s '%s'!\n", option, value);
             return FALSE;
         }
+    } else if (!strcmp(option, "Improve_Free_Nonzeros")) {
+        if (!strcmp(value, "no") || !strcmp(value, "0"))
+            pOptions->ImproveFreeNonzeros = FreeNonzerosNo;
+        else if (!strcmp(value, "yes") || ! strcmp(value, "1"))
+            pOptions->ImproveFreeNonzeros = FreeNonzerosYes;
+        else {
+            fprintf(stderr, "SetOption(): unknown %s '%s'!\n", option, value);
+            return FALSE;
+        }
     } else if (!strcmp(option, "SquareMatrix_DistributeVectorsEqual")) {
         if (!strcmp(value, "no") || !strcmp(value, "0"))
             pOptions->SquareMatrix_DistributeVectorsEqual = EqVecNo;
diff --git a/src/Options.h b/src/Options.h
index 9a7a191470bef6a8df1d8ff1ffca1f80c97b5be1..087f2160e5b8bdd5d0cfb8696147e1e50da6f07f 100644
--- a/src/Options.h
+++ b/src/Options.h
@@ -58,6 +58,7 @@ struct opts {
     enum {MetricLambda, MetricCut, MetricLambdaLambdaMinusOne} Metric;
     enum {FreeNetYes, FreeNetNo} DiscardFreeNets;
     enum {ZeroVolNo, ZeroVolYes} ZeroVolumeSearch;
+    enum {FreeNonzerosNo, FreeNonzerosYes} ImproveFreeNonzeros;
   
     /* Matrix options */
     enum {EqVecNo, EqVecYes} SquareMatrix_DistributeVectorsEqual;
diff --git a/src/SparseMatrix.c b/src/SparseMatrix.c
index e07a3bed6b8f857f11b77ed9b14640bb7723d276..930372211e320c4cfda7671f32f8563ac1e5496e 100644
--- a/src/SparseMatrix.c
+++ b/src/SparseMatrix.c
@@ -2080,6 +2080,135 @@ int SparseMatrixSymmetric2Full(struct sparsematrix *pM) {
    
     return TRUE;
 } /* end SparseMatrixSymmetric2Full */
+
+
+int SparseMatrixFull2Symmetric(struct sparsematrix *pM, char MMTypeCode) {
+
+    /* This function converts a full sparse matrix to the
+       symmetric, skew-symmetric, or hermitian representation.
+
+       This function is the 'inverse' of SparseMatrixSymmetric2Full().
+       
+       See SparseMatrixSymmetric2Full() for more information.
+    */
+  
+    long t, tt, Ndiag;
+    long start, NrNzEltsNew;
+    int q;
+    
+    if (!pM) {
+        fprintf(stderr, "SparseMatrixFull2Symmetric(): Null parameter!\n");
+        return FALSE;
+    }
+  
+    /* Check if matrix is general square */
+    if (pM->m != pM->n || pM->MMTypeCode[3] != 'G') {
+        fprintf(stderr, "SparseMatrixFull2Symmetric(): matrix is not general square!\n");
+        return FALSE;
+    }
+    if (MMTypeCode != 'S' && MMTypeCode != 'K' && MMTypeCode != 'H') {
+        fprintf(stderr, "SparseMatrixFull2Symmetric(): type code must equal S, K or H!\n");
+        return FALSE;
+    }
+  
+    /* Count the number of diagonal elements */
+    Ndiag = 0;
+    for (t = 0; t < pM->NrNzElts; t++)
+        if (pM->i[t] == pM->j[t])
+            Ndiag++;
+    NrNzEltsNew = (pM->NrNzElts + Ndiag)/2;
+
+    if (MMTypeCode == 'K' && Ndiag > 0) {
+        fprintf(stderr, "SparseMatrixFull2Symmetric(): matrix is not skew-symmetric!\n");
+        return FALSE;
+    }
+
+    /* Update Pstart */
+    if (pM->NrProcs >= 1 && pM->Pstart != NULL) {
+        tt = 0;
+        start = pM->Pstart[0];
+        for (q = 0; q < pM->NrProcs; q++) {
+            for (t = start; t < pM->Pstart[q+1]; t++) {
+                if (pM->i[t] >= pM->j[t])
+                    tt++;
+            }
+            
+            start = pM->Pstart[q+1];
+            pM->Pstart[q+1] = tt;
+        }
+    
+        if(tt != NrNzEltsNew) {
+            fprintf(stderr, "SparseMatrixFull2Symmetric(): Processor weights not rearranged correctly!\n");
+            return FALSE;
+        }
+    }
+    
+    /* Copy the entries */
+    tt = 0;
+    for (t = 0; t < pM->NrNzElts; t++) {
+        /* Copy entry t into tt */
+        if (pM->i[t] < pM->j[t]) {
+            continue;
+        }
+
+        pM->i[tt] = pM->i[t];
+        pM->j[tt] = pM->j[t];
+        if (pM->MMTypeCode[2] != 'P')
+            pM->ReValue[tt] = pM->ReValue[t];
+        if (pM->MMTypeCode[2] == 'C')
+            pM->ImValue[tt] = pM->ImValue[t];
+        tt++;
+    }
+    
+    if(tt != NrNzEltsNew) {
+        fprintf(stderr, "SparseMatrixFull2Symmetric(): Nonzeros not rearranged correctly!\n");
+        return FALSE;
+    }
+   
+    /* Allocate memory for the entries */
+    pM->i = (long *) realloc(pM->i, (NrNzEltsNew+1) * sizeof(long));
+    pM->j = (long *) realloc(pM->j, (NrNzEltsNew+1) * sizeof(long));
+    
+    if (pM->i == NULL || pM->j == NULL) {
+        fprintf(stderr, "SparseMatrixFull2Symmetric(): Unable to reallocate!\n");
+        return FALSE;
+    }
+    
+    if (pM->MMTypeCode[2] != 'P') {
+        pM->ReValue = (double *) realloc(pM->ReValue, (NrNzEltsNew+1)*sizeof(double));
+        
+        if (pM->ReValue == NULL) {
+            fprintf(stderr, "SparseMatrixFull2Symmetric(): Unable to reallocate!\n");
+            return FALSE;
+        }
+    }
+    
+    if (pM->MMTypeCode[2] == 'C') {
+        pM->ImValue = (double *) realloc(pM->ImValue, (NrNzEltsNew+1)*sizeof(double));
+        
+        if (pM->ImValue == NULL) {
+            fprintf(stderr, "SparseMatrixFull2Symmetric(): Unable to reallocate!\n");
+            return FALSE;
+        }
+    }
+  
+    /* Update the number of nonzero entries and the matrix type code */
+    pM->NrNzElts = NrNzEltsNew;
+    pM->MMTypeCode[3] = MMTypeCode;
+    switch(MMTypeCode) {
+        case 'S':
+            strcpy(pM->Symmetry, "symmetric");
+            break;
+        case 'K':
+            strcpy(pM->Symmetry, "skew-symmetric");
+            break;
+        case 'H':
+            strcpy(pM->Symmetry, "hermitian");
+            break;
+    }
+   
+    return TRUE;
+} /* end SparseMatrixFull2Symmetric */
   
 
 int SparseMatrixSymmetricLower2Random(struct sparsematrix *pM) {
diff --git a/src/SparseMatrix.h b/src/SparseMatrix.h
index 932d9910ee722de366cbb158a8ed974a4b846449..372530212c50ff6d4736c78ae45a7eaa750f2632 100644
--- a/src/SparseMatrix.h
+++ b/src/SparseMatrix.h
@@ -130,6 +130,7 @@ int AddDummiesToSparseMatrix(struct sparsematrix *pM);
 int RemoveDummiesFromSparseMatrix(struct sparsematrix *pM);
 
 int SparseMatrixSymmetric2Full(struct sparsematrix *pM);
+int SparseMatrixFull2Symmetric(struct sparsematrix *pM, char MMTypeCode);
 int SparseMatrixSymmetricLower2Random(struct sparsematrix *pM);
 int SparseMatrixSymmetricRandom2Lower(struct sparsematrix *pM);
 
diff --git a/tests/Makefile b/tests/Makefile
index b9ad9441557114f7ed308b996813db0f7e291315..46ff30e7cde13fb68eb31cefd3fca1fac6894283 100644
--- a/tests/Makefile
+++ b/tests/Makefile
@@ -24,10 +24,10 @@ test_ComputeWeight test_SplitMatrixSimple test_SplitMatrixKLFM \
 test_DistributeMatrixMondriaan \
 test_CreateNewBiPartHyperGraph test_DeleteBiPartHyperGraph \
 test_SparseMatrix2BiPartHyperGraph test_BiPartHyperGraph2SparseMatrix \
-test_SparseMatrixSymmetric2Full test_SparseMatrixSymmetricLower2Random \
-test_SparseMatrixSymmetricRandom2Lower test_AddDummiesToSparseMatrix \
-test_RemoveDummiesFromSparseMatrix test_MMSparseMatrixInit \
-test_MMDeleteSparseMatrix test_MMSparseMatrixAllocateMemory \
+test_SparseMatrixSymmetric2Full test_SparseMatrixFull2Symmetric \
+test_SparseMatrixSymmetricLower2Random test_SparseMatrixSymmetricRandom2Lower \
+test_AddDummiesToSparseMatrix test_RemoveDummiesFromSparseMatrix \
+test_MMSparseMatrixInit test_MMDeleteSparseMatrix test_MMSparseMatrixAllocateMemory \
 test_MMSparseMatrixFreeMemory test_MMSparseMatrixReadEntries \
 test_MMSparseMatrixReadPstart test_MMSparseMatrixReadWeights \
 test_MMSparseMatrixReadTail test_MMSparseMatrixReadHeader \
@@ -49,7 +49,8 @@ test_FindOptimalPathMatching \
 test_MatchATA \
 test_CoarsenGraph test_RunMLGraphPart \
 test_SparseMatrixToCRS_CCS test_DetectConnectedComponents \
-test_Heap test_TwoColorTree test_SubsetSum test_ZeroVolumeSearch
+test_Heap test_TwoColorTree test_SubsetSum test_ZeroVolumeSearch \
+test_FreeNonzeros
 
 .SECONDARY: ${TESTTARGETS:%=%.o}
 
@@ -272,6 +273,9 @@ test_BiPartHyperGraph2SparseMatrix: test_BiPartHyperGraph2SparseMatrix.o ${GRAPH
 test_SparseMatrixSymmetric2Full: test_SparseMatrixSymmetric2Full.o ${SPARSEMATRIX:%=%.o} ${MATALLOC:%=%.o}
 	make -r OBJDEPS='$^' $@.target
 
+test_SparseMatrixFull2Symmetric: test_SparseMatrixFull2Symmetric.o ${SPARSEMATRIX:%=%.o}
+	make -r OBJDEPS='$^' $@.target
+
 test_SparseMatrixSymmetricLower2Random: test_SparseMatrixSymmetricLower2Random.o ${SPARSEMATRIX:%=%.o} ${MATALLOC:%=%.o}
 	make -r OBJDEPS='$^' $@.target
 
@@ -377,6 +381,9 @@ test_SubsetSum: test_SubsetSum.o ${SUBSETSUM:%=%.o}
 test_ZeroVolumeSearch: testHelper_DisconnectedMatrix.o test_ZeroVolumeSearch.o ${SPARSEMATRIX:%=%.o} ${DISTRIBUTEVECLIB:%=%.o} ${ZEROVOLUMESEARCH:%=%.o}
 	make -r OBJDEPS='$^' $@.target
 
+test_FreeNonzeros: test_FreeNonzeros.o ${SPARSEMATRIX:%=%.o} ${FREENONZEROS:%=%.o}
+	make -r OBJDEPS='$^' $@.target
+
 test_GetParameters: test_GetParameters.o ${OPTIONS:%=%.o}
 	make -r OBJDEPS='$^' $@.target
 
diff --git a/tests/Mondriaan.defaults b/tests/Mondriaan.defaults
index a936f7b92977da17817d548cc5265a70b3470575..8ebb27ba78f8827df414a4a5bd534b57a9001f31 100644
--- a/tests/Mondriaan.defaults
+++ b/tests/Mondriaan.defaults
@@ -8,6 +8,7 @@ Partitioner                                    mondriaan
 Metric                                         lambda1
 Discard_Free_Nets                              yes
 Zero_Volume_Search                             no
+Improve_Free_Nonzeros                          no
 SquareMatrix_DistributeVectorsEqual            no 
 SquareMatrix_DistributeVectorsEqual_AddDummies yes 
 SymmetricMatrix_UseSingleEntry                 no
diff --git a/tests/runtest b/tests/runtest
index 8c666000d1d8dcf076bf1b8a3ba789ea49a2136e..543aa8bca4569cbb8c65cabda59b0fc63fcd98a1 100755
--- a/tests/runtest
+++ b/tests/runtest
@@ -137,6 +137,7 @@ test_Heap
 test_TwoColorTree
 test_SubsetSum
 test_ZeroVolumeSearch
+test_FreeNonzeros
 
 # Test of Graph
 test_CreateNewBiPartHyperGraph
@@ -146,6 +147,7 @@ test_BiPartHyperGraph2SparseMatrix
 
 # Test of SparseMatrix 
 test_SparseMatrixSymmetric2Full
+test_SparseMatrixFull2Symmetric
 test_SparseMatrixSymmetricLower2Random
 test_SparseMatrixSymmetricRandom2Lower
 test_AddDummiesToSparseMatrix
diff --git a/tests/test_FreeNonzeros.c b/tests/test_FreeNonzeros.c
new file mode 100644
index 0000000000000000000000000000000000000000..98babc0399edf4caf0bd4ea531758e5849910ee5
--- /dev/null
+++ b/tests/test_FreeNonzeros.c
@@ -0,0 +1,184 @@
+#include "FreeNonzeros.h"
+#include "DistributeMat.h"
+#include "DistributeVecLib.h"
+#include <math.h>
+
+void test_FreeNonzeros(int symmetric);
+
+int main(int argc, char **argv) {
+
+    printf("Test FreeNonzeros: ");
+    
+    SetRandomSeed(111);
+    /*struct timeval tv;
+    gettimeofday(&tv,NULL);
+    unsigned long time_in_micros = 1000000 * tv.tv_sec + tv.tv_usec;
+    SetRandomSeed(time_in_micros);*/
+    
+    test_FreeNonzeros(FALSE);
+    test_FreeNonzeros(TRUE);
+    
+    printf("OK\n");
+    exit(0);
+
+} /* end main */
+
+/**
+ * Compute the total communication volume of a matrix.
+ */
+long ComputeVolume(struct sparsematrix *pM, int symmetric) {
+    if(symmetric) {
+        SparseMatrixSymmetric2Full(pM);
+    }
+    
+    long tmp, ComVolV, ComVolU;
+    CalcCom(pM, NULL, ROW, &ComVolV, &tmp, &tmp, &tmp, &tmp);
+    CalcCom(pM, NULL, COL, &ComVolU, &tmp, &tmp, &tmp, &tmp);
+    
+    if(symmetric) {
+        SparseMatrixFull2Symmetric(pM, 'S');
+    }
+    
+    return ComVolV+ComVolU;
+
+} /* end ComputeVolume */
+
+
+/**
+ * Test ImproveFreeNonzeros*()
+ * 
+ * Input:
+ * symmetric: Whether the matrix should be symmetric
+ */
+void test_FreeNonzeros(int symmetric) {
+    
+    struct sparsematrix A;
+    struct sparsematrix *pA = &A;
+    struct opts options;
+    
+    options.SymmetricMatrix_UseSingleEntry = symmetric ? SingleEntYes : SingleEntNo;
+    
+    long t, p, i, j;
+    long P = 5, m = 20, n = 20;
+    
+    /* Set up matrix struct */
+    MMSparseMatrixInit(pA);
+    pA->m = m;
+    pA->n = n;
+    pA->NrNzElts = P*(m+n)+((m-P)*(n-P))/5;
+    pA->NrProcs = P; /* maximum number of parts */
+
+    pA->i = (long *)malloc(pA->NrNzElts*sizeof(long));
+    pA->j = (long *)malloc(pA->NrNzElts*sizeof(long));
+    pA->Pstart = (long *)malloc((P+1)*sizeof(long));
+    
+    if (pA->i == NULL || pA->j == NULL || pA->Pstart == NULL) {
+        printf("Error\n");
+        exit(1);
+    }
+    
+    pA->MMTypeCode[0]='D'; /* normal matrix */
+    pA->MMTypeCode[1]='C'; /* coordinate scheme */
+    pA->MMTypeCode[2]='P'; /* pattern only */
+    if(symmetric)
+        pA->MMTypeCode[3]='S'; /* symmetric */
+    else
+        pA->MMTypeCode[3]='G'; /* general, no symmetry */
+    
+    /* Fill matrix */
+    t = 0;
+    for(p=0; p<P; ++p) {
+        /* First, fill rows and columns such that we have many processors per row/column */
+        pA->Pstart[p] = t;
+        for(i=P; i<pA->m; ++i) {
+            if(Random1(0, 2) == 0) {
+                continue;
+            }
+            pA->i[t] = i;
+            pA->j[t] = p;
+            ++t;
+        }
+        if(!symmetric) {
+            for(j=P; j<pA->n; ++j) {
+                if(Random1(0, 2) == 0) {
+                    continue;
+                }
+                pA->i[t] = p;
+                pA->j[t] = j;
+                ++t;
+            }
+        }
+    }
+    
+    /* Add additional nonzeros to last partition */
+    i = j = 0;
+    while(TRUE) {
+        j += Random1(5, 10);
+        if(symmetric) {
+            while(j > i) {
+                j -= i;
+                ++i;
+            }
+        }
+        else {
+            while(j >= pA->n-P) {
+                j -= pA->n-P;
+                ++i;
+            }
+        }
+        if(i >= pA->m-P) {
+            break;
+        }
+        pA->i[t] = P+i;
+        pA->j[t] = P+j;
+        ++t;
+    }
+    pA->Pstart[P] = pA->NrNzElts = t;
+    
+    /* We should provide a procs array */
+    int *procs = (int *)malloc(P*sizeof(int));
+    if (procs == NULL) {
+        printf("Error\n");
+        exit(1);
+    }
+    for(p=0; p<P; ++p) {
+        procs[p] = 1;
+    }
+    
+    /* Compute statistics before applying algorithm */
+    long totalImbalanceBefore = 0, weight;
+    long avgNrNzElts = pA->NrNzElts / P;
+    for(p=0; p<P; ++p) {
+        weight = ComputeWeight(pA, pA->Pstart[p], pA->Pstart[p+1]-1, NULL, &options);
+        totalImbalanceBefore += abs(weight-avgNrNzElts);
+    }
+    long volumeBefore = ComputeVolume(pA, symmetric);
+    
+    /* Run algorithm */
+    ImproveFreeNonzeros(pA, &options, procs, 3, 4);
+    
+    /* Compute statistics after applying algorithm */
+    long totalImbalanceAfter = 0.0;
+    for(p=0; p<P; ++p) {
+        weight = ComputeWeight(pA, pA->Pstart[p], pA->Pstart[p+1]-1, NULL, &options);
+        totalImbalanceAfter += abs(weight-avgNrNzElts);
+    }
+    long volumeAfter = ComputeVolume(pA, symmetric);
+    
+    /* Check that the total imbalance has not increased */
+    if(totalImbalanceAfter > totalImbalanceBefore) {
+        printf("Error1\n");
+        exit(1);
+    }
+    
+    /* Check that the volume has not increased */
+    if(volumeAfter > volumeBefore) {
+        printf("Error2\n");
+        exit(1);
+    }
+    
+    free(procs);
+    MMDeleteSparseMatrix(pA);
+    
+} /* end test_FreeNonzeros */
+
diff --git a/tests/test_SparseMatrixFull2Symmetric.c b/tests/test_SparseMatrixFull2Symmetric.c
new file mode 100644
index 0000000000000000000000000000000000000000..747c3a1cb7ed28a726afceaa31ed8f29f04069b6
--- /dev/null
+++ b/tests/test_SparseMatrixFull2Symmetric.c
@@ -0,0 +1,111 @@
+#include "SparseMatrix.h"
+#include <math.h>
+
+void test_SparseMatrixFull2Symmetric();
+
+int main(int argc, char **argv) {
+
+    printf("Test SparseMatrixFull2Symmetric: ");
+    
+    SetRandomSeed(111);
+    /*struct timeval tv;
+    gettimeofday(&tv,NULL);
+    unsigned long time_in_micros = 1000000 * tv.tv_sec + tv.tv_usec;
+    SetRandomSeed(time_in_micros);*/
+    
+    test_SparseMatrixFull2Symmetric();
+    
+    printf("OK\n");
+    exit(0);
+
+} /* end main */
+
+
+/**
+ * Test SparseMatrixFull2Symmetric()
+ */
+void test_SparseMatrixFull2Symmetric() {
+    
+    struct sparsematrix A, B;
+    struct sparsematrix *pA = &A, *pB = &B;
+    
+    long t, p, i, j;
+    long P = 5, n = 20;
+    
+    /* Set up matrix struct */
+    MMSparseMatrixInit(pA);
+    MMSparseMatrixInit(pB);
+    pA->m = pB->m = n;
+    pA->n = pB->n = n;
+    pA->NrNzElts = pB->NrNzElts = (n*n+n)/2;
+    pA->NrProcs = pB->NrProcs = P; /* maximum number of parts */
+
+    pA->i = (long *)malloc(pA->NrNzElts*sizeof(long));
+    pA->j = (long *)malloc(pA->NrNzElts*sizeof(long));
+    pA->Pstart = (long *)malloc((P+1)*sizeof(long));
+    
+    pB->i = (long *)malloc(pA->NrNzElts*sizeof(long));
+    pB->j = (long *)malloc(pA->NrNzElts*sizeof(long));
+    pB->Pstart = (long *)malloc((P+1)*sizeof(long));
+    
+    if (pA->i == NULL || pA->j == NULL || pA->Pstart == NULL || pB->i == NULL || pB->j == NULL || pB->Pstart == NULL) {
+        printf("Error\n");
+        exit(1);
+    }
+    
+    pA->MMTypeCode[0]=pB->MMTypeCode[0]='D'; /* normal matrix */
+    pA->MMTypeCode[1]=pB->MMTypeCode[1]='C'; /* coordinate scheme */
+    pA->MMTypeCode[2]=pB->MMTypeCode[2]='P'; /* pattern only */
+    pA->MMTypeCode[3]=pB->MMTypeCode[3]='S'; /* symmetric */
+    
+    /* Fill matrix */
+    t = 0;
+    i = j = 0;
+    while(TRUE) {
+        j += Random1(1, 5);
+        while(j > i) {
+            j -= i;
+            ++i;
+        }
+        if(i >= pA->m) {
+            break;
+        }
+        pA->i[t] = pB->i[t] = i;
+        pA->j[t] = pB->j[t] = j;
+        ++t;
+    }
+    pA->NrNzElts = pB->NrNzElts = t;
+    for(p=0; p<P; ++p) {
+        pA->Pstart[p] = pB->Pstart[p] = (pA->NrNzElts/P)*p;
+    }
+    pA->Pstart[P] = pB->Pstart[P] = pA->NrNzElts;
+    
+    /* Run procedures */
+    SparseMatrixSymmetric2Full(pA);
+    SparseMatrixFull2Symmetric(pA, 'S');
+    
+    /* Check result */
+    if(pA->NrNzElts != pB->NrNzElts || pA->NrProcs != P) {
+        printf("Error\n");
+        exit(1);
+    }
+    
+    for(p=0; p<=P; ++p) {
+        if(pA->Pstart[p] != pB->Pstart[p]) {
+            printf("Error\n");
+            exit(1);
+        }
+    }
+    
+    for(t=0; t<pA->NrNzElts; ++t) {
+        if(pA->i[t] != pB->i[t] || pA->j[t] != pB->j[t]) {
+            printf("Error\n");
+            exit(1);
+        }
+    }
+    
+    MMDeleteSparseMatrix(pA);
+    MMDeleteSparseMatrix(pB);
+    
+} /* end test_SparseMatrixFull2Symmetric */
+
diff --git a/tools/Mondriaan.defaults b/tools/Mondriaan.defaults
index 1176264f23e369444a444e2988b83d9336f5a285..96001cbd815aae98cab6428707fbc365943dc364 100644
--- a/tools/Mondriaan.defaults
+++ b/tools/Mondriaan.defaults
@@ -7,6 +7,7 @@ SplitMethod                                    KLFM
 Partitioner                                    mondriaan 
 Metric                                         lambda1 
 Discard_Free_Nets                              yes 
+Improve_Free_Nonzeros                          no
 SquareMatrix_DistributeVectorsEqual            no 
 SquareMatrix_DistributeVectorsEqual_AddDummies yes 
 SymmetricMatrix_UseSingleEntry                 no