Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 14434)
+++ /issm/trunk-jpl/src/c/classes/FemModel.cpp	(revision 14435)
@@ -826,5 +826,5 @@
 		for (i=0;i<this->elements->Size();i++){
 			element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
-			element->CreateKMatrix(Kff_temp,NULL,NULL);
+			element->CreateKMatrix(Kff_temp,NULL);
 		}
 		for (i=0;i<this->loads->Size();i++){
@@ -845,5 +845,5 @@
 	for (i=0;i<this->elements->Size();i++){
 		element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
-		element->CreateKMatrix(Kff,Kfs,df);
+		element->CreateKMatrix(Kff,Kfs);
 	}
 	for (i=0;i<this->loads->Size();i++){
@@ -872,4 +872,10 @@
 			if(load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pf,kmax);
 		}
+	}
+
+	/*Create dof vector for stiffness matrix preconditioning*/
+	for (i=0;i<this->elements->Size();i++){
+		element=dynamic_cast<Element*>(this->elements->GetObjectByOffset(i));
+		element->CreateDVector(df);
 	}
 
Index: /issm/trunk-jpl/src/c/classes/matrix/ElementVector.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/matrix/ElementVector.cpp	(revision 14434)
+++ /issm/trunk-jpl/src/c/classes/matrix/ElementVector.cpp	(revision 14435)
@@ -190,5 +190,5 @@
 	IssmDouble* localvalues=NULL;
 
-	if(this->fsize && pf){
+	if(this->fsize){
 		/*first, retrieve values that are in the f-set from the g-set values vector: */
 		localvalues=xNew<IssmDouble>(this->fsize);
Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Elements/Element.h	(revision 14434)
+++ /issm/trunk-jpl/src/c/classes/objects/Elements/Element.h	(revision 14435)
@@ -30,5 +30,6 @@
 		virtual void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters)=0;
 		virtual void   SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int set1_enum,int set2_enum)=0;
-		virtual void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>*  Kfs,Vector<IssmDouble>* df)=0;
+		virtual void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>*  Kfs)=0;
+		virtual void   CreateDVector(Vector<IssmDouble>* df)=0;
 		virtual void   CreatePVector(Vector<IssmDouble>* pf)=0;
 		virtual void   CreateJacobianMatrix(Matrix<IssmDouble>* Jff)=0;
Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp	(revision 14434)
+++ /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp	(revision 14435)
@@ -405,9 +405,8 @@
 /*}}}*/
 /*FUNCTION Penta::CreateKMatrix {{{*/
-void  Penta::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,Vector<IssmDouble>* df){
+void  Penta::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){
 
 	/*retrieve parameters: */
 	ElementMatrix* Ke=NULL;
-	ElementVector* De=NULL;
 	int analysis_type;
 	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
@@ -424,5 +423,5 @@
 		#ifdef _HAVE_DIAGNOSTIC_
 		case DiagnosticHorizAnalysisEnum:
-			Ke=CreateKMatrixDiagnosticHoriz(); De=CreateDVectorDiagnosticHoriz();
+			Ke=CreateKMatrixDiagnosticHoriz();
 			break;
 		case AdjointHorizAnalysisEnum:
@@ -467,4 +466,61 @@
 		delete Ke;
 	}
