Index: /issm/trunk-jpl/src/c/classes/kriging/Observations.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/kriging/Observations.cpp	(revision 17327)
+++ /issm/trunk-jpl/src/c/classes/kriging/Observations.cpp	(revision 17328)
@@ -401,7 +401,8 @@
 	/*Solve the three linear systems*/
 #if _HAVE_GSL_
-	DenseGslSolve(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0
-	DenseGslSolve(&Ginv1, Gamma,ones,n_obs);   // Gamma^-1 ones
-	DenseGslSolve(&GinvZ, Gamma,obs,n_obs);    // Gamma^-1 Z
+	DenseGslTripleSolve(&GinvG0,&Ginv1,&GinvZ,Gamma,Gamma,ones,obs,n_obs);
+	//DenseGslSolve(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0
+	//DenseGslSolve(&Ginv1, Gamma,ones,n_obs);   // Gamma^-1 ones
+	//DenseGslSolve(&GinvZ, Gamma,obs,n_obs);    // Gamma^-1 Z
 #else
 	_error_("GSL is required");
Index: /issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp	(revision 17327)
+++ /issm/trunk-jpl/src/c/toolkits/gsl/DenseGslSolve.cpp	(revision 17328)
@@ -32,5 +32,5 @@
 }
 /*}}}*/
-void DenseGslSolve(/*output*/ IssmPDouble** px,/*stiffness matrix:*/ IssmPDouble* Kff, int Kff_M, int Kff_N, /*right hand side load vector: */ IssmPDouble* pf, int pf_M, Parameters* parameters){ /*{{{*/
+void DenseGslSolve(IssmPDouble** px,IssmPDouble* Kff,int Kff_M,int Kff_N,IssmPDouble* pf,int pf_M,Parameters* parameters){ /*{{{*/
 
 	/*Intermediary: */
@@ -44,4 +44,46 @@
 	/*allocate output pointers: */
 	*px=x;
+}
+/*}}}*/
+void DenseGslTripleSolve(IssmPDouble** pX1,IssmPDouble** pX2,IssmPDouble** pX3,IssmPDouble* A,IssmPDouble* B1,IssmPDouble* B2,IssmPDouble* B3,int n){/*{{{*/
+
+	/*Allocate output: */
+	IssmPDouble *X1 = xNew<IssmPDouble>(n);
+	IssmPDouble *X2 = xNew<IssmPDouble>(n);
+	IssmPDouble *X3 = xNew<IssmPDouble>(n);
+
+#ifdef _HAVE_GSL_
+	/*GSL Matrices and vectors: */
+	int s;
+
+	double* Acopy = xNew<double>(n*n);
+	xMemCpy(Acopy,A,n*n);
+
+	/*Initialize gsl matrices and vectors: */
+	gsl_matrix_view a = gsl_matrix_view_array(Acopy,n,n);
+	gsl_vector_view b1= gsl_vector_view_array(B1,n);
+	gsl_vector_view b2= gsl_vector_view_array(B2,n);
+	gsl_vector_view b3= gsl_vector_view_array(B3,n);
+	gsl_vector_view x1= gsl_vector_view_array(X1,n);
+	gsl_vector_view x2= gsl_vector_view_array(X2,n);
+	gsl_vector_view x3= gsl_vector_view_array(X3,n);
+
+	/*Run LU and solve: */
+	gsl_permutation* p = gsl_permutation_alloc (n);
+	gsl_linalg_LU_decomp(&a.matrix,p,&s);
+
+	gsl_linalg_LU_solve(&a.matrix,p,&b1.vector,&x1.vector);
+	gsl_linalg_LU_solve(&a.matrix,p,&b2.vector,&x2.vector);
+	gsl_linalg_LU_solve(&a.matrix,p,&b3.vector,&x3.vector);
+
+	/*Clean up and assign output pointer*/
+	xDelete(Acopy);
+	gsl_permutation_free(p);
+#endif
+
+	/*allocate output pointers: */
+	*pX1=X1;
+	*pX2=X2;
+	*pX3=X3;
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/toolkits/gsl/gslincludes.h
===================================================================
--- /issm/trunk-jpl/src/c/toolkits/gsl/gslincludes.h	(revision 17327)
+++ /issm/trunk-jpl/src/c/toolkits/gsl/gslincludes.h	(revision 17328)
@@ -21,5 +21,6 @@
 
 void DenseGslSolve(IssmPDouble** pX,IssmPDouble* A,IssmPDouble* B, int n);
-void DenseGslSolve(IssmDouble** px, IssmDouble* Kff,int Kff_M, int Kff_N, IssmDouble* pf, int pf_M, Parameters* parameters);
+void DenseGslSolve(IssmDouble** pX,IssmDouble* Kff,int Kff_M,int Kff_N,IssmDouble* pf,int pf_M,Parameters* parameters);
+void DenseGslTripleSolve(IssmPDouble** pX1,IssmPDouble** pX2,IssmPDouble** pX3,IssmPDouble* A,IssmPDouble* B1,IssmPDouble* B2,IssmPDouble* B3,int n);
 
 void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n);