diff --git a/AdaptiveTreeVisualization.vcxproj b/AdaptiveTreeVisualization.vcxproj
index 6e1765e..4300eb7 100644
--- a/AdaptiveTreeVisualization.vcxproj
+++ b/AdaptiveTreeVisualization.vcxproj
@@ -142,12 +142,15 @@
       <OpenMPSupport>true</OpenMPSupport>
       <AdditionalOptions>
       </AdditionalOptions>
+      <AdditionalIncludeDirectories>.</AdditionalIncludeDirectories>
     </ClCompile>
     <Link>
       <SubSystem>Console</SubSystem>
       <EnableCOMDATFolding>true</EnableCOMDATFolding>
       <OptimizeReferences>true</OptimizeReferences>
       <GenerateDebugInformation>true</GenerateDebugInformation>
+      <AdditionalLibraryDirectories>$(OutDir)</AdditionalLibraryDirectories>
+      <AdditionalDependencies>ZLIB.lib;JPEG.lib;PNG.lib;kernel32.lib;user32.lib;gdi32.lib;winspool.lib;comdlg32.lib;advapi32.lib;shell32.lib;ole32.lib;oleaut32.lib;uuid.lib;odbc32.lib;odbccp32.lib;%(AdditionalDependencies)</AdditionalDependencies>
     </Link>
   </ItemDefinitionGroup>
   <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
diff --git a/Makefile b/Makefile
index 6961475..f5f6c1b 100644
--- a/Makefile
+++ b/Makefile
@@ -126,10 +126,12 @@ make_dir:
 	$(MD) -p $(BIN)
 
 $(BIN)$(PR_TARGET): $(PR_OBJECTS)
-	$(CXX) -o $@ $(PR_OBJECTS) $(LFLAGS)
+	cd PNG  && make
+	$(CXX) -o $@ $(PR_OBJECTS) -L$(BIN) $(LFLAGS) -ljpeg -lmypng -lz
 
 $(BIN)$(SR_TARGET): $(SR_OBJECTS)
-	$(CXX) -o $@ $(SR_OBJECTS) $(LFLAGS)
+	cd PNG  && make
+	$(CXX) -o $@ $(SR_OBJECTS) -L$(BIN) $(LFLAGS) -ljpeg -lmypng -lz
 
 $(BIN)$(ST_TARGET): $(ST_OBJECTS)
 	$(CXX) -o $@ $(ST_OBJECTS) $(LFLAGS)
@@ -142,7 +144,8 @@ $(BIN)$(IS_TARGET): $(IS_OBJECTS)
 	$(CXX) -o $@ $(IS_OBJECTS) -L$(BIN) $(LFLAGS) -ljpeg -lmypng -lz
 
 $(BIN)$(AV_TARGET): $(AV_OBJECTS)
-	$(CXX) -o $@ $(AV_OBJECTS) $(LFLAGS)
+	cd PNG  && make
+	$(CXX) -o $@ $(AV_OBJECTS) -L$(BIN) $(LFLAGS) -ljpeg -lmypng -lz
 
 $(BIN)%.o: $(SRC)%.c
 	$(CC) -c -o $@ -I$(INCLUDE) $<
diff --git a/PoissonRecon.vcxproj b/PoissonRecon.vcxproj
index fa78d9c..b5acc2c 100644
--- a/PoissonRecon.vcxproj
+++ b/PoissonRecon.vcxproj
@@ -154,7 +154,7 @@
       <TargetEnvironment>X64</TargetEnvironment>
     </Midl>
     <ClCompile>
-      <AdditionalIncludeDirectories>%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
+      <AdditionalIncludeDirectories>.;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
       <PreprocessorDefinitions>WIN32;NDEBUG;_CONSOLE;_CRT_SECURE_NO_DEPRECATE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
       <RuntimeLibrary>MultiThreadedDLL</RuntimeLibrary>
       <PrecompiledHeader>
@@ -177,10 +177,11 @@
       <DataExecutionPrevention>
       </DataExecutionPrevention>
       <TargetMachine>MachineX64</TargetMachine>
-      <AdditionalDependencies>psapi.lib;kernel32.lib;user32.lib;gdi32.lib;winspool.lib;comdlg32.lib;advapi32.lib;shell32.lib;ole32.lib;oleaut32.lib;uuid.lib;odbc32.lib;odbccp32.lib;%(AdditionalDependencies)</AdditionalDependencies>
+      <AdditionalDependencies>ZLIB.lib;JPEG.lib;PNG.lib;psapi.lib;kernel32.lib;user32.lib;gdi32.lib;winspool.lib;comdlg32.lib;advapi32.lib;shell32.lib;ole32.lib;oleaut32.lib;uuid.lib;odbc32.lib;odbccp32.lib;%(AdditionalDependencies)</AdditionalDependencies>
       <IgnoreSpecificDefaultLibraries>
       </IgnoreSpecificDefaultLibraries>
       <OutputFile>$(OutDir)$(TargetName)$(TargetExt)</OutputFile>
+      <AdditionalLibraryDirectories>$(OutDir)</AdditionalLibraryDirectories>
     </Link>
   </ItemDefinitionGroup>
   <ItemGroup>
diff --git a/SSDRecon.vcxproj b/SSDRecon.vcxproj
index 3647f45..1a5bd3e 100644
--- a/SSDRecon.vcxproj
+++ b/SSDRecon.vcxproj
@@ -144,14 +144,16 @@
       <AdditionalOptions>
       </AdditionalOptions>
       <WholeProgramOptimization>true</WholeProgramOptimization>
+      <AdditionalIncludeDirectories>.</AdditionalIncludeDirectories>
     </ClCompile>
     <Link>
       <SubSystem>Console</SubSystem>
       <EnableCOMDATFolding>true</EnableCOMDATFolding>
       <OptimizeReferences>true</OptimizeReferences>
       <GenerateDebugInformation>true</GenerateDebugInformation>
-      <AdditionalDependencies>psapi.lib;kernel32.lib;user32.lib;gdi32.lib;winspool.lib;comdlg32.lib;advapi32.lib;shell32.lib;ole32.lib;oleaut32.lib;uuid.lib;odbc32.lib;odbccp32.lib;%(AdditionalDependencies)</AdditionalDependencies>
+      <AdditionalDependencies>ZLIB.lib;JPEG.lib;PNG.lib;psapi.lib;kernel32.lib;user32.lib;gdi32.lib;winspool.lib;comdlg32.lib;advapi32.lib;shell32.lib;ole32.lib;oleaut32.lib;uuid.lib;odbc32.lib;odbccp32.lib;%(AdditionalDependencies)</AdditionalDependencies>
       <LinkTimeCodeGeneration>UseLinkTimeCodeGeneration</LinkTimeCodeGeneration>
+      <AdditionalLibraryDirectories>$(OutDir)</AdditionalLibraryDirectories>
     </Link>
   </ItemDefinitionGroup>
   <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
diff --git a/Src/AdaptiveTreeVisualization.cpp b/Src/AdaptiveTreeVisualization.cpp
index 98d4ed8..ddf70d0 100644
--- a/Src/AdaptiveTreeVisualization.cpp
+++ b/Src/AdaptiveTreeVisualization.cpp
@@ -38,6 +38,7 @@ DAMAGE.
 #include "FEMTree.h"
 #include "Ply.h"
 #include "PointStreamData.h"
+#include "Image.h"
 
 cmdLineParameter< char* >
 	In( "in" ) ,
@@ -89,6 +90,65 @@ void ShowUsage( char* ex )
 	printf( "\t[--%s]\n" , Verbose.name );
 }
 