+}
+/*}}}*/
+/*FUNCTION Penta::CreateKMatrixPrognostic {{{*/
+ElementMatrix* Penta::CreateKMatrixPrognostic(void){
+
+	if (!IsOnBed()) return NULL;
+
+	/*Depth Averaging Vx and Vy*/
+	this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
+	this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
+
+	Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
+	ElementMatrix* Ke=tria->CreateKMatrixPrognostic();
+	delete tria->material; delete tria;
+
+	/*Delete Vx and Vy averaged*/
+	this->inputs->DeleteInput(VxAverageEnum);
+	this->inputs->DeleteInput(VyAverageEnum);
+
+	/*clean up and return*/
+	return Ke;
+}
+/*}}}*/
+/*FUNCTION Penta::CreateBasalMassMatrix{{{*/
+ElementMatrix* Penta::CreateBasalMassMatrix(void){
+
+	if (!IsOnBed()) return NULL;
+
+	Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
+	ElementMatrix* Ke=tria->CreateMassMatrix();
+	delete tria->material; delete tria;
+
+	/*clean up and return*/
+	return Ke;
+}
+/*}}}*/
+/*FUNCTION Penta::CreateDVector {{{*/
+void  Penta::CreateDVector(Vector<IssmDouble>* df){
+
+	/*retrieve parameters: */
+	ElementMatrix* Ke=NULL;
+	ElementVector* De=NULL;
+	int analysis_type;
+	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+
+	/*Checks in debugging {{{*/
+	_assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs);
+	/*}}}*/
+
+	switch(analysis_type){
+		#ifdef _HAVE_DIAGNOSTIC_
+		case DiagnosticHorizAnalysisEnum:
+			De=CreateDVectorDiagnosticHoriz();
+			break;
+		#endif
+	}
+
 	/*Add to global Vector*/
 	if(De){
@@ -472,38 +528,4 @@
 		delete De;
 	}
-}
-/*}}}*/
-/*FUNCTION Penta::CreateKMatrixPrognostic {{{*/
-ElementMatrix* Penta::CreateKMatrixPrognostic(void){
-
-	if (!IsOnBed()) return NULL;
-
-	/*Depth Averaging Vx and Vy*/
-	this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
-	this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
-
-	Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
-	ElementMatrix* Ke=tria->CreateKMatrixPrognostic();
-	delete tria->material; delete tria;
-
-	/*Delete Vx and Vy averaged*/
-	this->inputs->DeleteInput(VxAverageEnum);
-	this->inputs->DeleteInput(VyAverageEnum);
-
-	/*clean up and return*/
-	return Ke;
-}
-/*}}}*/
-/*FUNCTION Penta::CreateBasalMassMatrix{{{*/
-ElementMatrix* Penta::CreateBasalMassMatrix(void){
-
-	if (!IsOnBed()) return NULL;
-
-	Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
-	ElementMatrix* Ke=tria->CreateMassMatrix();
-	delete tria->material; delete tria;
-
-	/*clean up and return*/
-	return Ke;
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h	(revision 14434)
+++ /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h	(revision 14435)
@@ -81,5 +81,6 @@
 		void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
 		void   SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int set1_enum,int set2_enum);
-		void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,Vector<IssmDouble>* df);
+		void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
+		void   CreateDVector(Vector<IssmDouble>* df);
 		void   CreatePVector(Vector<IssmDouble>* pf);
 		void   CreateJacobianMatrix(Matrix<IssmDouble>* Jff);
Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp	(revision 14434)
+++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp	(revision 14435)
@@ -214,5 +214,5 @@
 /*}}}*/
 /*FUNCTION Tria::CreateKMatrix {{{*/
-void  Tria::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,Vector<IssmDouble>* df){
+void  Tria::CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs){
 
 	/*retreive parameters: */
@@ -561,4 +561,10 @@
 	delete gauss;
 	return Ke;
+}
+/*}}}*/
+/*FUNCTION Tria::CreateDVector {{{*/
+void  Tria::CreateDVector(Vector<IssmDouble>* df){
+
+	/*Nothing done yet*/
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h	(revision 14434)
+++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.h	(revision 14435)
@@ -77,5 +77,6 @@
 		void   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
 		void   SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int set1_enum,int set2_enum);
-		void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs,Vector<IssmDouble>* df);
+		void   CreateKMatrix(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
+		void   CreateDVector(Vector<IssmDouble>* df);
 		void   CreatePVector(Vector<IssmDouble>* pf);
 		void   CreateJacobianMatrix(Matrix<IssmDouble>* Jff);