Index: /issm/trunk-jpl/src/c/classes/kriging/Observations.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/kriging/Observations.cpp	(revision 17339)
+++ /issm/trunk-jpl/src/c/classes/kriging/Observations.cpp	(revision 17340)
@@ -354,13 +354,8 @@
 	IssmPDouble   prediction,error;
 	IssmPDouble   numerator,denominator,ratio;
-	IssmPDouble  *x            = NULL;
-	IssmPDouble  *y            = NULL;
-	IssmPDouble  *obs          = NULL;
-	IssmPDouble  *Gamma        = NULL;
-	IssmPDouble  *GinvG0       = NULL;
-	IssmPDouble  *Ginv1        = NULL;
-	IssmPDouble  *GinvZ        = NULL;
-	IssmPDouble  *gamma0       = NULL;
-	IssmPDouble  *ones         = NULL;
+	IssmPDouble  *x      = NULL;
+	IssmPDouble  *y      = NULL;
+	IssmPDouble  *obs    = NULL;
+	IssmPDouble  *Lambda = NULL;
 
 	/*Some checks*/
@@ -382,42 +377,36 @@
 
 	/*Allocate intermediary matrix and vectors*/
-	Gamma  = xNew<IssmPDouble>(n_obs*n_obs);
-	gamma0 = xNew<IssmPDouble>(n_obs);
-	ones   = xNew<IssmPDouble>(n_obs);
+	IssmPDouble* A = xNew<IssmPDouble>((n_obs+1)*(n_obs+1));
+	IssmPDouble* B = xNew<IssmPDouble>(n_obs+1);
 
 	/*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;
+			A[i*(n_obs+1)+j] = variogram->Covariance(x[i]-x[j],y[i]-y[j]);
+			A[j*(n_obs+1)+i] = A[i*(n_obs+1)+j];
+		}
+		A[i*(n_obs+1)+n_obs] = 1.;
+	}
+	for(i=0;i<n_obs;i++) A[n_obs*(n_obs+1)+i]=1.;
+	A[n_obs*(n_obs+1)+n_obs] = 0.;
 
 	/*Get semivariogram vector associated to this location*/
-	//for(i=0;i<n_obs;i++) gamma0[i] = variogram->SemiVariogram(x[i]-x_interp,y[i]-y_interp);
-	for(i=0;i<n_obs;i++) gamma0[i] = variogram->Covariance(x[i]-x_interp,y[i]-y_interp);
+	for(i=0;i<n_obs;i++) B[i] = variogram->Covariance(x[i]-x_interp,y[i]-y_interp);
+	B[n_obs] = 1.;
 
 	/*Solve the three linear systems*/
 #if _HAVE_GSL_
-	DenseGslTripleSolve(&GinvG0,&Ginv1,&GinvZ,Gamma,gamma0,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
+	DenseGslSolve(&Lambda,A,B,n_obs+1);    // Gamma^-1 Z
 #else
 	_error_("GSL is required");
 #endif
 
-	/*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;
-
+	/*Compute predictor*/
 	prediction = 0.;
-	error      = - numerator*numerator/denominator;
-	for(i=0;i<n_obs;i++) prediction += (gamma0[i]-ratio)*GinvZ[i];
-	for(i=0;i<n_obs;i++) error += gamma0[i]*GinvG0[i];
+	for(i=0;i<n_obs;i++) prediction += Lambda[i]*obs[i];
+
+	/*Compute error (GSLIB p15 eq II.14)*/
+	error = variogram->Covariance(0.,0.)*(1. - Lambda[n_obs]);;
+	for(i=0;i<n_obs;i++) error += -Lambda[i]*B[i];
 
 	/*clean-up*/
@@ -427,10 +416,7 @@
 	xDelete<IssmPDouble>(y);
 	xDelete<IssmPDouble>(obs);
-	xDelete<IssmPDouble>(Gamma);
-	xDelete<IssmPDouble>(gamma0);
-	xDelete<IssmPDouble>(ones);
-	xDelete<IssmPDouble>(GinvG0);
-	xDelete<IssmPDouble>(Ginv1);
-	xDelete<IssmPDouble>(GinvZ);
+	xDelete<IssmPDouble>(A);
+	xDelete<IssmPDouble>(B);
+	xDelete<IssmPDouble>(Lambda);
 
 }/*}}}*/