+template< typename Real , unsigned int Dim >
+bool WriteImage( const Real *values , int res , const char *fileName , bool verbose )
+{
+	if( Dim!=2 ) return false;
+	int resolution = 1;
+	for( int d=0 ; d<Dim ; d++ ) resolution *= res;
+
+	Real avg = 0;
+#pragma omp parallel for reduction( + : avg )
+	for( int i=0 ; i<resolution ; i++ ) avg += values[i];
+	avg /= (Real)resolution;
+
+	Real std = 0;
+#pragma omp parallel for reduction( + : std )
+	for( int i=0 ; i<resolution ; i++ ) std += ( values[i] - avg ) * ( values[i] - avg );
+	std = (Real)sqrt( std / resolution );
+
+	if( verbose ) printf( "Grid to image: [%.2f,%.2f] -> [0,255]\n" , avg - 2*std , avg + 2*std );
+
+	unsigned char *pixels = new unsigned char[ resolution*3 ];
+#pragma omp parallel for
+	for( int i=0 ; i<resolution ; i++ )
+	{
+		Real v = (Real)std::min< Real >( (Real)1. , std::max< Real >( (Real)-1. , ( values[i] - avg ) / (2*std ) ) );
+		v = (Real)( ( v + 1. ) / 2. * 256. );
+		unsigned char color = (unsigned char )std::min< Real >( (Real)255. , std::max< Real >( (Real)0. , v ) );
+		for( int c=0 ; c<3 ; c++ ) pixels[i*3+c ] = color;
+	}
+	bool success = ImageWriter::Write( fileName , pixels , res , res , 3 );
+	delete[] pixels;
+	return success;
+}
+
+template< typename Real , unsigned int Dim >
+void WriteGrid( const Real *values , int res , const char *fileName )
+{
+	int resolution = 1;
+	for( int d=0 ; d<Dim ; d++ ) resolution *= res;
+
+	char *ext = GetFileExtension( fileName );
+	FILE *fp = fopen( fileName , "wb" );
+	if( !fp ) ERROR_OUT( "Failed to open voxel file for writing: %s" , fileName );
+	else
+	{
+		fwrite( &res , sizeof(int) , 1 , fp );
+		if( typeid(Real)==typeid(float) ) fwrite( values , sizeof(float) , resolution , fp );
+		else
+		{
+			float *fValues = new float[resolution];
+			for( int i=0 ; i<resolution ; i++ ) fValues[i] = float( values[i] );
+			fwrite( fValues , sizeof(float) , resolution , fp );
+			delete[] fValues;
+		}
+		fclose( fp );
+		DeletePointer( values );
+	}
+	delete[] ext;
+}
+
 template< unsigned int Dim , class Real , unsigned int FEMSig >
 void _Execute( const FEMTree< Dim , Real >* tree , FILE* fp )
 {
@@ -100,39 +160,12 @@ void _Execute( const FEMTree< Dim , Real >* tree , FILE* fp )
 	// Output the grid
 	if( OutGrid.set )
 	{
-		FILE* _fp = fopen( OutGrid.value , "wb" );
-		if( !_fp ) WARN( "Failed to open grid file for writing: %s" , OutGrid.value );
-		else
-		{
-			int res = 0;
-			double t = Time();
-			Pointer( Real ) values = tree->template regularGridEvaluate< true >( coefficients , res , -1 , PrimalGrid.set );
-			if( Verbose.set ) printf( "Got grid: %.2f(s)\n" , Time()-t );
-
-			int cells = 1;
-			for( int d=0 ; d<Dim ; d++ ) cells *= res;
-			Pointer( float ) fValues = AllocPointer< float >( cells );
-
-			Real min , max;
-			min = max = values[0];
-			for( int i=0 ; i<cells ; i++ ) min = std::min< Real >( min , values[i] ) , max = std::max< Real >( max , values[i] );
-
-			int strides[Dim] , _strides[Dim] ; strides[0] = _strides[0] = 1;
-			int idx[Dim+1] , _idx[Dim+1] ; idx[0] = _idx[0] = 0;
-			for( int d=1 ; d<Dim ; d++ ) strides[d] = strides[d-1] * res , _strides[d] = _strides[d-1] * res;
-			WindowLoop< Dim >::Run
-			(
-				0 , res ,
-				[&]( int d , int i ){ _idx[d+1] = _idx[d] + _strides[d] * i , idx[d+1] = idx[d] + strides[d] * i; } ,
-				[&]( void ){ fValues[ _idx[Dim] ] = (float)values[ idx[Dim] ]; }
-			);
-			DeletePointer( values );
-
-			fwrite( &res , sizeof(int) , 1 , _fp );
-			fwrite( fValues , sizeof(float) , cells , _fp );
-			fclose( _fp );
-			FreePointer( fValues );
-		}
+		int res = 0;
+		double t = Time();
+		Pointer( Real ) values = tree->template regularGridEvaluate< true >( coefficients , res , -1 , PrimalGrid.set );
+		if( Verbose.set ) printf( "Got grid: %.2f(s)\n" , Time()-t );
+		if( !WriteImage< Real , Dim >( values , res , OutGrid.value , Verbose.set ) ) WriteGrid< Real , Dim >( values , res , OutGrid.value );
+		DeletePointer( values );
 	}
 
 	// Output the mesh
diff --git a/Src/EDTInHeat.cpp b/Src/EDTInHeat.cpp
index 79f0c0b..40a966a 100644
--- a/Src/EDTInHeat.cpp
+++ b/Src/EDTInHeat.cpp
@@ -567,12 +567,7 @@ int main( int argc , char* argv[] )
 	static const int Degree = DEFAULT_FEM_DEGREE;
 	static const BoundaryType BType = BOUNDARY_FREE;
 
-	WARN( "Compiled for degree-%d, boundary-%s, %s-precision _only_" , Degree , BoundaryNames[ BType ] , sizeof(DefaultFloatType)==4 ? "single" : "double" );
-#if 1
 	WARN( "Compiled for degree-%d, boundary-%s, %s-precision _only_" , Degree , BoundaryNames[ BType ] , sizeof(Real)==4 ? "single" : "double" );
-#else
-	WARN( "Compiled for degree-%d, boundary-%s, %s-precision _only_" , Degree , BoundaryNames[ BType ] , sizeof(DefaultFloatType)==4 ? "single" : "double" );
-#endif
 	if( BaseDepth.value>FullDepth.value )
 	{
 		if( BaseDepth.set ) WARN( "Base depth must be smaller than full depth: %d <= %d" , BaseDepth.value , FullDepth.value );
diff --git a/Src/FEMTree.IsoSurface.specialized.inl b/Src/FEMTree.IsoSurface.specialized.inl
index 9f6744e..4bc2f36 100644
--- a/Src/FEMTree.IsoSurface.specialized.inl
+++ b/Src/FEMTree.IsoSurface.specialized.inl
@@ -26,6 +26,8 @@ ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF S
 DAMAGE.
 */
 
+#include <sstream>
+#include <iomanip>
 #include <unordered_map>
 #include "MyMiscellany.h"
 #include "MarchingCubes.h"
@@ -1646,6 +1648,14 @@ public:
 		double cornersTime , verticesTime , edgesTime , surfaceTime;
 		double copyFinerTime , setTableTime;
 		IsoStats( void ) : cornersTime(0) , verticesTime(0) , edgesTime(0) , surfaceTime(0) , copyFinerTime(0) , setTableTime(0) {;}
+		std::string toString( void ) const
+		{
+			std::stringstream stream;
+			stream << "Corners / Vertices / Edges / Surface / Set Table / Copy Finer: ";
+			stream << std::fixed << std::setprecision(1) << cornersTime << " / " << verticesTime << " / " << edgesTime << " / " << surfaceTime << " / " << setTableTime << " / " << copyFinerTime;
+			stream << " (s)";
+			return stream.str();
+		}
 	};
 	template< typename Data , unsigned int ... FEMSigs , unsigned int WeightDegree , unsigned int DataSig >
 	static IsoStats Extract( UIntPack< FEMSigs ... > , UIntPack< WeightDegree > , UIntPack< DataSig > , const FEMTree< Dim , Real >& tree , const DensityEstimator< WeightDegree >* densityWeights , const SparseNodeData< ProjectiveData< Data , Real > , IsotropicUIntPack< Dim , DataSig > >* data , const DenseNodeData< Real , UIntPack< FEMSigs ... > >& coefficients , Real isoValue , CoredMeshData< Vertex >& mesh , std::function< void ( Vertex& , Point< Real , Dim > , Real , Data ) > SetVertex , bool nonLinearFit , bool addBarycenter , bool polygonMesh , bool flipOrientation )
diff --git a/Src/FEMTree.System.inl b/Src/FEMTree.System.inl
index fdc368b..316596e 100644
--- a/Src/FEMTree.System.inl
+++ b/Src/FEMTree.System.inl
@@ -242,7 +242,7 @@ int FEMTree< Dim , Real >::_solveFullSystemGS( UIntPack< FEMSigs ... > , const t
 	double bNorm=0 , inRNorm=0 , outRNorm=0;
 	if( depth>=0 )
 	{
-		SparseMatrix< Real > M;
+		SparseMatrix< Real , int , WindowSize< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >::Size > M;
 		double t = Time();
 		Pointer( Real ) D = AllocPointer< Real >( _sNodesEnd( depth ) - _sNodesBegin( depth ) );
 		Pointer( T ) _constraints = AllocPointer< T >( _sNodesSize( depth ) );
@@ -263,11 +263,11 @@ int FEMTree< Dim , Real >::_solveFullSystemGS( UIntPack< FEMSigs ... > , const t
 		Pointer( T ) X = GetPointer( &solution[0] + _sNodesBegin( depth ) , _sNodesSize( depth ) );
 		if( computeNorms )
 #pragma omp parallel for reduction( + : bNorm , inRNorm )
-			for( int j=0 ; j<M.rowNum ; j++ )
+			for( int j=0 ; j<M.rows() ; j++ )
 			{
 				T temp = {};
 				ConstPointer( MatrixEntry< Real > ) start = M[j];
-				ConstPointer( MatrixEntry< Real > ) end = start + (unsigned long long)M.rowSizes[j];
+				ConstPointer( MatrixEntry< Real > ) end = start + (unsigned long long)M.rowSize(j);
 				ConstPointer( MatrixEntry< Real > ) e;
 				for( e=start ; e!=end ; e++ ) temp += X[ e->N ] * e->Value;
 				bNorm += Dot( B[j] , B[j] );
@@ -282,11 +282,11 @@ int FEMTree< Dim , Real >::_solveFullSystemGS( UIntPack< FEMSigs ... > , const t
 
 		if( computeNorms )
 #pragma omp parallel for reduction( + : outRNorm )
-			for( int j=0 ; j<M.rowNum ; j++ )
+			for( int j=0 ; j<M.rows() ; j++ )
 			{
 				T temp = {};
 				ConstPointer( MatrixEntry< Real > ) start = M[j];
-				ConstPointer( MatrixEntry< Real > ) end = start + (unsigned long long)M.rowSizes[j];
+				ConstPointer( MatrixEntry< Real > ) end = start + (unsigned long long)M.rowSize(j);
 				ConstPointer( MatrixEntry< Real > ) e;
 				for( e=start ; e!=end ; e++ ) temp += X[ e->N ] * e->Value;
 				outRNorm += Dot( temp-B[j] , temp-B[j] );
@@ -408,7 +408,7 @@ int FEMTree< Dim , Real >::_solveSlicedSystemGS( UIntPack< FEMSigs ... > , const
 		// The number of in-core blocks over which we either solve or compute residuals
 		int matrixBlocks = std::max< int >( 1 , std::min< int >( solveBlocks+2*residualOffset , blockEnd-blockBegin ) );
 		// The list of matrices for each in-memory block
-		Pointer( SparseMatrix< Real > ) _M = NewPointer< SparseMatrix< Real > >( matrixBlocks );
+		Pointer( SparseMatrix< Real , int , WindowSize< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >::Size > ) _M = NewPointer< SparseMatrix< Real , int , WindowSize< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >::Size > >( matrixBlocks );
 		Pointer( Pointer( Real ) ) _D = AllocPointer< Pointer( Real ) >( matrixBlocks );
 		std::vector< Pointer( T ) > _constraints( matrixBlocks );
 		for( int i=0 ; i<matrixBlocks ; i++ ) _D[i] = NullPointer( Real ) , _constraints[i] = NullPointer( T );
@@ -421,6 +421,14 @@ int FEMTree< Dim , Real >::_solveSlicedSystemGS( UIntPack< FEMSigs ... > , const
 		// If we are solving forward we start in a block S with S mod ColorModulus = ColorModulus-1
 		// and end in a block E with E mod ColorModulus = 0
 		while( MOD( solveWindow.begin(!forward) , ColorModulus )!=( forward ? ColorModulus-1 : 0 ) ) solveWindow -= dir , residualWindow -= dir;
+		int maxBlockSize = 0;
+		BlockWindow _residualWindow = residualWindow;
+		for( ; _residualWindow.end(!forward)*dir<FullWindow.end(forward)*dir ; _residualWindow += dir )
+		{
+			int b = _residualWindow.begin(!forward);
+			if( FullWindow.inBlock( b ) ) maxBlockSize = std::max< int >( maxBlockSize , _sNodesEnd( depth , BlockLast( b ) ) - _sNodesBegin( depth , BlockFirst( b ) ) );
+		}
+		for( int i=0 ; i<matrixBlocks ; i++ ) _constraints[i] = AllocPointer< T >( maxBlockSize ) , _D[i] = AllocPointer< Real >( maxBlockSize );
 		for( ; residualWindow.end(!forward)*dir<FullWindow.end(forward)*dir ; residualWindow += dir , solveWindow += dir )
 		{
 			double t;
@@ -435,10 +443,6 @@ int FEMTree< Dim , Real >::_solveSlicedSystemGS( UIntPack< FEMSigs ... > , const
 					int b = residualBlock , _b = MOD( b , matrixBlocks );
 
 					t = Time();
-					FreePointer( _D[_b] );
-					_D[_b] = AllocPointer< Real >( _sNodesEnd( depth , BlockLast( b ) ) - _sNodesBegin( depth , BlockFirst( b ) ) );
-					FreePointer( _constraints[_b] );
-					_constraints[_b] = AllocPointer< T >( _sNodesEnd( depth , BlockLast( b ) ) - _sNodesBegin( depth , BlockFirst( b ) ) );
 					_getSliceMatrixAndProlongationConstraints( UIntPack< FEMSigs ... >() , F , _M[_b] , _D[_b] , bsData , depth , _sNodesBegin( depth , BlockFirst( b ) ) , _sNodesEnd( depth , BlockLast( b ) ) , prolongedSolution , _constraints[_b] , ccStencil , pcStencils , interpolationInfo... );
 #pragma omp parallel for
 					for( int i=_sNodesBegin( depth , BlockFirst( b ) ) ; i<_sNodesEnd( depth , BlockLast( b ) ) ; i++ ) _constraints[_b][ i - _sNodesBegin( depth , BlockFirst( b ) ) ] = constraints[i] - _constraints[_b][ i - _sNodesBegin( depth , BlockFirst( b ) ) ];
@@ -452,11 +456,11 @@ int FEMTree< Dim , Real >::_solveSlicedSystemGS( UIntPack< FEMSigs ... > , const
 						ConstPointer( T ) B = _constraints[_b];
 						ConstPointer( T ) X = XBlocks( depth , b , solution );
 #pragma omp parallel for reduction( + : bNorm , inRNorm )
-						for( int j=0 ; j<_M[_b].rowNum ; j++ )
+						for( int j=0 ; j<_M[_b].rows() ; j++ )
 						{
 							T temp = {};
 							ConstPointer( MatrixEntry< Real > ) start = _M[_b][j];
-							ConstPointer( MatrixEntry< Real > ) end = start + (unsigned long long)_M[_b].rowSizes[j];
+							ConstPointer( MatrixEntry< Real > ) end = start + (unsigned long long)_M[_b].rowSize(j);
 							ConstPointer( MatrixEntry< Real > ) e;
 							for( e=start ; e!=end ; e++ ) temp += X[ e->N ] * e->Value;
 							bNorm += Dot( B[j] , B[j] );
@@ -493,11 +497,11 @@ int FEMTree< Dim , Real >::_solveSlicedSystemGS( UIntPack< FEMSigs ... > , const
 					ConstPointer( T ) B = _constraints[_b];
 					ConstPointer( T ) X = XBlocks( depth , b , solution );
 #pragma omp parallel for reduction( + : outRNorm )
-					for( int j=0 ; j<_M[_b].rowNum ; j++ )
+					for( int j=0 ; j<_M[_b].rows() ; j++ )
 					{
 						T temp = {};
 						ConstPointer( MatrixEntry< Real > ) start = _M[_b][j];
-						ConstPointer( MatrixEntry< Real > ) end = start + (unsigned long long)_M[_b].rowSizes[j];
+						ConstPointer( MatrixEntry< Real > ) end = start + (unsigned long long)_M[_b].rowSize(j);
 						ConstPointer( MatrixEntry< Real > ) e;
 						for( e=start ; e!=end ; e++ ) temp += X[ e->N ] * e->Value;
 						outRNorm += Dot( temp-B[j] , temp-B[j] );
@@ -525,7 +529,7 @@ int FEMTree< Dim , Real >::_solveSystemCG( UIntPack< FEMSigs ... > , const typen
 	int iter = 0;
 	Pointer( T ) X = GetPointer( &solution[0] + _sNodesBegin(depth) , _sNodesSize(depth) );
 	ConstPointer( T ) B = GetPointer( &constraints[0] + _sNodesBegin(depth) , _sNodesSize(depth) );
-	SparseMatrix< Real > M;
+	SparseMatrix< Real , int , WindowSize< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >::Size > M;
 
 	double& systemTime = stats.systemTime;
 	double&  solveTime = stats. solveTime;
@@ -544,10 +548,10 @@ int FEMTree< Dim , Real >::_solveSystemCG( UIntPack< FEMSigs ... > , const typen
 	systemTime = Time()-systemTime;
 	solveTime = Time();
 	// Solve the linear system
-	accuracy = Real( accuracy / 100000 ) * M.rowNum;
+	accuracy = Real( accuracy / 100000 ) * M.rows();
 	int dims[] = { ( _BSplineEnd< FEMSigs >( depth ) - _BSplineBegin< FEMSigs >( depth ) ) ... };
 	int nonZeroRows = 0;
-	for( int i=0 ; i<M.rowNum ; i++ ) if( M.rowSizes[i] ) nonZeroRows++;
+	for( int i=0 ; i<M.rows() ; i++ ) if( M.rowSize(i) ) nonZeroRows++;
 	int totalDim = 1;
 	for( int d=0 ; d<Dim ; d++ ) totalDim *= dims[d];
 	BoundaryType bTypes[] = { FEMSignature< FEMSigs >::BType ... };
@@ -558,11 +562,11 @@ int FEMTree< Dim , Real >::_solveSystemCG( UIntPack< FEMSigs ... > , const typen
 	if( computeNorms )
 	{
 #pragma omp parallel for reduction( + : bNorm , inRNorm )
-		for( int j=0 ; j<M.rowNum ; j++ )
+		for( int j=0 ; j<M.rows() ; j++ )
 		{
 			T temp = {};
 			ConstPointer( MatrixEntry< Real > ) start = M[j];
-			ConstPointer( MatrixEntry< Real > ) end = start + (unsigned long long)M.rowSizes[j];
+			ConstPointer( MatrixEntry< Real > ) end = start + (unsigned long long)M.rowSize(j);
 			ConstPointer( MatrixEntry< Real > ) e;
 			for( e=start ; e!=end ; e++ ) temp += X[ e->N ] * e->Value;
 			bNorm += Dot( B[j] , B[j] );
@@ -574,33 +578,33 @@ int FEMTree< Dim , Real >::_solveSystemCG( UIntPack< FEMSigs ... > , const typen
 	struct SPDFunctor
 	{
 	protected:
-		const SparseMatrix< Real >& _M;
+		const SparseMatrix< Real , int , WindowSize< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >::Size >& _M;
 		bool _addDCTerm;
 	public:
-		SPDFunctor( const SparseMatrix< Real >& M , bool addDCTerm ) : _M(M) , _addDCTerm(addDCTerm){ }
+		SPDFunctor( const SparseMatrix< Real , int , WindowSize< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >::Size >& M , bool addDCTerm ) : _M(M) , _addDCTerm(addDCTerm){ }
 		void operator()( ConstPointer( T ) in , Pointer( T ) out ) const
 		{
 			_M.multiply( in , out );
 			if( _addDCTerm )
 			{
 				T average = {};
-				for( int i=0 ; i<_M.rowNum ; i++ ) average += in[i];
-				average /= _M.rowNum;
-				for( int i=0 ; i<_M.rowNum ; i++ ) out[i] += average;
+				for( int i=0 ; i<_M.rows() ; i++ ) average += in[i];
+				average /= _M.rows();
+				for( int i=0 ; i<_M.rows() ; i++ ) out[i] += average;
 			}
 		}
 	};
-	if( iters ) iter = SolveCG< SPDFunctor , T , Real >( SPDFunctor( M , addDCTerm ) , (int)M.rowNum , ( ConstPointer( T ) )B , iters , X , Real( accuracy ) , Dot );
+	if( iters ) iter = SolveCG< SPDFunctor , T , Real >( SPDFunctor( M , addDCTerm ) , (int)M.rows() , ( ConstPointer( T ) )B , iters , X , Real( accuracy ) , Dot );
 
 	solveTime = Time()-solveTime;
 	if( computeNorms )
 	{
 #pragma omp parallel for reduction( + : outRNorm )
-		for( int j=0 ; j<M.rowNum ; j++ )
+		for( int j=0 ; j<M.rows() ; j++ )
 		{
 			T temp = {};
 			ConstPointer( MatrixEntry< Real > ) start = M[j];
-			ConstPointer( MatrixEntry< Real > ) end = start + (unsigned long long)M.rowSizes[j];
+			ConstPointer( MatrixEntry< Real > ) end = start + (unsigned long long)M.rowSize(j);
 			ConstPointer( MatrixEntry< Real > ) e;
 			for( e=start ; e!=end ; e++ ) temp += X[ e->N ] * e->Value;
 			outRNorm += Dot( temp-B[j] , temp-B[j] );
@@ -657,11 +661,11 @@ void FEMTree< Dim , Real >::_solveRegularMG( UIntPack< FEMSigs ... > ,  typename
 		const SparseMatrix< Real , int >& _M = M.back();
 		ConstPointer( T ) _X = X.back();
 #pragma omp parallel for reduction( + : bNorm , inRNorm )
-		for( int j=0 ; j<_M.rowNum ; j++ )
+		for( int j=0 ; j<_M.rows() ; j++ )
 		{
 			T temp = {};
 			ConstPointer( MatrixEntry< Real > ) start = _M[j];
-			ConstPointer( MatrixEntry< Real > ) end = start + (unsigned long long)_M.rowSizes[j];
+			ConstPointer( MatrixEntry< Real > ) end = start + (unsigned long long)_M.rowSize(j);
 			ConstPointer( MatrixEntry< Real > ) e;
 			for( e=start ; e!=end ; e++ ) temp += _X[ e->N ] * e->Value;
 			bNorm += Dot( _B[j] , _B[j] );
@@ -689,24 +693,24 @@ void FEMTree< Dim , Real >::_solveRegularMG( UIntPack< FEMSigs ... > ,  typename
 			struct SPDFunctor
 			{
 			protected:
-				const SparseMatrix< Real >& _M;
+				const SparseMatrix< Real , int >& _M;
 				bool _addDCTerm;
 			public:
-				SPDFunctor( const SparseMatrix< Real >& M , bool addDCTerm ) : _M(M) , _addDCTerm(addDCTerm){ }
+				SPDFunctor( const SparseMatrix< Real , int  >& M , bool addDCTerm ) : _M(M) , _addDCTerm(addDCTerm){ }
 				void operator()( ConstPointer( T ) in , Pointer( T ) out ) const
 				{
 					_M.multiply( in , out );
 					if( _addDCTerm )
 					{
 						T average = {};
-						for( int i=0 ; i<_M.rowNum ; i++ ) average += in[i];
-						average /= _M.rowNum;
-						for( int i=0 ; i<_M.rowNum ; i++ ) out[i] += average;
+						for( int i=0 ; i<_M.rows() ; i++ ) average += in[i];
+						average /= _M.rows();
+						for( int i=0 ; i<_M.rows() ; i++ ) out[i] += average;
 					}
 				}
 			};
 			int nonZeroRows = 0;
-			for( int i=0 ; i<M[d].rowNum ; i++ ) if( M[d].rowSizes[i] ) nonZeroRows++;
+			for( int i=0 ; i<M[d].rows() ; i++ ) if( M[d].rowSize(i) ) nonZeroRows++;
 			int totalDim = 1;
 			int dims[] = { ( _BSplineEnd< FEMSigs >( depth ) - _BSplineBegin< FEMSigs >( depth ) ) ... };
 			for( int dd=0 ; dd<Dim ; dd++ ) totalDim *= dims[dd];
@@ -731,11 +735,11 @@ void FEMTree< Dim , Real >::_solveRegularMG( UIntPack< FEMSigs ... > ,  typename
 		const SparseMatrix< Real , int >& _M = M.back();
 		ConstPointer( T ) _X = X.back();
 #pragma omp parallel for reduction( + : outRNorm )
-		for( int j=0 ; j<_M.rowNum ; j++ )
+		for( int j=0 ; j<_M.rows() ; j++ )
 		{
 			T temp = {};
 			ConstPointer( MatrixEntry< Real > ) start = _M[j];
-			ConstPointer( MatrixEntry< Real > ) end = start + (unsigned long long)_M.rowSizes[j];
+			ConstPointer( MatrixEntry< Real > ) end = start + (unsigned long long)_M.rowSize(j);
 			ConstPointer( MatrixEntry< Real > ) e;
 			for( e=start ; e!=end ; e++ ) temp += _X[ e->N ] * e->Value;
 			outRNorm += Dot( temp-_B[j] , temp-_B[j] );
@@ -885,6 +889,70 @@ void FEMTree< Dim , Real >::_addPointValues( UIntPack< FEMSigs ... > , StaticWin
 	);
 }
 
+template< unsigned int Dim , class Real >
+template< typename T , unsigned int ... PointDs , unsigned int ... FEMSigs >
+T FEMTree< Dim , Real >::_setMatrixRowAndGetConstraintFromProlongation( UIntPack< FEMSigs ... > , const BaseSystem< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& F , const typename FEMTreeNode::template ConstNeighbors< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >& pNeighbors , const typename FEMTreeNode::template ConstNeighbors< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >& neighbors , size_t idx , SparseMatrix< Real , int , WindowSize< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >::Size > &M , int offset , const PCStencils< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& pcStencils , const CCStencil< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& ccStencil , const PointEvaluator< UIntPack< FEMSigs ... > , UIntPack< FEMSignature< FEMSigs >::Degree ... > >& bsData , ConstPointer( T ) prolongedSolution , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const
+{
+	T constraint ={};
+	typedef UIntPack< ( -BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapStart ) ... > OverlapRadii;
+	typedef UIntPack<    BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize    ... > OverlapSizes;
+
+	int count = 0;
+	const FEMTreeNode* node = neighbors.neighbors.data[ WindowIndex< OverlapSizes , OverlapRadii >::Index ];
+	Pointer( MatrixEntry< Real > ) row = M[idx];
+
+	LocalDepth d ; LocalOffset off;
+	_localDepthAndOffset( node , d , off );
+	if( d>0 && prolongedSolution )
+	{
+		int cIdx = (int)( node - node->parent->children );
+		constraint = _getConstraintFromProlongedSolution( UIntPack< FEMSigs ... >() , F , neighbors , pNeighbors , node , prolongedSolution , pcStencils.data[cIdx] , bsData , interpolationInfo... );
+	}
+
+	bool isInterior = BaseFEMIntegrator::IsInteriorlyOverlapped( UIntPack< FEMSignature< FEMSigs >::Degree ... >() , UIntPack< FEMSignature< FEMSigs >::Degree ... >() , d , off );
+
+	StaticWindow< Real , UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > > pointValues;
+	memset( pointValues.data , 0 , sizeof(Real)*WindowSize< OverlapSizes >::Size );
+	_addPointValues( UIntPack< FEMSigs ... >() , pointValues , neighbors , bsData , interpolationInfo ... );
+	int nodeIndex = node->nodeData.nodeIndex;
+	if( isInterior ) // General case, so try to make fast
+	{
+		const FEMTreeNode* const * _nodes = neighbors.neighbors.data;
+		ConstPointer( double ) _stencil = ccStencil.data;
+		Real* _values = pointValues.data;
+		row[count++] = MatrixEntry< Real >( nodeIndex-offset , (Real)( _values[ WindowIndex< OverlapSizes , OverlapRadii >::Index ] + _stencil[ WindowIndex< OverlapSizes , OverlapRadii >::Index ] ) );
+		for( int i=0 ; i<WindowSize< OverlapSizes >::Size ; i++ ) if( _isValidFEM1Node( _nodes[i] ) )
+		{
+			if( i!=WindowIndex< OverlapSizes , OverlapRadii >::Index ) row[count++] = MatrixEntry< Real >( _nodes[i]->nodeData.nodeIndex-offset , (Real)( _values[i] + _stencil[i] ) );
+		}
+	}
+	else
+	{
+		LocalDepth d ; LocalOffset off;
+		_localDepthAndOffset( node , d , off );
+		Real temp = (Real)F.ccIntegrate( off , off ) + pointValues.data[ WindowIndex< OverlapSizes , OverlapRadii >::Index ];
+
+		row[count++] = MatrixEntry< Real >( nodeIndex-offset , temp );
+		LocalOffset _off;
+		WindowLoop< Dim >::Run
+		(
+			ZeroUIntPack< Dim >() , OverlapSizes() ,
+			[&]( int d , int i ){ _off[d] = off[d] - (int)OverlapRadii::Values[d] + i; } ,
+			[&]( const FEMTreeNode* _node , Real pointValue )
+		{
+			if( node!=_node && FEMIntegrator::IsValidFEMNode( UIntPack< FEMSigs ... >() , d , _off ) )
+			{
+				Real temp = (Real)F.ccIntegrate( _off , off ) + pointValue;
+				if( _isValidFEM1Node( _node ) ) row[count++] = MatrixEntry< Real >( _node->nodeData.nodeIndex-offset , temp );
+			}
+		} ,
+			neighbors.neighbors() , pointValues()
+			);
+	}
+	M.setRowSize( idx , count );
+	return constraint;
+}
+
 template< unsigned int Dim , class Real >
 template< typename T , unsigned int ... PointDs , unsigned int ... FEMSigs >
 T FEMTree< Dim , Real >::_setMatrixRowAndGetConstraintFromProlongation( UIntPack< FEMSigs ... > , const BaseSystem< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& F , const typename FEMTreeNode::template ConstNeighbors< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >& pNeighbors , const typename FEMTreeNode::template ConstNeighbors< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >& neighbors , Pointer( MatrixEntry< Real > ) row , int offset , const PCStencils< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& pcStencils , const CCStencil< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& ccStencil , const PointEvaluator< UIntPack< FEMSigs ... > , UIntPack< FEMSignature< FEMSigs >::Degree ... > >& bsData , ConstPointer( T ) prolongedSolution , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const
@@ -1501,7 +1569,7 @@ CumulativeDerivativeValues< T , Dim , PointD > FEMTree< Dim , Real >::_finerFunc
 
 template< unsigned int Dim , class Real >
 template< unsigned int ... FEMSigs , typename T , unsigned int ... PointDs >
-int FEMTree< Dim , Real >::_getSliceMatrixAndProlongationConstraints( UIntPack< FEMSigs ... > , const typename BaseFEMIntegrator::template System< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& F , SparseMatrix< Real >& matrix , Pointer( Real ) diagonalR , const PointEvaluator< UIntPack< FEMSigs ... > , UIntPack< FEMSignature< FEMSigs >::Degree ... > >& bsData , LocalDepth depth , int nBegin , int nEnd , ConstPointer( T ) prolongedSolution , Pointer( T ) constraints , const CCStencil< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& ccStencil , const PCStencils< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& pcStencils , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const
+int FEMTree< Dim , Real >::_getSliceMatrixAndProlongationConstraints( UIntPack< FEMSigs ... > , const typename BaseFEMIntegrator::template System< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& F , SparseMatrix< Real , int , WindowSize< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >::Size >& matrix , Pointer( Real ) diagonalR , const PointEvaluator< UIntPack< FEMSigs ... > , UIntPack< FEMSignature< FEMSigs >::Degree ... > >& bsData , LocalDepth depth , int nBegin , int nEnd , ConstPointer( T ) prolongedSolution , Pointer( T ) constraints , const CCStencil< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& ccStencil , const PCStencils< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& pcStencils , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const
 {
 	typedef UIntPack< FEMSignature< FEMSigs >::Degree ... > FEMDegrees;
 	typedef UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > OverlapSizes;
@@ -1518,18 +1586,9 @@ int FEMTree< Dim , Real >::_getSliceMatrixAndProlongationConstraints( UIntPack<
 		// Get the matrix row size	
 		typename FEMTreeNode::template ConstNeighbors< OverlapSizes > neighbors , pNeighbors;
 		neighborKey.getNeighbors( OverlapRadii() , OverlapRadii() , node , pNeighbors , neighbors );
-#if defined( __GNUC__ ) && __GNUC__ < 5
-#warning "you've got me gcc version<5"
-		int count = _getMatrixRowSize( UIntPack< FEMSigs ... >() , neighbors );
-#else // !__GNUC__ || __GNUC__ >=5
-		int count = _getMatrixRowSize< FEMSigs ... >( neighbors );
-#endif // __GNUC__ || __GNUC__ < 4
-		// Allocate memory for the row
-		matrix.setRowSize( i , count );
-
 		// Set the row entries
-		if( constraints ) constraints[i] = _setMatrixRowAndGetConstraintFromProlongation( UIntPack< FEMSigs ... >() , F , pNeighbors , neighbors , matrix[i] , nBegin , pcStencils , ccStencil , bsData , prolongedSolution , interpolationInfo... );
-		else                               _setMatrixRowAndGetConstraintFromProlongation( UIntPack< FEMSigs ... >() , F , pNeighbors , neighbors , matrix[i] , nBegin , pcStencils , ccStencil , bsData , prolongedSolution , interpolationInfo... );
+		if( constraints ) constraints[i] = _setMatrixRowAndGetConstraintFromProlongation( UIntPack< FEMSigs ... >() , F , pNeighbors , neighbors , i , matrix , nBegin , pcStencils , ccStencil , bsData , prolongedSolution , interpolationInfo... );
+		else                               _setMatrixRowAndGetConstraintFromProlongation( UIntPack< FEMSigs ... >() , F , pNeighbors , neighbors , i , matrix , nBegin , pcStencils , ccStencil , bsData , prolongedSolution , interpolationInfo... );
 		if( diagonalR ) diagonalR[i] = (Real)1. / matrix[i][0].Value;
 	}
 	else if( constraints ) constraints[i] = T();
@@ -1547,12 +1606,12 @@ int FEMTree< Dim , Real >::_getSliceMatrixAndProlongationConstraints( UIntPack<
 
 template< unsigned int Dim , class Real >
 template< typename T , unsigned int ... PointDs , unsigned int ... FEMSigs >
-SparseMatrix< Real > FEMTree< Dim , Real >::systemMatrix( UIntPack< FEMSigs ... > , typename BaseFEMIntegrator::template System< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& F , LocalDepth depth , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const
+SparseMatrix< Real , int > FEMTree< Dim , Real >::systemMatrix( UIntPack< FEMSigs ... > , typename BaseFEMIntegrator::template System< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& F , LocalDepth depth , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const
 {
 	_setFEM1ValidityFlags( UIntPack< FEMSigs ... >() );
 	typedef typename BaseFEMIntegrator::template System< UIntPack< FEMSignature< FEMSigs >::Degree ... > > BaseSystem;
 	if( depth<0 || depth>_maxDepth ) ERROR_OUT( "System depth out of bounds: %d <= %d <= %d" , 0 , depth , _maxDepth );
-	SparseMatrix< Real > matrix;
+	SparseMatrix< Real , int > matrix;
 	F.init( depth );
 	PointEvaluator< UIntPack< FEMSigs ... > , UIntPack< FEMSignature< FEMSigs >::Degree ... > > bsData( depth );
 
@@ -1588,13 +1647,13 @@ SparseMatrix< Real > FEMTree< Dim , Real >::systemMatrix( UIntPack< FEMSigs ...
 
 template< unsigned int Dim , class Real >
 template< typename T , unsigned int ... PointDs , unsigned int ... FEMSigs >
-SparseMatrix< Real > FEMTree< Dim , Real >::prolongedSystemMatrix( UIntPack< FEMSigs ... > , typename BaseFEMIntegrator::template System< UIntPack<FEMSignature< FEMSigs >::Degree ... > >& F , LocalDepth highDepth , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const
+SparseMatrix< Real , int > FEMTree< Dim , Real >::prolongedSystemMatrix( UIntPack< FEMSigs ... > , typename BaseFEMIntegrator::template System< UIntPack<FEMSignature< FEMSigs >::Degree ... > >& F , LocalDepth highDepth , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const
 {
 	_setFEM1ValidityFlags( UIntPack< FEMSigs ... >() );
 	if( highDepth<=0 || highDepth>_maxDepth ) ERROR_OUT( "System depth out of bounds: %d < %d <= %d" , 0 , highDepth , _maxDepth );
 
 	LocalDepth lowDepth = highDepth-1;
-	SparseMatrix< Real > matrix;
+	SparseMatrix< Real , int > matrix;
 	F.init( highDepth );
 	PointEvaluator< UIntPack< FEMSigs ... > , UIntPack< FEMSignature< FEMSigs >::Degree ... > > bsData( highDepth );
 	typedef UIntPack< ( -BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapStart ) ... > OverlapRadii;
@@ -1625,9 +1684,9 @@ SparseMatrix< Real > FEMTree< Dim , Real >::prolongedSystemMatrix( UIntPack< FEM
 
 template< unsigned int Dim , class Real >
 template< unsigned int ... FEMSigs >
-SparseMatrix< Real > FEMTree< Dim , Real >::downSampleMatrix( UIntPack< FEMSigs ... > , LocalDepth highDepth ) const
+SparseMatrix< Real , int > FEMTree< Dim , Real >::downSampleMatrix( UIntPack< FEMSigs ... > , LocalDepth highDepth ) const
 {
-	SparseMatrix< Real > matrix;
+	SparseMatrix< Real , int > matrix;
 	_setFEM1ValidityFlags( UIntPack< FEMSigs ... >() );
 	typedef UIntPack< FEMSignature< FEMSigs >::Degree ... > FEMDegrees;
 	typedef UIntPack< BSplineSupportSizes< FEMSignature< FEMSigs >::Degree >::UpSampleSize ... > UpSampleSizes;
@@ -1714,56 +1773,56 @@ SparseMatrix< Real > FEMTree< Dim , Real >::downSampleMatrix( UIntPack< FEMSigs
 
 template< unsigned int Dim , class Real >
 template< typename T , unsigned int ... PointDs , unsigned int ... FEMSigs >
-SparseMatrix< Real > FEMTree< Dim , Real >::fullSystemMatrix( UIntPack< FEMSigs ... > , typename BaseFEMIntegrator::template System< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& F , LocalDepth depth , bool nonRefinableOnly , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const
+SparseMatrix< Real , int > FEMTree< Dim , Real >::fullSystemMatrix( UIntPack< FEMSigs ... > , typename BaseFEMIntegrator::template System< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& F , LocalDepth depth , bool nonRefinableOnly , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const
 {
-	SparseMatrix< Real > M;
-	std::vector< SparseMatrix< Real > >                  systemMatrices( depth+1 );
-	std::vector< SparseMatrix< Real > >         prolongedSystemMatrices( depth   );
-	std::vector< std::vector< SparseMatrix< Real > > > upSampleMatrices( depth-1 );
+	SparseMatrix< Real , int > M;
+	std::vector< SparseMatrix< Real , int > >                  systemMatrices( depth+1 );
+	std::vector< SparseMatrix< Real , int > >         prolongedSystemMatrices( depth   );
+	std::vector< std::vector< SparseMatrix< Real , int > > > upSampleMatrices( depth-1 );
 
 	for( int d=0 ; d<depth-1 ; d++ ) upSampleMatrices[d].resize( depth );
 	unsigned int size = _sNodesEnd( depth );
 	for( int d=0 ; d<=depth ; d++ )
 	{
-		SparseMatrix< Real >& M = systemMatrices[d];
+		SparseMatrix< Real , int >& M = systemMatrices[d];
 		M.resize( size );
-		SparseMatrix< Real > _M = systemMatrix< Real >( UIntPack< FEMSigs ... >() , F , d , interpolationInfo ... );
+		SparseMatrix< Real , int > _M = systemMatrix< Real >( UIntPack< FEMSigs ... >() , F , d , interpolationInfo ... );
 #pragma omp parallel for
 		for( int i=0 ; i<_M.rows() ; i++ )
 		{
-			M.setRowSize( i + _sNodesBegin(d) , _M.rowSizes[i] );
-			for( int j=0 ; j<_M.rowSizes[i] ; j++ ) M[i+_sNodesBegin(d)][j] = MatrixEntry< Real >( _M[i][j].N + _sNodesBegin(d) , _M[i][j].Value );
+			M.setRowSize( i + _sNodesBegin(d) , _M.rowSize(i) );
+			for( int j=0 ; j<_M.rowSize(i) ; j++ ) M[i+_sNodesBegin(d)][j] = MatrixEntry< Real >( _M[i][j].N + _sNodesBegin(d) , _M[i][j].Value );
 		}
 	}
 	for( int d=0 ; d<depth ; d++ )
 	{
-		SparseMatrix< Real >& M = prolongedSystemMatrices[d];
+		SparseMatrix< Real , int >& M = prolongedSystemMatrices[d];
 		M.resize( size );
-		SparseMatrix< Real > _M = prolongedSystemMatrix< Real >( UIntPack< FEMSigs ... >() , F , d+1 , interpolationInfo ... );
+		SparseMatrix< Real , int > _M = prolongedSystemMatrix< Real >( UIntPack< FEMSigs ... >() , F , d+1 , interpolationInfo ... );
 #pragma omp parallel for
 		for( int i=0 ; i<_M.rows() ; i++ )
 		{
-			M.setRowSize( i + _sNodesBegin(d+1) , _M.rowSizes[i] );
-			for( int j=0 ; j<_M.rowSizes[i] ; j++ ) M[i+_sNodesBegin(d+1)][j] = MatrixEntry< Real >( _M[i][j].N + _sNodesBegin(d) , _M[i][j].Value );
+			M.setRowSize( i + _sNodesBegin(d+1) , _M.rowSize(i) );
+			for( int j=0 ; j<_M.rowSize(i) ; j++ ) M[i+_sNodesBegin(d+1)][j] = MatrixEntry< Real >( _M[i][j].N + _sNodesBegin(d) , _M[i][j].Value );
 		}
 	}
 	for( int d=0 ; d<depth-1 ; d++ )
 	{
-		SparseMatrix< Real >& M = upSampleMatrices[d][d+1];
+		SparseMatrix< Real , int >& M = upSampleMatrices[d][d+1];
 		M.resize( size );
-		SparseMatrix< Real > _M = downSampleMatrix( UIntPack< FEMSigs ... >() , d+1 ).transpose( _sNodesSize( d+1 ) );
+		SparseMatrix< Real , int > _M = downSampleMatrix( UIntPack< FEMSigs ... >() , d+1 ).transpose( _sNodesSize( d+1 ) );
 #pragma omp parallel for
 		for( int i=0 ; i<_M.rows() ; i++ )
 		{
-			M.setRowSize( i + _sNodesBegin(d+1) , _M.rowSizes[i] );
-			for( int j=0 ; j<_M.rowSizes[i] ; j++ ) M[i+_sNodesBegin(d+1)][j] = MatrixEntry< Real >( _M[i][j].N + _sNodesBegin(d) , _M[i][j].Value );
+			M.setRowSize( i + _sNodesBegin(d+1) , _M.rowSize(i) );
+			for( int j=0 ; j<_M.rowSize(i) ; j++ ) M[i+_sNodesBegin(d+1)][j] = MatrixEntry< Real >( _M[i][j].N + _sNodesBegin(d) , _M[i][j].Value );
 		}
 		for( int dd=0 ; dd<d ; dd++ ) upSampleMatrices[dd][d+1] = upSampleMatrices[d][d+1] * upSampleMatrices[dd][d];
 	}
 
 	auto Matrix = [&]( int d1 , int d2 )
 	{
-		SparseMatrix< Real > _M;
+		SparseMatrix< Real , int > _M;
 		int _d1 = d1<d2 ? d1 : d2 , _d2 = d2<d1 ? d1 : d2;
 		if     ( _d1==_d2   ) _M =          systemMatrices[_d1];
 		else if( _d2==_d1+1 ) _M = prolongedSystemMatrices[_d2-1];
@@ -1777,13 +1836,13 @@ SparseMatrix< Real > FEMTree< Dim , Real >::fullSystemMatrix( UIntPack< FEMSigs
 		M += Matrix( d1 , d1 );
 		for( int d2=0 ; d2<=depth ; d2++ ) if( d1!=d2 )
 		{
-			SparseMatrix< Real > _M = Matrix( d1 , d2 );
+			SparseMatrix< Real , int > _M = Matrix( d1 , d2 );
 #pragma omp parallel for
-			for( int i=0 ; i<_M.rows() ; i++ ) if( _M.rowSizes[i] )
+			for( int i=0 ; i<_M.rows() ; i++ ) if( _M.rowSize(i) )
 			{
-				size_t oldSize = M.rowSizes[i];
-				M.resetRowSize( i , oldSize + _M.rowSizes[i] );
-				for( int j=0 ; j<_M.rowSizes[i] ; j++ ) M[i][oldSize+j] = _M[i][j];
+				size_t oldSize = M.rowSize(i);
+				M.resetRowSize( i , oldSize + _M.rowSize(i) );
+				for( int j=0 ; j<_M.rowSize(i) ; j++ ) M[i][oldSize+j] = _M[i][j];
 			}
 		}
 	}
@@ -1802,8 +1861,8 @@ SparseMatrix< Real > FEMTree< Dim , Real >::fullSystemMatrix( UIntPack< FEMSigs
 			else
 			{
 				int jj=0;
-				for( int j=0 ; j<M.rowSizes[i] ; j++ ) if( !( _isRefinableNode( _sNodes.treeNodes[ M[i][j].N ] ) && _localDepth( _sNodes.treeNodes[ M[i][j].N ] )<depth ) && _isValidFEM1Node( _sNodes.treeNodes[i] ) ) M[i][jj++] = M[i][j];
-				if( jj!=M.rowSizes[i] ) M.resetRowSize( i , jj );
+				for( int j=0 ; j<M.rowSize(i) ; j++ ) if( !( _isRefinableNode( _sNodes.treeNodes[ M[i][j].N ] ) && _localDepth( _sNodes.treeNodes[ M[i][j].N ] )<depth ) && _isValidFEM1Node( _sNodes.treeNodes[i] ) ) M[i][jj++] = M[i][j];
+				if( jj!=M.rowSize(i) ) M.resetRowSize( i , jj );
 			}
 	}
 	return M;
diff --git a/Src/FEMTree.h b/Src/FEMTree.h
index ea5e8d6..767403c 100644
--- a/Src/FEMTree.h
+++ b/Src/FEMTree.h
@@ -42,7 +42,7 @@ DAMAGE.
 #ifndef FEM_TREE_INCLUDED
 #define FEM_TREE_INCLUDED
 
-#define VERSION "10.05"
+#define VERSION "10.06"
 #define MEMORY_ALLOCATOR_BLOCK_SIZE 1<<12
 
 #define NEW_CODE
@@ -55,6 +55,7 @@ DAMAGE.
 #include "RegularTree.h"
 #include "SparseMatrix.h"
 #include <functional>
+#include <string>
 
 
 template< unsigned int Dim , class Real > class FEMTree;
@@ -1940,6 +1941,8 @@ protected:
 	template< unsigned int ... FEMSigs >
 	int _getMatrixRowSize( const typename FEMTreeNode::template ConstNeighbors< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >& neighbors ) const;
 #endif // __GNUC__ || __GNUC__ < 4
+	template< typename T , unsigned int ... PointDs , unsigned int ... FEMSigs >
+	T _setMatrixRowAndGetConstraintFromProlongation( UIntPack< FEMSigs ... > , const BaseSystem< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& F , const typename FEMTreeNode::template ConstNeighbors< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >& pNeighbors , const typename FEMTreeNode::template ConstNeighbors< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >& neighbors , size_t idx , SparseMatrix< Real , int , WindowSize< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >::Size > &M , int offset , const PCStencils< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& pcStencils , const CCStencil< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& ccStencil , const PointEvaluator< UIntPack< FEMSigs ... > , UIntPack< FEMSignature< FEMSigs >::Degree ... > >& bsData , ConstPointer( T ) prolongedSolution , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const;
 	template< typename T , unsigned int ... PointDs , unsigned int ... FEMSigs >
 	T _setMatrixRowAndGetConstraintFromProlongation( UIntPack< FEMSigs ... > , const BaseSystem< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& F , const typename FEMTreeNode::template ConstNeighbors< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >& pNeighbors , const typename FEMTreeNode::template ConstNeighbors< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >& neighbors , Pointer( MatrixEntry< Real > ) row , int offset , const PCStencils< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& pcStencils , const CCStencil< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& ccStencil , const PointEvaluator< UIntPack< FEMSigs ... > , UIntPack< FEMSignature< FEMSigs >::Degree ... > >& bsData , ConstPointer( T ) prolongedSolution , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const;
 	template< typename T , unsigned int ... PointDs , unsigned int ... FEMSigs >
@@ -1974,7 +1977,7 @@ protected:
 	CumulativeDerivativeValues< T , Dim , PointD >   _finerFunctionValues( UIntPack< FEMSigs ... > , Point< Real , Dim > p , const ConstPointSupportKey< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& neighborKey , const FEMTreeNode* node , const PointEvaluator< UIntPack< FEMSigs ... > , UIntPack< FEMSignature< FEMSigs >::Degree ... > >& bsData , ConstPointer( T ) coefficients ) const;
 
 	template< unsigned int ... FEMSigs , typename T , unsigned int ... PointDs >
-	int _getSliceMatrixAndProlongationConstraints( UIntPack< FEMSigs ... > , const BaseSystem< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& F , SparseMatrix< Real >& matrix , Pointer( Real ) diagonalR , const PointEvaluator< UIntPack< FEMSigs ... > , UIntPack< FEMSignature< FEMSigs >::Degree ... > >& bsData , LocalDepth depth , int nBegin , int nEnd , ConstPointer( T ) prolongedSolution , Pointer( T ) constraints , const CCStencil < UIntPack< FEMSignature< FEMSigs >::Degree ... > >& ccStencil , const PCStencils< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& pcStencils , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const;
+	int _getSliceMatrixAndProlongationConstraints( UIntPack< FEMSigs ... > , const BaseSystem< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& F , SparseMatrix< Real , int , WindowSize< UIntPack< BSplineOverlapSizes< FEMSignature< FEMSigs >::Degree >::OverlapSize ... > >::Size >& matrix , Pointer( Real ) diagonalR , const PointEvaluator< UIntPack< FEMSigs ... > , UIntPack< FEMSignature< FEMSigs >::Degree ... > >& bsData , LocalDepth depth , int nBegin , int nEnd , ConstPointer( T ) prolongedSolution , Pointer( T ) constraints , const CCStencil < UIntPack< FEMSignature< FEMSigs >::Degree ... > >& ccStencil , const PCStencils< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& pcStencils , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const;
 
 	// Down samples constraints @(depth) to constraints @(depth-1)
 	template< class C , unsigned ... Degrees , unsigned int ... FEMSigs > void _downSample( UIntPack< FEMSigs ... > , typename BaseFEMIntegrator::template RestrictionProlongation< UIntPack< Degrees ... > >& RP , LocalDepth highDepth , Pointer( C ) constraints ) const;
@@ -2422,14 +2425,14 @@ public:
 	}
 
 	template< typename T , unsigned int ... PointDs , unsigned int ... FEMSigs >
-	SparseMatrix< Real > systemMatrix( UIntPack< FEMSigs ... > , typename BaseFEMIntegrator::System< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& F , LocalDepth depth , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const;
+	SparseMatrix< Real , int > systemMatrix( UIntPack< FEMSigs ... > , typename BaseFEMIntegrator::System< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& F , LocalDepth depth , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const;
 	template< typename T , unsigned int ... PointDs , unsigned int ... FEMSigs >
-	SparseMatrix< Real > prolongedSystemMatrix( UIntPack< FEMSigs ... > , typename BaseFEMIntegrator::System< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& F , LocalDepth highDepth , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const;
+	SparseMatrix< Real , int > prolongedSystemMatrix( UIntPack< FEMSigs ... > , typename BaseFEMIntegrator::System< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& F , LocalDepth highDepth , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const;
 	template< unsigned int ... FEMSigs >
-	SparseMatrix< Real > downSampleMatrix( UIntPack< FEMSigs ... > , LocalDepth highDepth ) const;
+	SparseMatrix< Real , int > downSampleMatrix( UIntPack< FEMSigs ... > , LocalDepth highDepth ) const;
 
 	template< typename T , unsigned int ... PointDs , unsigned int ... FEMSigs >
-	SparseMatrix< Real > fullSystemMatrix( UIntPack< FEMSigs ... > , typename BaseFEMIntegrator::System< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& F , LocalDepth depth , bool nonRefinableOnly , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const;
+	SparseMatrix< Real , int > fullSystemMatrix( UIntPack< FEMSigs ... > , typename BaseFEMIntegrator::System< UIntPack< FEMSignature< FEMSigs >::Degree ... > >& F , LocalDepth depth , bool nonRefinableOnly , const InterpolationInfo< T , PointDs >* ... interpolationInfo ) const;
 
 	struct SolverInfo
 	{
@@ -2543,7 +2546,10 @@ template< unsigned int Dim , class Real > double FEMTree< Dim , Real >::_LocalMe
 template< unsigned int Dim , class Real , class Vertex >
 struct IsoSurfaceExtractor
 {
-	struct IsoStats{};
+	struct IsoStats
+	{
+		std::string toString( void ) const { return std::string( "Iso-surface extraction not supported for dimension %d" , Dim ); }
+	};
 	template< typename Data , unsigned int ... FEMSigs , unsigned int WeightDegree , unsigned int DataSig >
 	static IsoStats Extract
 	(
diff --git a/Src/FEMTree.inl b/Src/FEMTree.inl
index bf00596..b07a188 100644
--- a/Src/FEMTree.inl
+++ b/Src/FEMTree.inl
@@ -482,7 +482,8 @@ void FEMTree< Dim , Real >::finalizeForMultigrid( LocalDepth fullDepth , const H
 	for( LocalDepth d=_maxDepth-1 ; d>=0 ; d-- )
 	{
 		std::vector< FEMTreeNode* > nodes;
-		for( FEMTreeNode* node=_tree->nextNode() ; node ; node=_tree->nextNode( node ) ) if( _localDepth( node )==d && IsActiveNode< Dim >( node->children ) ) nodes.push_back( node );
+		auto NodeTerminationLambda = [&]( const FEMTreeNode *node ){ return _localDepth( node )==d; };
+		for( FEMTreeNode* node=_tree->nextNode( NodeTerminationLambda , NULL ) ; node ; node=_tree->nextNode( NodeTerminationLambda , node ) ) if( _localDepth( node )==d && IsActiveNode< Dim >( node->children ) ) nodes.push_back( node );
 #pragma omp parallel for
 		for( int i=0 ; i<nodes.size() ; i++ )
 		{
@@ -505,6 +506,7 @@ void FEMTree< Dim , Real >::finalizeForMultigrid( LocalDepth fullDepth , const H
 template< unsigned int Dim , class Real >
 void FEMTree< Dim , Real >::_setSpaceValidityFlags( void ) const
 {
+#pragma omp parallel for
 	for( int i=0 ; i<_sNodes.size() ; i++ )
 	{
 		const unsigned char MASK = ~( FEMTreeNodeData::SPACE_FLAG );
@@ -606,12 +608,17 @@ template< unsigned int Dim , class Real >
 template< class HasDataFunctor >
 void FEMTree< Dim , Real >::_clipTree( const HasDataFunctor& f , LocalDepth fullDepth )
 {
-	for( FEMTreeNode* temp=_tree->nextNode() ; temp ; temp=_tree->nextNode(temp) ) if( temp->children && _localDepth( temp )>=fullDepth )
-	{
-		bool hasData = false;
-		for( int c=0 ; c<(1<<Dim) && !hasData ; c++ ) hasData |= f( temp->children + c );
-		for( int c=0 ; c<(1<<Dim) ; c++ ) SetGhostFlag< Dim >( temp->children+c , !hasData );
-	}
+	std::vector< FEMTreeNode * > nodes;
+	auto NodeTerminationLambda = [&]( const FEMTreeNode *node ){ return _localDepth( node )==fullDepth; };
+	for( FEMTreeNode* temp=_tree->nextNode( NodeTerminationLambda , NULL ) ; temp ; temp=_tree->nextNode( NodeTerminationLambda , temp ) ) if( _localDepth( temp )==fullDepth ) nodes.push_back( temp );
+#pragma omp parallel for
+	for( int i=0 ; i<nodes.size() ; i++ )
+		for( FEMTreeNode* node=nodes[i]->nextNode() ; node ; node=nodes[i]->nextNode(node) ) if( node->children )
+		{
+			bool hasData = false;
+			for( int c=0 ; c<(1<<Dim) && !hasData ; c++ ) hasData |= f( node->children + c );
+			for( int c=0 ; c<(1<<Dim) ; c++ ) SetGhostFlag< Dim >( node->children+c , !hasData );
+		}
 }
 
 template< unsigned int Dim , class Real >
diff --git a/Src/PNG.inl b/Src/PNG.inl
index da5daae..d024f78 100644
--- a/Src/PNG.inl
+++ b/Src/PNG.inl
@@ -99,7 +99,6 @@ inline bool PNGReader::GetInfo( const char* fileName , unsigned int& width , uns
 	return true;
 }
 
-#if 1
 PNGWriter::PNGWriter( const char* fileName , unsigned int width , unsigned int height , unsigned int channels , unsigned int quality )
 {
 	_currentRow = 0;
@@ -148,52 +147,3 @@ unsigned int PNGWriter::nextRows( const unsigned char* rows , unsigned int rowNu
 	for( unsigned int r=0 ; r<rowNum ; r++ ) png_write_row( _png_ptr , (png_bytep)( rows + r * 3 * sizeof( unsigned char ) * _png_ptr->width ) );
 	return _currentRow += rowNum;
 }
-#else
-
-void PNGWriteColor( const char* fileName , const unsigned char* pixels , int width , int height )
-{
-	FILE* fp = fopen( fileName , "wb" );
-	if( !fp ) ERROR_OUT( "Failed to open file for writing: %s" , fileName );
-	png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING,0,0,0);
-	if(!png_ptr)	return;
-	png_infop info_ptr = png_create_info_struct(png_ptr);
-	if(!info_ptr)	return;
-	png_init_io(png_ptr, fp);
-	// turn off compression or set another filter
-	// png_set_filter(png_ptr, 0, PNG_FILTER_NONE);
-	png_set_IHDR(png_ptr, info_ptr, width , height ,
-		8,PNG_COLOR_TYPE_RGB,
-		PNG_INTERLACE_NONE,
-		PNG_COMPRESSION_TYPE_DEFAULT,
-		PNG_FILTER_TYPE_DEFAULT);
-	if (0) {                    // high-level write
-		std::vector<unsigned char> matrix( width * height * 3 );
-		std::vector<png_bytep> row_pointers( height );
-		for(int y=0;y<height;y++)
-		{
-			row_pointers[y]=&matrix[y*width*3];
-			unsigned char* buf=&matrix[y*width*3];
-			for(int x=0;x<width;x++)
-				for(int z=0;z<3;z++)
-					*buf++ = pixels[ (y*width+x)*3 + z ];
-		}
-		png_set_rows(png_ptr, info_ptr, &row_pointers[0]);
-		int png_transforms=0;
-		png_write_png(png_ptr, info_ptr, png_transforms, NULL);
-	} else {                    // low-level write
-		png_write_info(png_ptr, info_ptr);
-		// png_set_filler(png_ptr, 0, PNG_FILLER_AFTER);
-		//  but no way to provide GRAY data with RGBA fill, so pack each row
-		std::vector<unsigned char> buffer(width*3);
-		for(int y=0;y<height;y++)
-		{
-			unsigned char* buf=&buffer[0];
-			for(int x=0;x<width;x++)
-				for(int z=0;z<3;z++)
-					*buf++ = pixels[ (y*width+x)*3+z ];
-			png_bytep row_pointer=&buffer[0];
-			png_write_row(png_ptr, row_pointer);
-		}
-	}
-}
-#endif
\ No newline at end of file
diff --git a/Src/PoissonRecon.cpp b/Src/PoissonRecon.cpp
index 2276188..f0d4146 100644
--- a/Src/PoissonRecon.cpp
+++ b/Src/PoissonRecon.cpp
@@ -48,6 +48,7 @@ DAMAGE.
 #include "FEMTree.h"
 #include "Ply.h"
 #include "PointStreamData.h"
+#include "Image.h"
 
 MessageWriter messageWriter;
 
@@ -339,7 +340,8 @@ void ExtractMesh( UIntPack< FEMSigs ... > , std::tuple< SampleData ... > , FEMTr
 	else isoStats = IsoSurfaceExtractor< Dim , Real , Vertex >::template Extract< TotalPointSampleData >( Sigs() , UIntPack< WEIGHT_DEGREE >() , UIntPack< DataSig >() , tree , density , NULL , solution , isoValue , *mesh , SetVertex , !LinearFit.set , !NonManifold.set , PolygonMesh.set , false );
 #endif // __GNUC__ || __GNUC__ < 4
 	messageWriter( "Vertices / Polygons: %d / %d\n" , mesh->outOfCorePointCount()+mesh->inCorePoints.size() , mesh->polygonCount() );
-	messageWriter( "Corners / Vertices / Edges / Surface / Set Table / Copy Finer: %.1f / %.1f / %.1f / %.1f / %.1f / %.1f (s)\n" , isoStats.cornersTime , isoStats.verticesTime , isoStats.edgesTime , isoStats.surfaceTime , isoStats.setTableTime , isoStats.copyFinerTime );
+	std::string isoStatsString = isoStats.toString() + std::string( "\n" );
+	messageWriter( isoStatsString.c_str() );
 	if( PolygonMesh.set ) profiler.dumpOutput2( comments , "#         Got polygons:" );
 	else                  profiler.dumpOutput2( comments , "#        Got triangles:" );
 
@@ -350,6 +352,65 @@ void ExtractMesh( UIntPack< FEMSigs ... > , std::tuple< SampleData ... > , FEMTr
 	delete mesh;
 }
 
+template< typename Real , unsigned int Dim >
+bool WriteImage( const Real *values , int res , const char *fileName , bool verbose )
+{
+	if( Dim!=2 ) return false;
+	int resolution = 1;
+	for( int d=0 ; d<Dim ; d++ ) resolution *= res;
+
+	Real avg = 0;
+#pragma omp parallel for reduction( + : avg )
+	for( int i=0 ; i<resolution ; i++ ) avg += values[i];
+	avg /= (Real)resolution;
+
+	Real std = 0;
+#pragma omp parallel for reduction( + : std )
+	for( int i=0 ; i<resolution ; i++ ) std += ( values[i] - avg ) * ( values[i] - avg );
+	std = (Real)sqrt( std / resolution );
+
+	if( verbose ) printf( "Grid to image: [%.2f,%.2f] -> [0,255]\n" , avg - 2*std , avg + 2*std );
+
+	unsigned char *pixels = new unsigned char[ resolution*3 ];
+#pragma omp parallel for
+	for( int i=0 ; i<resolution ; i++ )
+	{
+		Real v = (Real)std::min< Real >( (Real)1. , std::max< Real >( (Real)-1. , ( values[i] - avg ) / (2*std ) ) );
+		v = (Real)( ( v + 1. ) / 2. * 256. );
+		unsigned char color = (unsigned char )std::min< Real >( (Real)255. , std::max< Real >( (Real)0. , v ) );
+		for( int c=0 ; c<3 ; c++ ) pixels[i*3+c ] = color;
+	}
+	bool success = ImageWriter::Write( fileName , pixels , res , res , 3 );
+	delete[] pixels;
+	return success;
+}
+
+template< typename Real , unsigned int Dim >
+void WriteGrid( const Real *values , int res , const char *fileName )
+{
+	int resolution = 1;
+	for( int d=0 ; d<Dim ; d++ ) resolution *= res;
+
+	char *ext = GetFileExtension( fileName );
+	FILE *fp = fopen( fileName , "wb" );
+	if( !fp ) ERROR_OUT( "Failed to open voxel file for writing: %s" , fileName );
+	else
+	{
+		fwrite( &res , sizeof(int) , 1 , fp );
+		if( typeid(Real)==typeid(float) ) fwrite( values , sizeof(float) , resolution , fp );
+		else
+		{
+			float *fValues = new float[resolution];
+			for( int i=0 ; i<resolution ; i++ ) fValues[i] = float( values[i] );
+			fwrite( fValues , sizeof(float) , resolution , fp );
+			delete[] fValues;
+		}
+		fclose( fp );
+		DeletePointer( values );
+	}
+	delete[] ext;
+}
+
 template< class Real , typename ... SampleData , unsigned int ... FEMSigs >
 void Execute( int argc , char* argv[] , UIntPack< FEMSigs ... > )
 {
@@ -601,28 +662,16 @@ void Execute( int argc , char* argv[] , UIntPack< FEMSigs ... > )
 
 	if( VoxelGrid.set )
 	{
-		FILE* fp = fopen( VoxelGrid.value , "wb" );
-		if( !fp ) WARN( "Failed to open voxel file for writing: %s" , VoxelGrid.value );
-		else
-		{
-			int res = 0;
-			profiler.start();
-			Pointer( Real ) values = tree.template regularGridEvaluate< true >( solution , res , -1 , PrimalVoxel.set );
+		int res = 0;
+		profiler.start();
+		Pointer( Real ) values = tree.template regularGridEvaluate< true >( solution , res , -1 , PrimalVoxel.set );
+		int resolution = 1;
+		for( int d=0 ; d<Dim ; d++ ) resolution *= res;
 #pragma omp parallel for
-			for( int i=0 ; i<res*res*res ; i++ ) values[i] -= isoValue;
-			profiler.dumpOutput( "Got voxel grid:" );
-			fwrite( &res , sizeof(int) , 1 , fp );
-			if( typeid(Real)==typeid(float) ) fwrite( values , sizeof(float) , res*res*res , fp );
-			else
-			{
-				float *fValues = new float[res*res*res];
-				for( int i=0 ; i<res*res*res ; i++ ) fValues[i] = float( values[i] );
-				fwrite( fValues , sizeof(float) , res*res*res , fp );
-				delete[] fValues;
-			}
-			fclose( fp );
-			DeletePointer( values );
-		}
+		for( int i=0 ; i<resolution ; i++ ) values[i] -= isoValue;
+		profiler.dumpOutput( "Got voxel grid:" );
+		if( !WriteImage< Real , DIMENSION >( values , res , VoxelGrid.value , Verbose.set ) ) WriteGrid< Real , DIMENSION >( values , res , VoxelGrid.value );
+		DeletePointer( values );
 	}
 
 	if( Out.set )
diff --git a/Src/RegularTree.h b/Src/RegularTree.h
index f7c5892..5e1a1c2 100644
--- a/Src/RegularTree.h
+++ b/Src/RegularTree.h
@@ -78,6 +78,13 @@ public:
 
 	const RegularTreeNode* nextLeaf( const RegularTreeNode* currentLeaf=NULL ) const;
 	RegularTreeNode* nextLeaf( RegularTreeNode* currentLeaf=NULL );
+
+	// This lambda takes a RegularTreeNode* as an argument and returns true if we do not need to traverse the tree beyond the specified node.
+	template< typename NodeTerminationLambda >
+	const RegularTreeNode* nextNode( NodeTerminationLambda &ntl , const RegularTreeNode* currentNode ) const;
+	template< typename NodeTerminationLambda >
+	RegularTreeNode* nextNode( NodeTerminationLambda &ntl , RegularTreeNode* currentNode );
+
 	const RegularTreeNode* nextNode( const RegularTreeNode* currentNode=NULL ) const;
 	RegularTreeNode* nextNode( RegularTreeNode* currentNode=NULL );
 	const RegularTreeNode* nextBranch( const RegularTreeNode* current ) const;
diff --git a/Src/RegularTree.inl b/Src/RegularTree.inl
index 2382344..bafec9a 100644
--- a/Src/RegularTree.inl
+++ b/Src/RegularTree.inl
@@ -311,6 +311,23 @@ RegularTreeNode< Dim , NodeData , DepthAndOffsetType >* RegularTreeNode< Dim , N
 	else return temp->nextLeaf();
 }
 
+template< unsigned int Dim , class NodeData , class DepthAndOffsetType >
+template< typename NodeTerminationLambda >
+const RegularTreeNode< Dim , NodeData , DepthAndOffsetType >* RegularTreeNode< Dim , NodeData , DepthAndOffsetType >::nextNode( NodeTerminationLambda &ntl , const RegularTreeNode *current ) const
+{
+	if( !current ) return this;
+	else if( current->children && !ntl(current) ) return current->children;
+	else return nextBranch( current );
+}
+template< unsigned int Dim , class NodeData , class DepthAndOffsetType >
+template< typename NodeTerminationLambda >
+RegularTreeNode< Dim , NodeData , DepthAndOffsetType >* RegularTreeNode< Dim , NodeData , DepthAndOffsetType >::nextNode( NodeTerminationLambda &ntl , RegularTreeNode* current )
+{
+	if( !current ) return this;
+	else if( current->children && !ntl(current) ) return current->children;
+	else return nextBranch( current );
+}
+
 template< unsigned int Dim , class NodeData , class DepthAndOffsetType >
 const RegularTreeNode< Dim , NodeData , DepthAndOffsetType >* RegularTreeNode< Dim , NodeData , DepthAndOffsetType >::nextNode( const RegularTreeNode* current ) const
 {
diff --git a/Src/SSDRecon.cpp b/Src/SSDRecon.cpp
index 22817dc..9845c10 100644
--- a/Src/SSDRecon.cpp
+++ b/Src/SSDRecon.cpp
@@ -48,6 +48,7 @@ DAMAGE.
 #include "FEMTree.h"
 #include "Ply.h"
 #include "PointStreamData.h"
+#include "Image.h"
 
 MessageWriter messageWriter;
 
@@ -292,7 +293,11 @@ template< unsigned int Dim , typename Real , typename TotalPointSampleData >
 struct SystemDual
 {
 	CumulativeDerivativeValues< Real , Dim , 1 > weight;
-	SystemDual( Real v , Real g ) : weight( v , g , g , g ) { }
+	SystemDual( Real v , Real g )
+	{
+		weight[0] = v;
+		for( int d=0 ; d<Dim ; d++ ) weight[d+1] = g;
+	}
 	CumulativeDerivativeValues< Real , Dim , 1 > operator()( Point< Real , Dim > p , const TotalPointSampleData& data , const CumulativeDerivativeValues< Real , Dim , 1 >& dValues ) const
 	{
 		return dValues * weight;
@@ -360,7 +365,8 @@ void ExtractMesh( UIntPack< FEMSigs ... > , std::tuple< SampleData ... > , FEMTr
 	else isoStats = IsoSurfaceExtractor< Dim , Real , Vertex >::template Extract< TotalPointSampleData >( Sigs() , UIntPack< WEIGHT_DEGREE >() , UIntPack< DataSig >() , tree , density , NULL , solution , isoValue , *mesh , SetVertex , NonLinearFit.set , !NonManifold.set , PolygonMesh.set , false );
 #endif // __GNUC__ || __GNUC__ < 4
 	messageWriter( "Vertices / Polygons: %d / %d\n" , mesh->outOfCorePointCount()+mesh->inCorePoints.size() , mesh->polygonCount() );
-	messageWriter( "Corners / Vertices / Edges / Surface / Set Table / Copy Finer: %.1f / %.1f / %.1f / %.1f / %.1f / %.1f (s)\n" , isoStats.cornersTime , isoStats.verticesTime , isoStats.edgesTime , isoStats.surfaceTime , isoStats.setTableTime , isoStats.copyFinerTime );
+	std::string isoStatsString = isoStats.toString() + std::string( "\n" );
+	messageWriter( isoStatsString.c_str() );
 	if( PolygonMesh.set ) profiler.dumpOutput2( comments , "#         Got polygons:" );
 	else                  profiler.dumpOutput2( comments , "#        Got triangles:" );
 
@@ -370,6 +376,65 @@ void ExtractMesh( UIntPack< FEMSigs ... > , std::tuple< SampleData ... > , FEMTr
 	delete mesh;
 }
 
+template< typename Real , unsigned int Dim >
+bool WriteImage( const Real *values , int res , const char *fileName , bool verbose )
+{
+	if( Dim!=2 ) return false;
+	int resolution = 1;
+	for( int d=0 ; d<Dim ; d++ ) resolution *= res;
+
+	Real avg = 0;
+#pragma omp parallel for reduction( + : avg )
+	for( int i=0 ; i<resolution ; i++ ) avg += values[i];
+	avg /= (Real)resolution;
+
+	Real std = 0;
+#pragma omp parallel for reduction( + : std )
+	for( int i=0 ; i<resolution ; i++ ) std += ( values[i] - avg ) * ( values[i] - avg );
+	std = (Real)sqrt( std / resolution );
+
+	if( verbose ) printf( "Grid to image: [%.2f,%.2f] -> [0,255]\n" , avg - 2*std , avg + 2*std );
+
+	unsigned char *pixels = new unsigned char[ resolution*3 ];
+#pragma omp parallel for
+	for( int i=0 ; i<resolution ; i++ )
+	{
+		Real v = (Real)std::min< Real >( (Real)1. , std::max< Real >( (Real)-1. , ( values[i] - avg ) / (2*std ) ) );
+		v = (Real)( ( v + 1. ) / 2. * 256. );
+		unsigned char color = (unsigned char )std::min< Real >( (Real)255. , std::max< Real >( (Real)0. , v ) );
+		for( int c=0 ; c<3 ; c++ ) pixels[i*3+c ] = color;
+	}
+	bool success = ImageWriter::Write( fileName , pixels , res , res , 3 );
+	delete[] pixels;
+	return success;
+}
+
+template< typename Real , unsigned int Dim >
+void WriteGrid( const Real *values , int res , const char *fileName )
+{
+	int resolution = 1;
+	for( int d=0 ; d<Dim ; d++ ) resolution *= res;
+
+	char *ext = GetFileExtension( fileName );
+	FILE *fp = fopen( fileName , "wb" );
+	if( !fp ) ERROR_OUT( "Failed to open voxel file for writing: %s" , fileName );
+	else
+	{
+		fwrite( &res , sizeof(int) , 1 , fp );
+		if( typeid(Real)==typeid(float) ) fwrite( values , sizeof(float) , resolution , fp );
+		else
+		{
+			float *fValues = new float[resolution];
+			for( int i=0 ; i<resolution ; i++ ) fValues[i] = float( values[i] );
+			fwrite( fValues , sizeof(float) , resolution , fp );
+			delete[] fValues;
+		}
+		fclose( fp );
+		DeletePointer( values );
+	}
+	delete[] ext;
+}
+
 template< class Real , typename ... SampleData , unsigned int ... FEMSigs >
 void Execute( int argc , char* argv[] , UIntPack< FEMSigs ... > )
 {
@@ -602,28 +667,16 @@ void Execute( int argc , char* argv[] , UIntPack< FEMSigs ... > )
 
 	if( VoxelGrid.set )
 	{
-		FILE* fp = fopen( VoxelGrid.value , "wb" );
-		if( !fp ) WARN( "Failed to open voxel file for writing: %s" , VoxelGrid.value );
-		else
-		{
-			int res = 0;
-			profiler.start();
-			Pointer( Real ) values = tree.template regularGridEvaluate< true >( solution , res , -1 , PrimalVoxel.set );
+		int res = 0;
+		profiler.start();
+		Pointer( Real ) values = tree.template regularGridEvaluate< true >( solution , res , -1 , PrimalVoxel.set );
+		int resolution = 1;
+		for( int d=0 ; d<Dim ; d++ ) resolution *= res;
 #pragma omp parallel for
-			for( int i=0 ; i<res*res*res ; i++ ) values[i] -= isoValue;
-			profiler.dumpOutput( "Got voxel grid:" );
-			fwrite( &res , sizeof(int) , 1 , fp );
-			if( typeid(Real)==typeid(float) ) fwrite( values , sizeof(float) , res*res*res , fp );
-			else
-			{
-				float *fValues = new float[res*res*res];
-				for( int i=0 ; i<res*res*res ; i++ ) fValues[i] = float( values[i] );
-				fwrite( fValues , sizeof(float) , res*res*res , fp );
-				delete[] fValues;
-			}
-			fclose( fp );
-			DeletePointer( values );
-		}
+		for( int i=0 ; i<resolution ; i++ ) values[i] -= isoValue;
+		profiler.dumpOutput( "Got voxel grid:" );
+		if( !WriteImage< Real , DIMENSION >( values , res , VoxelGrid.value , Verbose.set ) ) WriteGrid< Real , DIMENSION >( values , res , VoxelGrid.value );
+		DeletePointer( values );
 	}
 
 	if( Out.set )
diff --git a/Src/SparseMatrix.h b/Src/SparseMatrix.h
index bc157e1..8c6bbf8 100644
--- a/Src/SparseMatrix.h
+++ b/Src/SparseMatrix.h
@@ -31,9 +31,11 @@ DAMAGE.
 #include "SparseMatrixInterface.h"
 #include "Array.h"
 
-template< class T , class IndexType=int > class SparseMatrix : public SparseMatrixInterface< T , ConstPointer( MatrixEntry< T , IndexType > ) >
+template< class T , class IndexType , size_t MaxRowSize=0 > class SparseMatrix;
+
+template< class T , class IndexType > class SparseMatrix< T , IndexType , 0 > : public SparseMatrixInterface< T , ConstPointer( MatrixEntry< T , IndexType > ) >
 {
-	template< class T2 , class IndexType2 > friend class SparseMatrix;
+	template< class T2 , class IndexType2 , size_t MaxRowSize2 > friend class SparseMatrix;
 	Pointer( Pointer( MatrixEntry< T , IndexType > ) ) _entries;
 public:
 	static void Swap( SparseMatrix& M1 , SparseMatrix& M2 )
@@ -52,17 +54,17 @@ public:
 	SparseMatrix( const SparseMatrix& M );
 	SparseMatrix( SparseMatrix&& M );
 	template< class T2 , class IndexType2 >
-	SparseMatrix( const SparseMatrix< T2 , IndexType2 >& M );
+	SparseMatrix( const SparseMatrix< T2 , IndexType2 , 0 >& M );
 	~SparseMatrix();
 	SparseMatrix& operator = ( SparseMatrix&& M );
 	SparseMatrix< T , IndexType >& operator = ( const SparseMatrix< T , IndexType >& M );
 	template< class T2 , class IndexType2 >
-	SparseMatrix< T , IndexType >& operator = ( const SparseMatrix< T2 , IndexType2 >& M );
+	SparseMatrix< T , IndexType , 0 >& operator = ( const SparseMatrix< T2 , IndexType2 , 0 >& M );
 
 	template< class T2 > void operator()( const T2* in , T2* out ) const;
 
 	template< class T2 , class IndexType2 >
-	SparseMatrix< T , IndexType >& copy( const SparseMatrix< T2 , IndexType2 >& M );
+	SparseMatrix< T , IndexType , 0 >& copy( const SparseMatrix< T2 , IndexType2 , 0 >& M );
 
 	inline ConstPointer( MatrixEntry< T , IndexType > ) begin( size_t row ) const { return _entries[row]; }
 	inline ConstPointer( MatrixEntry< T , IndexType > ) end  ( size_t row ) const { return _entries[row] + (unsigned long long)rowSizes[row]; }
@@ -100,5 +102,57 @@ public:
 	template< class const_iterator >
 	static SparseMatrix Transpose( const SparseMatrixInterface< T , const_iterator >& At , size_t outRows , T (*TransposeFunction)( const T& )=NULL );
 };
+
+template< class T , class IndexType , size_t MaxRowSize > class SparseMatrix : public SparseMatrixInterface< T , ConstPointer( MatrixEntry< T , IndexType > ) >
+{
+	template< class T2 , class IndexType2 > friend class _SparseMatrix;
+	Pointer( MatrixEntry< T , IndexType > ) _entries;
+	size_t _rowNum;
+	Pointer( size_t ) _rowSizes;
+	size_t _maxRows;
+public:
+	static void Swap( SparseMatrix& M1 , SparseMatrix& M2 )
+	{
+		std::swap( M1._rowNum , M2._rowNum );
+		std::swap( M1._rowSizes , M2._rowSizes );
+		std::swap( M1._entries , M2._entries );
+	}
+	typedef SparseMatrixInterface< T , ConstPointer( MatrixEntry< T , IndexType > ) > Interface;
+	typedef ConstPointer( MatrixEntry< T , IndexType > ) RowIterator;
+
+	SparseMatrix( void );
+	SparseMatrix( const SparseMatrix& M );
+	SparseMatrix( SparseMatrix&& M );
+	template< class T2 , class IndexType2 >
+	SparseMatrix( const SparseMatrix< T2 , IndexType2 , MaxRowSize >& M );
+	SparseMatrix& operator = ( SparseMatrix&& M );
+	SparseMatrix< T , IndexType , MaxRowSize >& operator = ( const SparseMatrix< T , IndexType , MaxRowSize >& M );
+	template< class T2 , class IndexType2 >
+	SparseMatrix< T , IndexType , MaxRowSize >& operator = ( const SparseMatrix< T2 , IndexType2 , MaxRowSize >& M );
+	~SparseMatrix( void );
+
+	template< class T2 > void operator()( const T2* in , T2* out ) const;
+
+	inline ConstPointer( MatrixEntry< T , IndexType > ) begin( size_t row ) const { return _entries + MaxRowSize * row; }
+	inline ConstPointer( MatrixEntry< T , IndexType > ) end  ( size_t row ) const { return _entries + MaxRowSize * row + (unsigned long long)_rowSizes[row]; }
+	inline size_t rows                              ( void )       const { return _rowNum; }
+	inline size_t rowSize                           ( size_t idx ) const { return _rowSizes[idx]; }
+
+	SparseMatrix( size_t rowNum );
+	void resize	( size_t rowNum );
+	void setRowSize( size_t row , size_t rowSize );
+	void resetRowSize( size_t row , size_t rowSize );
+	inline      Pointer( MatrixEntry< T , IndexType > ) operator[] ( size_t idx )       { return _entries + MaxRowSize * idx; }
+	inline ConstPointer( MatrixEntry< T , IndexType > ) operator[] ( size_t idx ) const { return _entries + MaxRowSize * idx; }
+
+	// With copy move, these should be well-behaved from a memory perspective
+	SparseMatrix  operator *  ( T s ) const;
+	SparseMatrix  operator /  ( T s ) const;
+	SparseMatrix& operator *= ( T s );
+	SparseMatrix& operator /= ( T s );
+
+	Pointer( T ) operator * ( const Pointer( T ) in ) const;
+};
+
 #include "SparseMatrix.inl"
 #endif /* __SPARSEMATRIX_HPP */
diff --git a/Src/SparseMatrix.inl b/Src/SparseMatrix.inl
index 44d2b3d..d912c1d 100644
--- a/Src/SparseMatrix.inl
+++ b/Src/SparseMatrix.inl
@@ -30,15 +30,11 @@ DAMAGE.
 #include <complex>
 #include <unordered_map>
 
-///////////////////
-//  SparseMatrix //
-///////////////////
-///////////////////////////////////////
-// SparseMatrix Methods and Memebers //
-///////////////////////////////////////
-
+///////////////////////////////////////////////////////////////
+//  SparseMatrix (unconstrained max row size specialization) //
+///////////////////////////////////////////////////////////////
 template< class T , class IndexType >
-SparseMatrix< T , IndexType >::SparseMatrix( void )
+SparseMatrix< T , IndexType , 0 >::SparseMatrix( void )
 {
 	rowSizes = NullPointer( size_t );
 	rowNum = 0;
@@ -46,7 +42,7 @@ SparseMatrix< T , IndexType >::SparseMatrix( void )
 }
 
 template< class T , class IndexType >
-SparseMatrix< T , IndexType >::SparseMatrix( size_t rowNum )
+SparseMatrix< T , IndexType , 0 >::SparseMatrix( size_t rowNum )
 {
 	this->rowNum = 0;
 	rowSizes = NullPointer( size_t );
@@ -54,7 +50,7 @@ SparseMatrix< T , IndexType >::SparseMatrix( size_t rowNum )
 	resize( rowNum );
 }
 template< class T , class IndexType >
-SparseMatrix< T , IndexType >::SparseMatrix( const SparseMatrix& M )
+SparseMatrix< T , IndexType , 0 >::SparseMatrix( const SparseMatrix& M )
 {
 	rowSizes = NullPointer( size_t );
 	rowNum = 0;
@@ -67,7 +63,7 @@ SparseMatrix< T , IndexType >::SparseMatrix( const SparseMatrix& M )
 	}
 }
 template< class T , class IndexType >
-SparseMatrix< T , IndexType >::SparseMatrix( SparseMatrix&& M )
+SparseMatrix< T , IndexType , 0 >::SparseMatrix( SparseMatrix&& M )
 {
 	rowSizes = NullPointer( size_t );
 	rowNum = 0;
@@ -77,7 +73,7 @@ SparseMatrix< T , IndexType >::SparseMatrix( SparseMatrix&& M )
 }
 template< class T , class IndexType >
 template< class T2 , class IndexType2 >
-SparseMatrix< T , IndexType >::SparseMatrix( const SparseMatrix< T2 , IndexType2 >& M )
+SparseMatrix< T , IndexType , 0 >::SparseMatrix( const SparseMatrix< T2 , IndexType2 , 0 >& M )
 {
 	rowSizes = NullPointer( size_t );
 	rowNum = 0;
@@ -92,7 +88,7 @@ SparseMatrix< T , IndexType >::SparseMatrix( const SparseMatrix< T2 , IndexType2
 
 template< class T , class IndexType >
 template< class T2 , class IndexType2 >
-SparseMatrix< T , IndexType >& SparseMatrix< T , IndexType >::copy( const SparseMatrix< T2 , IndexType2 >& M  )
+SparseMatrix< T , IndexType , 0 >& SparseMatrix< T , IndexType , 0 >::copy( const SparseMatrix< T2 , IndexType2 , 0 >& M  )
 {
 	resize( M.rowNum );
 	for ( int i=0 ; i<rowNum ; i++)
@@ -107,14 +103,14 @@ SparseMatrix< T , IndexType >& SparseMatrix< T , IndexType >::copy( const Sparse
 	return *this;
 }
 template< class T , class IndexType >
-SparseMatrix< T , IndexType >& SparseMatrix< T , IndexType >::operator = ( SparseMatrix< T , IndexType >&& M )
+SparseMatrix< T , IndexType , 0 >& SparseMatrix< T , IndexType , 0 >::operator = ( SparseMatrix< T , IndexType , 0 >&& M )
 {
 	Swap( *this , M );
 	return *this;
 }
 
 template< class T , class IndexType >
-SparseMatrix< T , IndexType >& SparseMatrix< T , IndexType >::operator = ( const SparseMatrix< T , IndexType >& M )
+SparseMatrix< T , IndexType , 0 >& SparseMatrix< T , IndexType , 0 >::operator = ( const SparseMatrix< T , IndexType , 0 >& M )
 {
 	resize( M.rowNum );
 	for( int i=0 ; i<rowNum ; i++ )
@@ -126,7 +122,7 @@ SparseMatrix< T , IndexType >& SparseMatrix< T , IndexType >::operator = ( const
 }
 template< class T , class IndexType >
 template< class T2 , class IndexType2 >
-SparseMatrix< T , IndexType >& SparseMatrix< T , IndexType >::operator = (const SparseMatrix< T2 , IndexType2 >& M)
+SparseMatrix< T , IndexType , 0 >& SparseMatrix< T , IndexType , 0 >::operator = (const SparseMatrix< T2 , IndexType2 , 0 >& M)
 {
 	resize( M.rowNum );
 	for( int i=0 ; i<rowNum ; i++ )
@@ -139,13 +135,12 @@ SparseMatrix< T , IndexType >& SparseMatrix< T , IndexType >::operator = (const
 
 template< class T , class IndexType >
 template< class T2 >
-void SparseMatrix< T , IndexType >::operator() ( const T2* in , T2* out ) const { Interface::multiply( in , out ); }
+void SparseMatrix< T , IndexType , 0 >::operator() ( const T2* in , T2* out ) const { Interface::multiply( in , out ); }
 
-
-template< class T , class IndexType > SparseMatrix< T , IndexType >::~SparseMatrix( void ) { resize( 0 ); }
+template< class T , class IndexType > SparseMatrix< T , IndexType , 0 >::~SparseMatrix( void ) { resize( 0 ); }
 
 template< class T , class IndexType >
-void SparseMatrix< T , IndexType >::resize( size_t r )
+void SparseMatrix< T , IndexType , 0 >::resize( size_t r )
 {
 	if( rowNum>0 )
 	{
@@ -163,7 +158,7 @@ void SparseMatrix< T , IndexType >::resize( size_t r )
 }
 
 template< class T , class IndexType >
-void SparseMatrix< T , IndexType >::setRowSize( size_t row , size_t count )
+void SparseMatrix< T , IndexType , 0 >::setRowSize( size_t row , size_t count )
 {
 	if( row>=0 && row<rowNum )
 	{
@@ -178,7 +173,7 @@ void SparseMatrix< T , IndexType >::setRowSize( size_t row , size_t count )
 	else ERROR_OUT( "Row is out of bounds: 0 <= %d < %d" , (int)row , (int)rowNum );
 }
 template< class T , class IndexType >
-void SparseMatrix< T , IndexType >::resetRowSize( size_t row , size_t count )
+void SparseMatrix< T , IndexType , 0 >::resetRowSize( size_t row , size_t count )
 {
 	if( row>=0 && row<rowNum )
 	{
@@ -191,7 +186,7 @@ void SparseMatrix< T , IndexType >::resetRowSize( size_t row , size_t count )
 }
 
 template< class T , class IndexType >
-SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::Identity( size_t dim )
+SparseMatrix< T , IndexType , 0 > SparseMatrix< T , IndexType , 0 >::Identity( size_t dim )
 {
 	SparseMatrix I;
 	I.resize( dim );
@@ -199,52 +194,52 @@ SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::Identity( size_t di
 	return I;
 }
 template< class T , class IndexType >
-SparseMatrix< T , IndexType >& SparseMatrix< T , IndexType >::operator *= ( T s )
+SparseMatrix< T , IndexType , 0 >& SparseMatrix< T , IndexType , 0 >::operator *= ( T s )
 {
 #pragma omp parallel for
 	for( int i=0 ; i<rowNum ; i++ ) for( int j=0 ; j<rowSizes[i] ; j++ ) _entries[i][j].Value *= s;
 	return *this;
 }
 template< class T , class IndexType >
-SparseMatrix< T , IndexType >& SparseMatrix< T , IndexType >::operator /= ( T s ){ return (*this) * ( (T)1./s ); }
+SparseMatrix< T , IndexType , 0 >& SparseMatrix< T , IndexType , 0 >::operator /= ( T s ){ return (*this) * ( (T)1./s ); }
 template< class T , class IndexType >
-SparseMatrix< T , IndexType >& SparseMatrix< T , IndexType >::operator *= ( const SparseMatrix< T , IndexType >& B )
+SparseMatrix< T , IndexType , 0 >& SparseMatrix< T , IndexType , 0 >::operator *= ( const SparseMatrix< T , IndexType , 0 >& B )
 {
 	SparseMatrix foo = (*this) * B;
 	(*this) = foo;
 	return *this;
 }
 template< class T , class IndexType >
-SparseMatrix< T , IndexType >& SparseMatrix< T , IndexType >::operator += ( const SparseMatrix< T , IndexType >& B )
+SparseMatrix< T , IndexType , 0 >& SparseMatrix< T , IndexType , 0 >::operator += ( const SparseMatrix< T , IndexType , 0 >& B )
 {
 	SparseMatrix foo = (*this) + B;
 	(*this) = foo;
 	return *this;
 }
 template< class T , class IndexType >
-SparseMatrix< T , IndexType >& SparseMatrix< T , IndexType >::operator -= ( const SparseMatrix< T , IndexType >& B )
+SparseMatrix< T , IndexType , 0 >& SparseMatrix< T , IndexType , 0 >::operator -= ( const SparseMatrix< T , IndexType , 0 >& B )
 {
 	SparseMatrix foo = (*this) - B;
 	(*this) = foo;
 	return *this;
 }
 template< class T , class IndexType >
-SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::operator * ( T s ) const
+SparseMatrix< T , IndexType , 0 > SparseMatrix< T , IndexType , 0 >::operator * ( T s ) const
 {
 	SparseMatrix out = (*this);
 	return out *= s;
 }
 template< class T , class IndexType >
-SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::operator / ( T s ) const { return (*this) * ( (T)1. / s ); }
+SparseMatrix< T , IndexType , 0 > SparseMatrix< T , IndexType , 0 >::operator / ( T s ) const { return (*this) * ( (T)1. / s ); }
 template< class T , class IndexType >
-Pointer( T ) SparseMatrix< T , IndexType >::operator * ( const Pointer( T ) in ) const
+Pointer( T ) SparseMatrix< T , IndexType , 0 >::operator * ( const Pointer( T ) in ) const
 {
 	Pointer( T ) out = AllocPointer< T >( rowNum );
 	MultiplyParallel( in , out , omp_get_num_procs() , 0 );
 	return out;
 }
 template< class T , class IndexType >
-SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::operator * ( const SparseMatrix< T , IndexType >& B ) const
+SparseMatrix< T , IndexType , 0 > SparseMatrix< T , IndexType , 0 >::operator * ( const SparseMatrix< T , IndexType , 0 >& B ) const
 {
 	SparseMatrix out;
 	const SparseMatrix& A = *this;
@@ -279,7 +274,7 @@ SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::operator * ( const
 	return out;
 }
 template< class T , class IndexType >
-SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::operator + ( const SparseMatrix< T , IndexType >& B ) const
+SparseMatrix< T , IndexType , 0 > SparseMatrix< T , IndexType , 0 >::operator + ( const SparseMatrix< T , IndexType , 0 >& B ) const
 {
 	const SparseMatrix& A = *this;
 	size_t rowNum = std::max< size_t >( A.rowNum , B.rowNum );
@@ -313,7 +308,7 @@ SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::operator + ( const
 	return out;
 }
 template< class T , class IndexType >
-SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::operator - ( const SparseMatrix< T , IndexType >& B ) const
+SparseMatrix< T , IndexType , 0 > SparseMatrix< T , IndexType , 0 >::operator - ( const SparseMatrix< T , IndexType , 0 >& B ) const
 {
 	const SparseMatrix& A = *this;
 	size_t rowNum = std::max< size_t >( A.rowNum , B.rowNum );
@@ -348,7 +343,7 @@ SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::operator - ( const
 }
 
 template< class T , class IndexType >
-SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::transpose( T (*TransposeFunction)( const T& ) ) const
+SparseMatrix< T , IndexType , 0 > SparseMatrix< T , IndexType , 0 >::transpose( T (*TransposeFunction)( const T& ) ) const
 {
 	SparseMatrix A;
 	const SparseMatrix& At = *this;
@@ -382,7 +377,7 @@ SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::transpose( T (*Tran
 	return A;
 }
 template< class T , class IndexType >
-SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::transpose( size_t aRows , T (*TransposeFunction)( const T& ) ) const
+SparseMatrix< T , IndexType , 0 > SparseMatrix< T , IndexType , 0 >::transpose( size_t aRows , T (*TransposeFunction)( const T& ) ) const
 {
 	SparseMatrix A;
 	const SparseMatrix& At = *this;
@@ -421,7 +416,7 @@ SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::transpose( size_t a
 
 template< class T , class IndexType >
 template< class A_const_iterator , class B_const_iterator >
-SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::Multiply( const SparseMatrixInterface< T , A_const_iterator >& A , const SparseMatrixInterface< T , B_const_iterator >& B )
+SparseMatrix< T , IndexType , 0 > SparseMatrix< T , IndexType , 0 >::Multiply( const SparseMatrixInterface< T , A_const_iterator >& A , const SparseMatrixInterface< T , B_const_iterator >& B )
 {
 	SparseMatrix M;
 	size_t aCols = 0 , aRows = A.rows();
@@ -458,9 +453,9 @@ SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::Multiply( const Spa
 }
 template< class T , class IndexType >
 template< class const_iterator >
-SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::Transpose( const SparseMatrixInterface< T , const_iterator >& At , T (*TransposeFunction)( const T& ) )
+SparseMatrix< T , IndexType , 0 > SparseMatrix< T , IndexType , 0 >::Transpose( const SparseMatrixInterface< T , const_iterator >& At , T (*TransposeFunction)( const T& ) )
 {
-	SparseMatrix< T , IndexType > A;
+	SparseMatrix< T , IndexType , 0 > A;
 	size_t aRows = 0 , aCols = At.rows();
 	for( size_t i=0 ; i<At.rows() ; i++ ) for( const_iterator iter=At.begin(i) ; iter!=At.end(i) ; iter++ ) if( aRows<=iter->N ) aRows = iter->N+1;
 
@@ -490,9 +485,9 @@ SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::Transpose( const Sp
 }
 template< class T , class IndexType >
 template< class const_iterator >
-SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::Transpose( const SparseMatrixInterface< T , const_iterator >& At , size_t outRows , T (*TransposeFunction)( const T& ) )
+SparseMatrix< T , IndexType , 0 > SparseMatrix< T , IndexType , 0 >::Transpose( const SparseMatrixInterface< T , const_iterator >& At , size_t outRows , T (*TransposeFunction)( const T& ) )
 {
-	SparseMatrix< T , IndexType > A;
+	SparseMatrix< T , IndexType , 0 > A;
 	size_t _aRows = 0 , aCols = At.rows() , aRows = outRows;
 	for( size_t i=0 ; i<At.rows() ; i++ ) for( const_iterator iter=At.begin(i) ; iter!=At.end(i) ; iter++ ) if( aCols<=iter->N ) _aRows = iter->N+1;
 	if( _aRows>aRows ) ERROR_OUT( "Prescribed output dimension too low: %d < %zu" , aRows , _aRows );
@@ -521,3 +516,148 @@ SparseMatrix< T , IndexType > SparseMatrix< T , IndexType >::Transpose( const Sp
 		}
 	return true;
 }
+
+///////////////////////////////////////////
+//  SparseMatrix (bounded max row size ) //
+///////////////////////////////////////////
+
+template< class T , class IndexType , size_t MaxRowSize >
+SparseMatrix< T , IndexType , MaxRowSize >::SparseMatrix( void )
+{
+	_rowSizes = NullPointer( size_t );
+	_rowNum = 0;
+	_entries = NullPointer( MatrixEntry< T , IndexType > );
+	_maxRows = 0;
+}
+
+template< class T , class IndexType , size_t MaxRowSize >
+SparseMatrix< T , IndexType , MaxRowSize >::SparseMatrix( size_t rowNum ) : SparseMatrix()
+{
+	resize( rowNum );
+}
+template< class T , class IndexType , size_t MaxRowSize >
+SparseMatrix< T , IndexType , MaxRowSize >::SparseMatrix( const SparseMatrix& M ) : SparseMatrix()
+{
+	resize( M._rowNum );
+	for( int i=0 ; i<_rowNum ; i++ )
+	{
+		_rowSizes[i] = M._rowSizes[i];
+		for( int j=0 ; j<_rowSizes[i] ; j++ ) _entries[ i + MaxRowSize*j ] = M._rowEntries[ i + MaxRowSize*j ];
+	}
+}
+template< class T , class IndexType , size_t MaxRowSize >
+SparseMatrix< T , IndexType , MaxRowSize >::SparseMatrix( SparseMatrix&& M ) : SparseMatrix()
+{
+	Swap( *this , M );
+}
+template< class T , class IndexType , size_t MaxRowSize >
+template< class T2 , class IndexType2 >
+SparseMatrix< T , IndexType , MaxRowSize >::SparseMatrix( const SparseMatrix< T2 , IndexType2 , MaxRowSize >& M ) : SparseMatrix()
+{
+	resize( M._rowNum );
+	for( int i=0 ; i<_rowNum ; i++ )
+	{
+		_rowSizes[i] = M._rowSizes[i];
+		for( int j=0 ; j<_rowSizes[i] ; j++ ) _entries[ i + MaxRowSize*j ] = MatrixEntry< T , IndexType >( M._rowEntries[i][j].N , T( M._entries[ i + MaxRowSize*j ].Value ) );
+	}
+}
+
+template< class T , class IndexType , size_t MaxRowSize >
+SparseMatrix< T , IndexType , MaxRowSize >& SparseMatrix< T , IndexType , MaxRowSize >::operator = ( SparseMatrix< T , IndexType , MaxRowSize >&& M )
+{
+	Swap( *this , M );
+	return *this;
+}
+
+template< class T , class IndexType , size_t MaxRowSize >
+SparseMatrix< T , IndexType , MaxRowSize >& SparseMatrix< T , IndexType , MaxRowSize >::operator = ( const SparseMatrix< T , IndexType , MaxRowSize >& M )
+{
+	resize( M._rowNum );
+	for( int i=0 ; i<_rowNum ; i++ )
+	{
+		_rowSizes[i] = M._rowSizes[i];
+		for( int j=0 ; j<_rowSizes[i] ; j++ ) _entries[ i + MaxRowSize*j ] = M._entries[ i + MaxRowSize*j ];
+	}
+	return *this;
+}
+
+template< class T , class IndexType , size_t MaxRowSize >
+template< class T2 , class IndexType2 >
+SparseMatrix< T , IndexType , MaxRowSize >& SparseMatrix< T , IndexType , MaxRowSize >::operator = ( const SparseMatrix< T2 , IndexType2 , MaxRowSize >& M )
+{
+	resize( M._rowNum );
+	for( int i=0 ; i<_rowNum ; i++ )
+	{
+		_rowSizes[i] = M._rowSizes[i];
+		for( int j=0 ; j<_rowSizes[i] ; j++ ) _entries[ i + MaxRowSize*j ] = MatrixEntry< T , IndexType >( M._entries[ i + MaxRowSize*j ].N , T( M._entries[ i + MaxRowSize*j ].Value ) );
+	}
+	return *this;
+}
+
+template< class T , class IndexType , size_t MaxRowSize >
+template< class T2 >
+void SparseMatrix< T , IndexType , MaxRowSize >::operator() ( const T2* in , T2* out ) const { Interface::multiply( in , out ); }
+
+template< class T , class IndexType , size_t MaxRowSize > SparseMatrix< T , IndexType , MaxRowSize >::~SparseMatrix( void ) { resize( 0 ); }
+
+template< class T , class IndexType , size_t MaxRowSize >
+void SparseMatrix< T , IndexType , MaxRowSize >::resize( size_t rowNum )
+{
+	_rowNum = rowNum;
+	if( !rowNum || rowNum>_maxRows )
+	{
+		FreePointer( _rowSizes );
+		FreePointer( _entries );
+
+		if( rowNum )
+		{
+			_rowSizes = AllocPointer< size_t >( rowNum ) , memset( _rowSizes , 0 , sizeof(size_t)*rowNum );
+			_entries = AllocPointer< MatrixEntry< T , IndexType > >( rowNum * MaxRowSize );
+			_maxRows = rowNum;
+		}
+	}
+}
+
+template< class T , class IndexType , size_t MaxRowSize >
+void SparseMatrix< T , IndexType , MaxRowSize >::setRowSize( size_t row , size_t rowSize )
+{
+	if( row>=_rowNum ) ERROR_OUT( "Row is out of bounds: 0 <= %d < %d" , (int)row , (int)_rowNum );
+	else if( rowSize>MaxRowSize ) ERROR_OUT( "Row size larger than max row size: %d < %d" , (int)rowSize , (int)MaxRowSize );
+	else _rowSizes[row] = rowSize;
+}
+template< class T , class IndexType , size_t MaxRowSize >
+void SparseMatrix< T , IndexType , MaxRowSize >::resetRowSize( size_t row , size_t rowSize )
+{
+	if( row>=_rowNum ) ERROR_OUT( "Row is out of bounds: 0 <= %d < %d" , (int)row , (int)_rowNum );
+	else if( rowSize>MaxRowSize ) ERROR_OUT( "Row size larger than max row size: %d < %d" , (int)rowSize , (int)MaxRowSize );
+	else _rowSizes[row] = rowSize;
+}
+
+template< class T , class IndexType , size_t MaxRowSize >
+SparseMatrix< T , IndexType , MaxRowSize >& SparseMatrix< T , IndexType , MaxRowSize >::operator *= ( T s )
+{
+#pragma omp parallel for
+	for( int i=0 ; i<_rowNum*MaxRowSize ; i++ ) _entries[i].Value *= s;
+	return *this;
+}
+
+template< class T , class IndexType , size_t MaxRowSize >
+SparseMatrix< T , IndexType , MaxRowSize >& SparseMatrix< T , IndexType , MaxRowSize >::operator /= ( T s ){ return (*this) * ( (T)1./s ); }
+
+template< class T , class IndexType , size_t MaxRowSize >
+SparseMatrix< T , IndexType , MaxRowSize > SparseMatrix< T , IndexType , MaxRowSize >::operator * ( T s ) const
+{
+	SparseMatrix out = (*this);
+	return out *= s;
+}
+
+template< class T , class IndexType , size_t MaxRowSize >
+SparseMatrix< T , IndexType , MaxRowSize > SparseMatrix< T , IndexType , MaxRowSize >::operator / ( T s ) const { return (*this) * ( (T)1. / s ); }
+
+template< class T , class IndexType , size_t MaxRowSize >
+Pointer( T ) SparseMatrix< T , IndexType , MaxRowSize >::operator * ( const Pointer( T ) in ) const
+{
+	Pointer( T ) out = AllocPointer< T >( _rowNum );
+	MultiplyParallel( in , out , omp_get_num_procs() , 0 );
+	return out;
+}