Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp	(revision 12182)
+++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp	(revision 12183)
@@ -9,6 +9,7 @@
 #include "../../objects/objects.h"
 #include "../modules.h"
+#include <gsl/gsl_linalg.h>
 
-int Krigingx(double** ppredictions,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp){
+int Krigingx(double** ppredictions,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp,Options* options){
 
 	/*output*/
@@ -49,7 +50,7 @@
 
 		/*Solve the three linear systems*/
-		MumpsSolve(&GinvG0,Gamma,gamma0,n_obs);       // Gamma^-1 gamma0
-		MumpsSolve(&Ginv1, Gamma,ones,n_obs);         // Gamma^-1 ones
-		MumpsSolve(&GinvZ, Gamma,observations,n_obs); // Gamma^-1 Z
+		GslSolve(&GinvG0,Gamma,gamma0,n_obs);       // Gamma^-1 gamma0
+		GslSolve(&Ginv1, Gamma,ones,n_obs);         // Gamma^-1 ones
+		GslSolve(&GinvZ, Gamma,observations,n_obs); // Gamma^-1 Z
 
 		/*Prepare predictor*/
@@ -103,56 +104,42 @@
 }
 
-void MumpsSolve(double** Xdouble,double* Adouble,double* Bdouble,int n){
-#ifdef _HAVE_PETSC_
+void GslSolve(double** pX,double* A,double* B,int n){
+#ifdef _HAVE_GSL_
 
-	/*Intermediaries*/
-	Mat A = NULL;
-	Vec B = NULL;
-	Vec X = NULL;
-	KSP ksp = NULL;
-	PC  pc  = NULL;
+	/*GSL Matrices and vectors: */
+	int              s;
+	gsl_matrix_view  a;
+	gsl_vector_view  b;
+	gsl_vector      *x = NULL;
+	gsl_permutation *p = NULL;
 
-	/*Create parameters and select MUMPS solver*/
-	PetscOptionsSetFromOptions();
-	PetscOptionsClear();
-	PetscOptionsInsertMultipleString("-mat_type mpiaij \
-				-ksp_type preonly \
-				-pc_type lu \
-				-pc_factor_mat_solver_package mumps \
-				-mat_mumps_icntl_14 120 \
-				-pc_factor_shift_positive_definite true");
+	/*A will be modified by LU decomposition. Use copy*/
+	double* Acopy = (double*)xmalloc(n*n*sizeof(double));
+	memcpy(Acopy,A,n*n*sizeof(double));
 
-	/*Create indexing and ones*/
-	int* idx=(int*)xmalloc(n*sizeof(int));
-	for(int i=0;i<n;i++) idx[i]=i;
+	/*Initialize gsl matrices and vectors: */
+	a = gsl_matrix_view_array (Acopy,n,n);
+	b = gsl_vector_view_array (B,n);
+	x = gsl_vector_alloc (n);
 
-	/*Convert matrix to Dense PETSc matrix*/
-	MatCreateSeqDense(PETSC_COMM_SELF,n,n,NULL,&A);
-	MatSetValues(A,n,idx,n,idx,Adouble,INSERT_VALUES);
-	MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); 
-	MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
+	/*Run LU and solve: */
+	p = gsl_permutation_alloc (n);
+	gsl_linalg_LU_decomp (&a.matrix, p, &s);
+	gsl_linalg_LU_solve (&a.matrix, p, &b.vector, x);
 
-	/*Convert vector to PETSc vector*/
-	VecCreateSeq(PETSC_COMM_SELF,n,&B);
-	VecSetValues(B,n,idx,Bdouble,INSERT_VALUES);
-	VecAssemblyBegin(B);
-	VecAssemblyEnd(B);
+	//printf ("x = \n");
+	//gsl_vector_fprintf (stdout, x, "%g");
 
-	/*Allocate output*/
-	VecCreateSeq(PETSC_COMM_SELF,n,&X);
+	/*Copy result*/
+	double* X = (double*)xmalloc(n*sizeof(double));
+	memcpy(X,gsl_vector_ptr(x,0),n*sizeof(double));
 
-	/*Go Solve*/
-	KSPCreate(MPI_COMM_WORLD,&ksp);
-	KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);
-	KSPSetFromOptions(ksp);
-	KSPGetPC(ksp,&pc);
-	PCFactorSetMatSolverPackage(pc,MATSOLVERMUMPS);
-	KSPSolve(ksp,B,X);
-	KSPFree(&ksp);
-
-	/*Serialize vector*/
-	VecToMPISerial(Xdouble,X);
+	/*Clean up and assign output pointer*/
+	xfree((void**)&Acopy);
+	gsl_permutation_free(p);
+	gsl_vector_free(x);
+	*pX=X;
 #else
-	error_("PETSc support required");
+	error_("GSL support required");
 #endif
 }
Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h	(revision 12182)
+++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h	(revision 12183)
@@ -9,7 +9,7 @@
 #include "../../toolkits/toolkits.h"
 
-int Krigingx(double** ppredictions,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp);
+int  Krigingx(double** ppredictions,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp,Options* options);
 void SemiVariogram(double* gamma,double x1,double y1,double x2,double y2);
-void MumpsSolve(double** X,double* A,double* B,int n); //Solve X for the linear system AX = B
+void GslSolve(double** pX,double* A,double* B,int n);
 
 #endif /* _KRIGINGX_H */