Index: /issm/trunk-jpl/src/c/Container/Observations.cpp
===================================================================
--- /issm/trunk-jpl/src/c/Container/Observations.cpp	(revision 13260)
+++ /issm/trunk-jpl/src/c/Container/Observations.cpp	(revision 13261)
@@ -36,12 +36,12 @@
 }
 /*}}}*/
-/*FUNCTION Observations::Observations(IssmDouble* observations_list,IssmDouble* x,IssmDouble* y,int n,Options* options){{{*/
-Observations::Observations(IssmDouble* observations_list,IssmDouble* x,IssmDouble* y,int n,Options* options){
+/*FUNCTION Observations::Observations(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options){{{*/
+Observations::Observations(IssmPDouble* observations_list,IssmPDouble* x,IssmPDouble* y,int n,Options* options){
 
 	/*Intermediaries*/
 	int          i,j,maxdepth,level,counter,index;
 	int          xi,yi;
-	IssmDouble       xmin,xmax,ymin,ymax;
-	IssmDouble       offset,minlength,minspacing,mintrimming,maxtrimming;
+	IssmPDouble       xmin,xmax,ymin,ymax;
+	IssmPDouble       offset,minlength,minspacing,mintrimming,maxtrimming;
 	Observation *observation = NULL;
 
@@ -66,9 +66,9 @@
 		options->Get(&minlength,"boxlength");
 		if(minlength<=0)_error_("boxlength should be a positive number");
-		maxdepth=reCast<int,IssmDouble>(log(max(xmax-xmin,ymax-ymin)/minlength +1)/log(2.0));
+		maxdepth=reCast<int,IssmPDouble>(log(max(xmax-xmin,ymax-ymin)/minlength +1)/log(2.0));
 	}
 	else{
 		maxdepth = 30;
-		minlength=max(xmax-xmin,ymax-ymin)/IssmDouble((1L<<maxdepth)-1);
+		minlength=max(xmax-xmin,ymax-ymin)/IssmPDouble((1L<<maxdepth)-1);
 	}
 
@@ -118,10 +118,10 @@
 /*Methods*/
 /*FUNCTION Observations::ClosestObservation{{{*/
-void Observations::ClosestObservation(IssmDouble *px,IssmDouble *py,IssmDouble *pobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius){
+void Observations::ClosestObservation(IssmPDouble *px,IssmPDouble *py,IssmPDouble *pobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){
 
 	/*Output and Intermediaries*/
 	bool         stop;
 	int          nobs,i,index;
-	IssmDouble       h2,hmin2,radius2;
+	IssmPDouble  h2,hmin2,radius2;
 	int         *indices      = NULL;
 	Observation *observation  = NULL;
@@ -166,17 +166,17 @@
 
 }/*}}}*/
-/*FUNCTION Observations::ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata){{{*/
-void Observations::ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int maxdata){
+/*FUNCTION Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata){{{*/
+void Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int maxdata){
 
 	/*Output and Intermediaries*/
 	bool         stop;
 	int          nobs,tempnobs,i,j,k,n,counter;
-	IssmDouble       h2,radius2;
+	IssmPDouble  h2,radius2;
 	int         *indices      = NULL;
 	int         *tempindices  = NULL;
-	IssmDouble      *dists        = NULL;
-	IssmDouble      *x            = NULL;
-	IssmDouble      *y            = NULL;
-	IssmDouble      *obs          = NULL;
+	IssmPDouble *dists        = NULL;
+	IssmPDouble *x            = NULL;
+	IssmPDouble *y            = NULL;
+	IssmPDouble *obs          = NULL;
 	Observation *observation  = NULL;
 
@@ -191,5 +191,5 @@
 	if(tempnobs){
 		indices = xNew<int>(tempnobs);
-		dists   = xNew<IssmDouble>(tempnobs);
+		dists   = xNew<IssmPDouble>(tempnobs);
 	}
 	nobs = 0;
@@ -226,12 +226,12 @@
 		}
 	}  
