Index: /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp	(revision 12420)
+++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp	(revision 12421)
@@ -56,6 +56,29 @@
 		observations->Variomap(predictions,x_interp,n_interp);
 	}
-	else if(strcmp(output,"prediction")==0){
-
+	else if(strcmp(output,"delaunay")==0){
+		int nobs,nel;
+		double *x     = NULL;
+		double *y     = NULL;
+		double *data  = NULL;
+		int    *index = NULL;
+		double *indexd= NULL;
+
+		observations->ObservationList(&x,&y,&data,&nobs);
+
+		printf("Generation Delaunay Triangulation\n");
+		BamgTriangulatex(&index,&nel,x,y,nobs);
+		indexd =(double*)xcalloc(nel*3,sizeof(double));
+		for(int i=0;i<nel*3;i++) indexd[i]=(double)index[i];
+		xfree((void**)&index);
+
+		printf("Interpolating\n");
+		xfree((void**)&predictions);
+		InterpFromMeshToMesh2dx(&predictions,indexd,x,y,nobs,nel,data,nobs,1,x_interp,y_interp,n_interp,NULL,0,new DataSet());
+		xfree((void**)&x);
+		xfree((void**)&y);
+		xfree((void**)&index);
+		xfree((void**)&data);
+	}
+	else if(strcmp(output,"nearestneighbor")==0){
 		/*initialize thread parameters: */
 		gate.n_interp     = n_interp;
@@ -72,4 +95,43 @@
 
 		/*launch the thread manager with Krigingxt as a core: */
+		LaunchThread(NearestNeighbort,(void*)&gate,num);
+		printf("\r      interpolation progress:  100.00%%\n");
+		xfree((void**)&gate.percent);
+	}
+	else if(strcmp(output,"idw")==0){ //Inverse distance weighting
+		/*initialize thread parameters: */
+		gate.n_interp     = n_interp;
+		gate.x_interp     = x_interp;
+		gate.y_interp     = y_interp;
+		gate.radius       = radius;
+		gate.mindata      = mindata;
+		gate.maxdata      = maxdata;
+		gate.variogram    = variogram;
+		gate.observations = observations;
+		gate.predictions  = predictions;
+		gate.error        = error;
+		gate.percent      = (double*)xcalloc(num,sizeof(double));
+
+		/*launch the thread manager with Krigingxt as a core: */
+		LaunchThread(idwt,(void*)&gate,num);
+		printf("\r      interpolation progress:  100.00%%\n");
+		xfree((void**)&gate.percent);
+	}
+	else if(strcmp(output,"prediction")==0){
+
+		/*initialize thread parameters: */
+		gate.n_interp     = n_interp;
+		gate.x_interp     = x_interp;
+		gate.y_interp     = y_interp;
+		gate.radius       = radius;
+		gate.mindata      = mindata;
+		gate.maxdata      = maxdata;
+		gate.variogram    = variogram;
+		gate.observations = observations;
+		gate.predictions  = predictions;
+		gate.error        = error;
+		gate.percent      = (double*)xcalloc(num,sizeof(double));
+
+		/*launch the thread manager with Krigingxt as a core: */
 		LaunchThread(Krigingxt,(void*)&gate,num);
 		printf("\r      interpolation progress:  100.00%%\n");
@@ -118,15 +180,55 @@
 
 	/*Intermediaries*/
+	double        localpercent;
+
+	/*partition loop across threads: */
+	PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
+	for(int idx=i0;idx<i1;idx++){
+
+		/*Print info*/
+		percent[my_thread]=double(idx-i0)/double(i1-i0)*100.;
+		localpercent=percent[0];
+		for(int i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]);
+		printf("\r      interpolation progress: %5.2lf%%",localpercent);
+
+		/*Kriging interpolation*/
+		observations->InterpolationKriging(&predictions[idx],&error[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,variogram);
+	}
+
+	return NULL;
+}/*}}}*/
+/*FUNCTION NearestNeighbort{{{*/
+void* NearestNeighbort(void* vpthread_handle){
+
+	/*gate variables :*/
+	KrigingxThreadStruct *gate        = NULL;
+	pthread_handle       *handle      = NULL;
+	int my_thread;
+	int num_threads;
+	int i0,i1;
+
+	/*recover handle and gate: */
+	handle      = (pthread_handle*)vpthread_handle;
+	gate        = (KrigingxThreadStruct*)handle->gate;
+	my_thread   = handle->id;
+	num_threads = handle->num;
+
+	/*recover parameters :*/
+	int           n_interp     = gate->n_interp;
+	double       *x_interp     = gate->x_interp;
+	double       *y_interp     = gate->y_interp;
+	double        radius       = gate->radius;
+	int           mindata      = gate->mindata;
+	int           maxdata      = gate->maxdata;
+	Variogram    *variogram    = gate->variogram;
+	Observations *observations = gate->observations;
+	double       *predictions  = gate->predictions;
+	double       *error        = gate->error;
+	double       *percent      = gate->percent;
+
+	/*Intermediaries*/
 	int           i,j,n_obs;
-	double        numerator,denominator,ratio,localpercent;
-	double       *x            = NULL;
-	double       *y            = NULL;
-	double       *obs          = NULL;
-	double       *Gamma        = NULL;
-	double       *GinvG0       = NULL;
-	double       *Ginv1        = NULL;
-	double       *GinvZ        = NULL;
-	double       *gamma0       = NULL;
-	double       *ones         = NULL;
+	double        x,y,obs;
+	double        localpercent;
 
 	/*partition loop across threads: */
@@ -140,58 +242,75 @@
 		printf("\r      interpolation progress: %5.2lf%%",localpercent);
 
-		/*Get list of observations for current point*/
+		/*Get closest observation to current point*/
+		observations->ClosestObservation(&x,&y,&obs,x_interp[idx],y_interp[idx],radius);
+		predictions[idx] = obs; 
+	}
+
+	return NULL;
+}/*}}}*/
+/*FUNCTION idwt{{{*/
+void* idwt(void* vpthread_handle){
+
+	/*gate variables :*/
+	KrigingxThreadStruct *gate        = NULL;
+	pthread_handle       *handle      = NULL;
+	int my_thread;
+	int num_threads;
+	int i0,i1;
+
+	/*recover handle and gate: */
+	handle      = (pthread_handle*)vpthread_handle;
+	gate        = (KrigingxThreadStruct*)handle->gate;
+	my_thread   = handle->id;
+	num_threads = handle->num;
+
+	/*recover parameters :*/
+	int           n_interp     = gate->n_interp;
+	double       *x_interp     = gate->x_interp;
+	double       *y_interp     = gate->y_interp;
+	double        radius       = gate->radius;
+	int           mindata      = gate->mindata;
+	int           maxdata      = gate->maxdata;
+	Variogram    *variogram    = gate->variogram;
+	Observations *observations = gate->observations;
+	double       *predictions  = gate->predictions;
+	double       *error        = gate->error;
+	double       *percent      = gate->percent;
+
+	/*Intermediaries*/
+	int           i,j,n_obs;
+	double        localpercent;
+	double        numerator,denominator,h,weight;
+	double *x     = NULL;
+	double *y     = NULL;
+	double *obs   = NULL;
+	double  power = 2.;
+
+	/*partition loop across threads: */
+	PartitionRange(&i0,&i1,n_interp,num_threads,my_thread);
+	for(int idx=i0;idx<i1;idx++){
+
+		/*Print info*/
+		percent[my_thread]=double(idx-i0)/double(i1-i0)*100.;
+		localpercent=percent[0];
+		for(i=1;i<num_threads;i++) localpercent=min(localpercent,percent[i]);
+		printf("\r      interpolation progress: %5.2lf%%",localpercent);
+
+		/*Get closest observation to current point*/
 		observations->ObservationList(&x,&y,&obs,&n_obs,x_interp[idx],y_interp[idx],radius,maxdata);
-		if(n_obs<mindata){
-			predictions[idx] = -999.0; 
-			error[idx]       = -999.0; 
-			continue;
+		if(n_obs){
+			numerator   = 0;
+			denominator = 0;
+			for(j=0;j<n_obs;j++){
+				h      = sqrt( (x[j]-x_interp[idx])*(x[j]-x_interp[idx]) + (y[j]-y_interp[idx])*(y[j]-y_interp[idx]));
+				weight = 1./pow(h,power);
+				numerator   += weight*obs[j];
+				denominator += weight;
+			}
+			predictions[idx] = numerator/denominator; 
 		}
-
-		/*Allocate intermediary matrix and vectors*/
-
-		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++){
-				//Gamma[i*n_obs+j] = variogram->SemiVariogram(x[i]-x[j],y[i]-y[j]);
-				Gamma[i*n_obs+j] = variogram->Covariance(x[i]-x[j],y[i]-y[j]);
-				Gamma[j*n_obs+i] = Gamma[i*n_obs+j];
-			}
+		else{
+			predictions[idx] = UNDEF; 
 		}
-		for(i=0;i<n_obs;i++) ones[i]=1;
-
-		/*Get semivariogram vector associated to this location*/
-		//for(i=0;i<n_obs;i++) gamma0[i] = variogram->SemiVariogram(x[i]-x_interp[idx],y[i]-y_interp[idx]);
-		for(i=0;i<n_obs;i++) gamma0[i] = variogram->Covariance(x[i]-x_interp[idx],y[i]-y_interp[idx]);
-
-		/*Solve the three linear systems*/
-		GslSolve(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0
-		GslSolve(&Ginv1, Gamma,ones,n_obs);   // Gamma^-1 ones
-		GslSolve(&GinvZ, Gamma,obs,n_obs);    // Gamma^-1 Z
-
-		/*Prepare predictor*/
-		numerator=-1.; 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.;
-		error[idx]       = - numerator*numerator/denominator;
-		for(i=0;i<n_obs;i++) predictions[idx] += (gamma0[i]-ratio)*GinvZ[i];
-		for(i=0;i<n_obs;i++) error[idx] += gamma0[i]*GinvG0[i];
-
-		/*clean-up*/
-		xfree((void**)&x);
-		xfree((void**)&y);
-		xfree((void**)&obs);
-		xfree((void**)&Gamma);
-		xfree((void**)&gamma0);
-		xfree((void**)&ones);
-		xfree((void**)&GinvG0);
-		xfree((void**)&Ginv1);
-		xfree((void**)&GinvZ);
 	}
 
@@ -219,42 +338,2 @@
 	*pvariogram = variogram;
 }/*}}}*/
-void GslSolve(double** pX,double* A,double* B,int n){/*{{{*/
-#ifdef _HAVE_GSL_
-
-		/*GSL Matrices and vectors: */
-		int              s;
-		gsl_matrix_view  a;
-		gsl_vector_view  b;
-		gsl_vector      *x = NULL;
-		gsl_permutation *p = NULL;
-
-		/*A will be modified by LU decomposition. Use copy*/
-		double* Acopy = (double*)xmalloc(n*n*sizeof(double));
-		memcpy(Acopy,A,n*n*sizeof(double));
-
-		/*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);
-
-		/*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);
-
-		//printf ("x = \n");
-		//gsl_vector_fprintf (stdout, x, "%g");
-
-		/*Copy result*/
-		double* X = (double*)xmalloc(n*sizeof(double));
-		memcpy(X,gsl_vector_ptr(x,0),n*sizeof(double));
-
-		/*Clean up and assign output pointer*/
-		xfree((void**)&Acopy);
-		gsl_permutation_free(p);
-		gsl_vector_free(x);
-		*pX=X;
-#else
-		_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 12420)
+++ /issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.h	(revision 12421)
@@ -15,5 +15,4 @@
 int  pKrigingx(double** ppredictions,double **perror,double* x, double* y, double* observations, int n_obs,double* x_interp,double* y_interp,int n_interp,Options* options);
 void ProcessVariogram(Variogram **pvariogram,Options* options);
-void GslSolve(double** pX,double* A,double* B,int n);
 
 /*threading: */
@@ -33,3 +32,5 @@
 
 void* Krigingxt(void*);
+void* NearestNeighbort(void*);
+void* idwt(void*);
 #endif /* _KRIGINGX_H */
Index: /issm/trunk-jpl/src/c/modules/Krigingx/pKrigingx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Krigingx/pKrigingx.cpp	(revision 12420)
+++ /issm/trunk-jpl/src/c/modules/Krigingx/pKrigingx.cpp	(revision 12421)
@@ -56,78 +56,10 @@
 	else if(strcmp(output,"prediction")==0){
 
-		/*Intermediaries*/
-		int           i,j,n_obs;
-		double        numerator,denominator,ratio,localpercent;
-		double       *x            = NULL;
-		double       *y            = NULL;
-		double       *obs          = NULL;
-		double       *Gamma        = NULL;
-		double       *GinvG0       = NULL;
-		double       *Ginv1        = NULL;
-		double       *GinvZ        = NULL;
-		double       *gamma0       = NULL;
-		double       *ones         = NULL;
-
 		/*partition loop across threads: */
 		for(int idx=my_rank;idx<n_interp;idx+=num_procs){
-
-			/*Print info*/
 			_printf_(true,"\r      interpolation progress: %5.2lf%%",double(idx)/double(n_interp)*100);
-
-			/*Get list of observations for current point*/
-			observations->ObservationList(&x,&y,&obs,&n_obs,x_interp[idx],y_interp[idx],radius,maxdata);
-			if(n_obs<mindata){
-				predictions[idx] = -999.0; 
-				error[idx]       = -999.0; 
-				continue;
-			}
-
-			/*Allocate intermediary matrix and vectors*/
-
-			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++){
-					//Gamma[i*n_obs+j] = variogram->SemiVariogram(x[i]-x[j],y[i]-y[j]);
-					Gamma[i*n_obs+j] = variogram->Covariance(x[i]-x[j],y[i]-y[j]);
-					Gamma[j*n_obs+i] = Gamma[i*n_obs+j];
-				}
-			}
-			for(i=0;i<n_obs;i++) ones[i]=1;
-
-			/*Get semivariogram vector associated to this location*/
-			//for(i=0;i<n_obs;i++) gamma0[i] = variogram->SemiVariogram(x[i]-x_interp[idx],y[i]-y_interp[idx]);
-			for(i=0;i<n_obs;i++) gamma0[i] = variogram->Covariance(x[i]-x_interp[idx],y[i]-y_interp[idx]);
-
-			/*Solve the three linear systems*/
-			GslSolve(&GinvG0,Gamma,gamma0,n_obs); // Gamma^-1 gamma0
-			GslSolve(&Ginv1, Gamma,ones,n_obs);   // Gamma^-1 ones
-			GslSolve(&GinvZ, Gamma,obs,n_obs);    // Gamma^-1 Z
-
-			/*Prepare predictor*/
-			numerator=-1.; 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.;
-			error[idx]       = - numerator*numerator/denominator;
-			for(i=0;i<n_obs;i++) predictions[idx] += (gamma0[i]-ratio)*GinvZ[i];
-			for(i=0;i<n_obs;i++) error[idx] += gamma0[i]*GinvG0[i];
-
-			/*clean-up*/
-			xfree((void**)&x);
-			xfree((void**)&y);
-			xfree((void**)&obs);
-			xfree((void**)&Gamma);
-			xfree((void**)&gamma0);
-			xfree((void**)&ones);
-			xfree((void**)&GinvG0);
-			xfree((void**)&Ginv1);
-			xfree((void**)&GinvZ);
+			observations->InterpolationKriging(&predictions[idx],&error[idx],x_interp[idx],y_interp[idx],radius,mindata,maxdata,variogram);
 		}
+		_printf_(true,"\r      interpolation progress: %5.2lf%%\n",100.);
 
 #ifdef _HAVE_MPI_
@@ -139,5 +71,4 @@
 		xfree((void**)&predictions); predictions=sumpredictions;
 #endif
-
 	}
 	else{