Index: /issm/trunk/src/c/solvers/solver_adjoint_linear.cpp
===================================================================
--- /issm/trunk/src/c/solvers/solver_adjoint_linear.cpp	(revision 6000)
+++ /issm/trunk/src/c/solvers/solver_adjoint_linear.cpp	(revision 6001)
@@ -13,12 +13,23 @@
 
 	/*intermediary: */
-	Mat Kgg = NULL, Kff=NULL, Kfs=NULL;
-	Vec ug  = NULL, uf=NULL;
-	Vec pg  = NULL, pf=NULL;
+	Mat  Kgg = NULL, Kff = NULL, Kfs = NULL;
+	Vec  ug  = NULL, uf  = NULL;
+	Vec  pg  = NULL, pf  = NULL;
+	bool kffpartitioning = false;
 
-	SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+	/*Recover parameters: */
+	femmodel->parameters->FindParam(&kffpartitioning,KffEnum);
+
+	if(kffpartitioning){
+		SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		Reduceloadx(pf, Kfs, femmodel->ys,true); MatFree(&Kfs); //true means spc = 0
+	}
+	else{
+		SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg);
+		Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets,femmodel->parameters,true);VecFree(&pg); MatFree(&Kfs);//true means spc = 0
+	}
 
 	Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg);
-	Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets,femmodel->parameters,true);VecFree(&pg); MatFree(&Kfs);//true means spc = 0
 
 	Solverx(&uf, Kff, pf, NULL, femmodel->parameters); MatFree(&Kff); VecFree(&pf);
Index: /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp
===================================================================
--- /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp	(revision 6000)
+++ /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp	(revision 6001)
@@ -24,5 +24,5 @@
 	/*parameters:*/
 	int  verbose=0;
-	bool kff=false;
+	bool kffpartitioning=false;
 	int min_mechanical_constraints;
 	int max_nonlinear_iterations;
@@ -30,5 +30,5 @@
 	/*Recover parameters: */
 	femmodel->parameters->FindParam(&verbose,VerboseEnum);
-	femmodel->parameters->FindParam(&kff,KffEnum);
+	femmodel->parameters->FindParam(&kffpartitioning,KffEnum);
 	femmodel->parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum);
 	femmodel->parameters->FindParam(&max_nonlinear_iterations,MaxNonlinearIterationsEnum);
@@ -54,7 +54,7 @@
 		VecFree(&old_uf);old_uf=uf;
 
-		if(kff){
+		if(kffpartitioning){
 			SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters);
-			Reduceloadx(pf, Kfs, femmodel->ys, femmodel->parameters); MatFree(&Kfs);
+			Reduceloadx(pf, Kfs, femmodel->ys); MatFree(&Kfs);
 		}
 		else{
Index: /issm/trunk/src/c/solvers/solver_linear.cpp
===================================================================
--- /issm/trunk/src/c/solvers/solver_linear.cpp	(revision 6000)
+++ /issm/trunk/src/c/solvers/solver_linear.cpp	(revision 6001)
@@ -11,15 +11,23 @@
 
 	/*intermediary: */
-	Mat Kgg = NULL, Kff = NULL, Kfs   = NULL;
-	Vec ug  = NULL, uf  = NULL;
-	Vec pg  = NULL, pf  = NULL;
+	Mat  Kgg = NULL, Kff = NULL, Kfs   = NULL;
+	Vec  ug  = NULL, uf  = NULL;
+	Vec  pg  = NULL, pf  = NULL;
+	bool kffpartitioning;
 
-	SystemMatricesx(&Kgg,NULL, NULL, &pg,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+	/*Recover parameters: */
+	femmodel->parameters->FindParam(&kffpartitioning,KffEnum);
 
-	Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg);
-	Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets,femmodel->parameters);VecFree(&pg); MatFree(&Kfs);
+	if(kffpartitioning){
+		SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		Reduceloadx(pf, Kfs, femmodel->ys); MatFree(&Kfs);
+	}
+	else{
+		SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+		Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg);
+		Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets,femmodel->parameters);VecFree(&pg); MatFree(&Kfs);
+	}
 
 	Solverx(&uf, Kff, pf, NULL, femmodel->parameters); MatFree(&Kff); VecFree(&pf);
-
 	Mergesolutionfromftogx(&ug, uf,femmodel->ys,femmodel->nodesets,femmodel->parameters);VecFree(&uf);
 
Index: /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp
===================================================================
--- /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp	(revision 6000)
+++ /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp	(revision 6001)
@@ -34,8 +34,9 @@
 	int verbose=0;
 	bool lowmem=0;
+	bool kffpartitioning;
 
 	/*Recover parameters: */
 	kflag=1; pflag=1;
-
+	femmodel->parameters->FindParam(&kffpartitioning,KffEnum);
 	femmodel->parameters->FindParam(&verbose,VerboseEnum);
 	femmodel->parameters->FindParam(&lowmem,LowmemEnum);
@@ -50,12 +51,13 @@
 	for(;;){
 
-
-		SystemMatricesx(&Kgg,NULL, NULL, &pg,NULL, &melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
-
-		Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg);
-		Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets,femmodel->parameters);
-
-		VecFree(&pg); 
-		MatFree(&Kfs);
+		if(kffpartitioning){
+			SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf,&melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+			Reduceloadx(pf, Kfs, femmodel->ys); MatFree(&Kfs);
+		}
+		else{
+			SystemMatricesx(&Kgg,NULL, NULL, &pg,NULL, &melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+			Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg);
+			Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets,femmodel->parameters);VecFree(&pg); MatFree(&Kfs);
+		}
 
 		VecFree(&tf);