Index: /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 24678)
+++ /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp	(revision 24679)
@@ -8,5 +8,5 @@
 #include "../AllocateSystemMatricesx/AllocateSystemMatricesx.h"
 
-void SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax,FemModel* femmodel){
+void SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax,FemModel* femmodel, bool isAllocated){
 
 	/*intermediary: */
@@ -64,5 +64,13 @@
 
 	/*Allocate stiffness matrices and load vector*/
-	AllocateSystemMatricesx(&Kff,&Kfs,&df,&pf,femmodel);
+	if(isAllocated) {
+		Kff  = *pKff;
+		Kfs  = *pKfs;
+		pf   = *ppf;
+		df   = *pdf;
+	}
+	else {
+		AllocateSystemMatricesx(&Kff,&Kfs,&df,&pf,femmodel);
+	}
 
 	/*Display size*/
Index: /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.h	(revision 24678)
+++ /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.h	(revision 24679)
@@ -9,5 +9,5 @@
 
 /* local prototypes: */
-void SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax,FemModel* femmodel);
+void SystemMatricesx(Matrix<IssmDouble>** pKff, Matrix<IssmDouble>** pKfs, Vector<IssmDouble>** ppf, Vector<IssmDouble>** pdf, IssmDouble* pkmax,FemModel* femmodel, bool isAllocated=false);
 
 #endif  /* _SYSTEMMATRICESX_H */
Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp
===================================================================
--- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp	(revision 24678)
+++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_nonlinear.cpp	(revision 24679)
@@ -58,4 +58,7 @@
 	InputUpdateFromSolutionx(femmodel,ug);
 
+	// allocate the matrices once and reuce them per iteration
+	AllocateSystemMatricesx(&Kff,&Kfs,&df,&pf,femmodel);
+
 	for(;;){
 
@@ -64,7 +67,7 @@
 		delete ug;
 
-		SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel);
+		SystemMatricesx(&Kff,&Kfs,&pf,&df,NULL,femmodel, true);
 		CreateNodalConstraintsx(&ys,femmodel->nodes);
-		Reduceloadx(pf, Kfs, ys); delete Kfs;
+		Reduceloadx(pf, Kfs, ys);
 		femmodel->profiler->Start(SOLVER);
 		Solverx(&uf, Kff, pf, old_uf, df, femmodel->parameters);
@@ -73,5 +76,5 @@
 		Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete ys;
 
-		convergence(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs); delete Kff; delete pf; delete df;
+		convergence(&converged,Kff,pf,uf,old_uf,eps_res,eps_rel,eps_abs);
 		InputUpdateFromConstantx(femmodel,converged,ConvergedEnum);
 		InputUpdateFromSolutionx(femmodel,ug);
@@ -102,5 +105,15 @@
 			break;
 		}
+
+		// Set the matrix entries to zero if we do an other iteration
+		Kff->SetZero();
+		Kfs->SetZero();
+		df->Set(0);
+		pf->Set(0);
 	}
+
+	// delete matrices after the iteration loop
+	delete Kff; delete pf; delete df;
+	delete Kfs;
 
 	if(VerboseConvergence()) _printf0_("\n   total number of iterations: " << count << "\n");
Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmAbsMat.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/IssmAbsMat.h	(revision 24678)
+++ /issm/trunk-jpl/src/c/toolkits/issm/IssmAbsMat.h	(revision 24679)
@@ -43,4 +43,5 @@
 		virtual void SetValues(int m,int* idxm,int n,int* idxn,doubletype* values,InsMode mode)=0;
 		virtual void Convert(MatrixType type)=0;
+		virtual void SetZero(void){};
 		#ifndef _HAVE_WRAPPERS_
 		virtual IssmAbsVec<IssmDouble>* Solve(IssmAbsVec<IssmDouble>* pf, Parameters* parameters)=0;
Index: /issm/trunk-jpl/src/c/toolkits/issm/IssmMat.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/issm/IssmMat.h	(revision 24678)
+++ /issm/trunk-jpl/src/c/toolkits/issm/IssmMat.h	(revision 24679)
@@ -239,4 +239,7 @@
 			matrix->convert(type);
 		}/*}}}*/
+		void SetZero(void){/*{{{*/
+			matrix->SetZero();
+		}/*}}}*/
 		#ifndef _HAVE_WRAPPERS_
 		IssmVec<doubletype>* Solve(IssmVec<doubletype>* pf, Parameters* parameters){ /*{{{*/
Index: /issm/trunk-jpl/src/c/toolkits/objects/Matrix.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/objects/Matrix.h	(revision 24678)
+++ /issm/trunk-jpl/src/c/toolkits/objects/Matrix.h	(revision 24679)
@@ -307,5 +307,18 @@
 		}
 		/*}}}*/
-
+		/*
+		* sets all values to 0 but keeps the structure of a sparse matrix
+		*/
+		void SetZero(void) {/*{{{*/
+			if(type==PetscMatType){
+				#ifdef _HAVE_PETSC_
+				this->pmatrix->SetZero();
+				#endif
+			}
+			else{
+				this->imatrix->SetZero();
+			}
+		}
+		/*}}}*/
 };
 
Index: /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.cpp	(revision 24678)
+++ /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.cpp	(revision 24679)
@@ -193,2 +193,6 @@
 }
 /*}}}*/
+void PetscMat::SetZero(void){/*{{{*/
+	MatZeroEntries(this->matrix);
+}
+/*}}}*/
Index: /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.h	(revision 24678)
+++ /issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscMat.h	(revision 24679)
@@ -52,4 +52,5 @@
 		void SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode);
 		void Convert(MatrixType type);
+		void SetZero(void);
 };