-	xDelete<IssmDouble>(dists);
+	xDelete<IssmPDouble>(dists);
 	xDelete<int>(tempindices);
 
 	if(nobs){
 		/*Allocate vectors*/
-		x   = xNew<IssmDouble>(nobs);
-		y   = xNew<IssmDouble>(nobs);
-		obs = xNew<IssmDouble>(nobs);
+		x   = xNew<IssmPDouble>(nobs);
+		y   = xNew<IssmPDouble>(nobs);
+		obs = xNew<IssmPDouble>(nobs);
 
 		/*Loop over all observations and fill in x, y and obs*/
@@ -249,12 +249,12 @@
 	*pnobs=nobs;
 }/*}}}*/
-/*FUNCTION Observations::ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs){{{*/
-void Observations::ObservationList(IssmDouble **px,IssmDouble **py,IssmDouble **pobs,int* pnobs){
+/*FUNCTION Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs){{{*/
+void Observations::ObservationList(IssmPDouble **px,IssmPDouble **py,IssmPDouble **pobs,int* pnobs){
 
 	/*Output and Intermediaries*/
 	int          nobs;
-	IssmDouble      *x            = NULL;
-	IssmDouble      *y            = NULL;
-	IssmDouble      *obs          = NULL;
+	IssmPDouble *x            = NULL;
+	IssmPDouble *y            = NULL;
+	IssmPDouble *obs          = NULL;
 	Observation *observation  = NULL;
 
@@ -262,7 +262,7 @@
 
 	if(nobs){
-		x   = xNew<IssmDouble>(nobs);
-		y   = xNew<IssmDouble>(nobs);
-		obs = xNew<IssmDouble>(nobs);
+		x   = xNew<IssmPDouble>(nobs);
+		y   = xNew<IssmPDouble>(nobs);
+		obs = xNew<IssmPDouble>(nobs);
 		for(int i=0;i<this->Size();i++){
 			observation=(Observation*)this->GetObjectByOffset(i);
@@ -278,13 +278,13 @@
 }/*}}}*/
 /*FUNCTION Observations::InterpolationIDW{{{*/
-void Observations::InterpolationIDW(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,IssmDouble power){
+void Observations::InterpolationIDW(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,IssmPDouble power){
 
 	/*Intermediaries*/
-	int    i,n_obs;
-	IssmDouble prediction;
-	IssmDouble numerator,denominator,h,weight;
-	IssmDouble *x   = NULL;
-	IssmDouble *y   = NULL;
-	IssmDouble *obs = NULL;
+	int         i,n_obs;
+	IssmPDouble prediction;
+	IssmPDouble numerator,denominator,h,weight;
+	IssmPDouble *x   = NULL;
+	IssmPDouble *y   = NULL;
+	IssmPDouble *obs = NULL;
 
 	/*Some checks*/
@@ -322,24 +322,24 @@
 	/*clean-up*/
 	*pprediction = prediction;
-	xDelete<IssmDouble>(x);
-	xDelete<IssmDouble>(y);
-	xDelete<IssmDouble>(obs);
+	xDelete<IssmPDouble>(x);
+	xDelete<IssmPDouble>(y);
+	xDelete<IssmPDouble>(obs);
 }/*}}}*/
 /*FUNCTION Observations::InterpolationKriging{{{*/
-void Observations::InterpolationKriging(IssmDouble *pprediction,IssmDouble *perror,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius,int mindata,int maxdata,Variogram* variogram){
+void Observations::InterpolationKriging(IssmPDouble *pprediction,IssmPDouble *perror,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius,int mindata,int maxdata,Variogram* variogram){
 
 	/*Intermediaries*/
 	int           i,j,n_obs;
-	IssmDouble        prediction,error;
-	IssmDouble        numerator,denominator,ratio;
-	IssmDouble       *x            = NULL;
-	IssmDouble       *y            = NULL;
-	IssmDouble       *obs          = NULL;
-	IssmDouble       *Gamma        = NULL;
-	IssmDouble       *GinvG0       = NULL;
-	IssmDouble       *Ginv1        = NULL;
-	IssmDouble       *GinvZ        = NULL;
-	IssmDouble       *gamma0       = NULL;
-	IssmDouble       *ones         = NULL;
+	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;
 
 	/*Some checks*/
@@ -361,7 +361,7 @@
 
 	/*Allocate intermediary matrix and vectors*/
-	Gamma  = xNew<IssmDouble>(n_obs*n_obs);
-	gamma0 = xNew<IssmDouble>(n_obs);
-	ones   = xNew<IssmDouble>(n_obs);
+	Gamma  = xNew<IssmPDouble>(n_obs*n_obs);
+	gamma0 = xNew<IssmPDouble>(n_obs);
+	ones   = xNew<IssmPDouble>(n_obs);
 
 	/*First: Create semivariogram matrix for observations*/
@@ -402,20 +402,20 @@
 	*pprediction = prediction;
 	*perror = error;
-	xDelete<IssmDouble>(x);
-	xDelete<IssmDouble>(y);
-	xDelete<IssmDouble>(obs);
-	xDelete<IssmDouble>(Gamma);
-	xDelete<IssmDouble>(gamma0);
-	xDelete<IssmDouble>(ones);
-	xDelete<IssmDouble>(GinvG0);
-	xDelete<IssmDouble>(Ginv1);
-	xDelete<IssmDouble>(GinvZ);
+	xDelete<IssmPDouble>(x);
+	xDelete<IssmPDouble>(y);
+	xDelete<IssmPDouble>(obs);
+	xDelete<IssmPDouble>(Gamma);
+	xDelete<IssmPDouble>(gamma0);
+	xDelete<IssmPDouble>(ones);
+	xDelete<IssmPDouble>(GinvG0);
+	xDelete<IssmPDouble>(Ginv1);
+	xDelete<IssmPDouble>(GinvZ);
 
 }/*}}}*/
 /*FUNCTION Observations::InterpolationNearestNeighbor{{{*/
-void Observations::InterpolationNearestNeighbor(IssmDouble *pprediction,IssmDouble x_interp,IssmDouble y_interp,IssmDouble radius){
+void Observations::InterpolationNearestNeighbor(IssmPDouble *pprediction,IssmPDouble x_interp,IssmPDouble y_interp,IssmPDouble radius){
 
 	/*Intermediaries*/
-	IssmDouble        x,y,obs;
+	IssmPDouble x,y,obs;
 
 	/*Get clostest observation*/
@@ -426,5 +426,5 @@
 }/*}}}*/
 /*FUNCTION Observations::QuadtreeColoring{{{*/
-void Observations::QuadtreeColoring(IssmDouble* A,IssmDouble *x,IssmDouble *y,int n){
+void Observations::QuadtreeColoring(IssmPDouble* A,IssmPDouble *x,IssmPDouble *y,int n){
 
 	int xi,yi,level;
@@ -433,18 +433,18 @@
 		this->quadtree->IntergerCoordinates(&xi,&yi,x[i],y[i]);
 		this->quadtree->QuadtreeDepth(&level,xi,yi);
-		A[i]=(IssmDouble)level;
+		A[i]=(IssmPDouble)level;
 	}
 
 }/*}}}*/
 /*FUNCTION Observations::Variomap{{{*/
-void Observations::Variomap(IssmDouble* gamma,IssmDouble *x,int n){
+void Observations::Variomap(IssmPDouble* gamma,IssmPDouble *x,int n){
 
 	/*Output and Intermediaries*/
 	int          i,j,k;
-	IssmDouble       distance;
+	IssmPDouble  distance;
 	Observation *observation1 = NULL;
 	Observation *observation2 = NULL;
 
-	IssmDouble *counter = xNew<IssmDouble>(n);
+	IssmPDouble *counter = xNew<IssmPDouble>(n);
 	for(j=0;j<n;j++) counter[j] = 0.0;
 	for(j=0;j<n;j++) gamma[j]   = 0.0;
@@ -475,4 +475,4 @@
 
 	/*Assign output pointer*/
-	xDelete<IssmDouble>(counter);
-}/*}}}*/
+	xDelete<IssmPDouble>(counter);
+}/*}}}*/
Index: /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h	(revision 13260)
+++ /issm/trunk-jpl/src/c/modules/Solverx/Solverx.h	(revision 13261)
@@ -24,4 +24,5 @@
 
 void SolverxSeq(SeqVec<IssmDouble>** puf,SeqMat<IssmDouble>* Kff, SeqVec<IssmDouble>* pf,Parameters* parameters);
+void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n);
 void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n);
 
Index: /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp	(revision 13260)
+++ /issm/trunk-jpl/src/c/modules/Solverx/SolverxSeq.cpp	(revision 13261)
@@ -52,19 +52,30 @@
 
 }/*}}}*/
+void SolverxSeq(IssmPDouble **pX, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
+
+	/*Allocate output*/
+	double* X = xNew<double>(n); 
+
+	/*Solve*/
+	SolverxSeq(X,A,B,n);
+
+	/*Assign output pointer*/
+	*pX=X;
+}
+/*}}}*/
+
 #ifdef _HAVE_ADOLC_
-
-int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y) {
+int EDF_for_solverx(int n, IssmPDouble *x, int m, IssmPDouble *y){/*{{{*/
     if(m*(m+1)!=n)_error_("Stiffness matrix should be square!");
     SolverxSeq(y,x, x+m*m, m);
     return 0;
-}
-
-void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){//{{{
+}/*}}}*/
+void SolverxSeq(IssmDouble *X,IssmDouble *A,IssmDouble *B,int n, Parameters* parameters){/*{{{*/
 	// pack inputs to conform to the EDF-prescribed interface
-        IssmDouble*  adoubleEDFin=xNew<IssmDouble>(n*(n+1)); // packed inputs, i.e. matrix and right hand side
-        for(int i=0; i<n*n;i++)adoubleEDFin[i]    =A[i]; // pack matrix
-        for(int i=0; i<n;  i++)adoubleEDFin[i+n*n]=B[i]; // pack the right hand side
+        IssmDouble*  adoubleEDFin=xNew<IssmDouble>(n*(n+1));  // packed inputs, i.e. matrix and right hand side
+        for(int i=0; i<n*n;i++)adoubleEDFin[i]    =A[i];      // pack matrix
+        for(int i=0; i<n;  i++)adoubleEDFin[i+n*n]=B[i];      // pack the right hand side
         IssmPDouble* pdoubleEDFin=xNew<IssmPDouble>(n*(n+1)); // provide space to transfer inputs during call_ext_fct
-	IssmPDouble* pdoubleEDFout=xNew<IssmPDouble>(n);       // provide space to transfer outputs during call_ext_fct
+	IssmPDouble* pdoubleEDFout=xNew<IssmPDouble>(n);           // provide space to transfer outputs during call_ext_fct
 	// call the wrapped solver through the registry entry we retrieve from parameters
 	call_ext_fct(dynamic_cast<GenericParam<Adolc_edf> * >(parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p,
@@ -77,6 +88,6 @@
 /*}}}*/
 #endif
-void SolverxSeq(IssmPDouble * X, IssmPDouble * A, IssmPDouble * B,int n){ //{{{
-	#ifdef _HAVE_GSL_
+void SolverxSeq(IssmPDouble *X, IssmPDouble *A, IssmPDouble *B,int n){ /*{{{*/
+#ifdef _HAVE_GSL_
 	/*GSL Matrices and vectors: */
 	int              s;
@@ -109,5 +120,5 @@
 	gsl_permutation_free(p);
 	gsl_vector_free(x);
-	#endif
+#endif
 }
 /*}}}*/