Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp	(revision 12173)
+++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp	(revision 12174)
@@ -15,8 +15,132 @@
 	double *predictions = NULL;
 
-	/*Allocate output*/
-	predictions=(double*)xmalloc(n_interp*sizeof(double));
+	/*Intermediaries*/
+	int     i,j;
+	double  numerator,denominator,ratio;
+	double *Gamma     = NULL;
+	double *GinvG0    = NULL;
+	double *Ginv1     = NULL;
+	double *GinvZ     = NULL;
+	double *gamma0    = NULL;
+	double *ones      = NULL;
 
-	/*Assign output pointer*/
+	/*Memory allocation*/
+	predictions =(double*)xmalloc(n_interp*sizeof(double));
+	Gamma       =(double*)xmalloc(n_obs*n_obs*sizeof(double));
+	gamma0      =(double*)xmalloc(n_obs*sizeof(double));
+	ones        =(double*)xmalloc(n_obs*sizeof(double));
+
+	/*First: Create semivariogram matrix for observations*/
+	for(i=0;i<n_obs;i++){
+		for(j=0;j<=i;j++){
+			SemiVariogram(&Gamma[i*n_obs+j],x[i],y[i],x[j],y[j]);
+			Gamma[j*n_obs+i] = Gamma[i*n_obs+j];
+		}
+	}
+	for(i=0;i<n_obs;i++) ones[i]=1;
+
+	/*Loop over all interpolations*/
+	for(int idx=0;idx<n_interp;idx++){
+
+		/*Get semivariogram vector associated to this location*/
+		for(i=0;i<n_obs;i++) SemiVariogram(&gamma0[i],x[i],y[i],x_interp[idx],y_interp[idx]);
+
+		/*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
+
+		/*Prepare predictor*/
+		numerator=0.; denominator=0.;
+		for(i=0;i<n_obs;i++) numerator  +=GinvG0[i];
+		for(i=0;i<n_obs;i++) denominator+=Ginv1[i];
+		ratio=numerator/denominator;
+
+		predictions[idx]=0;
+		for(i=0;i<n_obs;i++) predictions[idx] += (gamma0[i]-ratio)*GinvZ[i];
+
+		/*clean-up*/
+		xfree((void**)&GinvG0);
+		xfree((void**)&Ginv1);
+		xfree((void**)&GinvZ);
+	}
+
+	/*clean-up and Assign output pointer*/
+	xfree((void**)&Gamma);
+	xfree((void**)&gamma0);
+	xfree((void**)&ones);
 	*ppredictions=predictions;
 }
+
+void SemiVariogram(double* gamma,double x1,double y1,double x2,double y2){
+
+	/*Calculate distance*/
+	double r=sqrt(pow(x1-x2,2.)+pow(y1-y2,2.));
+
+	/*Switch between variogram models*/
+	switch(1){
+		case 1:{ /*Exponential*/
+			double c0=0.2;
+			double c1=0.8;
+			double a =1;
+			*gamma = c0 + c1*(1-exp(-r/a));
+			return;}
+		default:
+			_error_("Not implemented yet");
+	}
+}
+
+void MumpsSolve(double** Xdouble,double* Adouble,double* Bdouble,int n){
+#ifdef _HAVE_PETSC_
+
+	/*Intermediaries*/
+	Mat A = NULL;
+	Vec B = NULL;
+	Vec X = NULL;
+	KSP ksp = NULL;
+	PC  pc  = 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");
+
+	/*Create indexing and ones*/
+	int* idx=(int*)xmalloc(n*sizeof(int));
+	for(int i=0;i<n;i++) idx[i]=i;
+
+	/*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);
+
+	/*Convert vector to PETSc vector*/
+	VecCreateSeq(PETSC_COMM_SELF,n,&B);
+	VecSetValues(B,n,idx,Bdouble,INSERT_VALUES);
+	VecAssemblyBegin(B);
+	VecAssemblyEnd(B);
+
+	/*Allocate output*/
+	VecCreateSeq(PETSC_COMM_SELF,n,&X);
+
+	/*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);
+#else
+	error_("PETSc support required");
+#endif
+}
Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h	(revision 12173)
+++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h	(revision 12174)
@@ -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);
+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
 
 #endif /* _KRIGINGX_H */
-