Index: /issm/trunk-jpl/src/m/contrib/gravity/vfsa.cpp
===================================================================
--- /issm/trunk-jpl/src/m/contrib/gravity/vfsa.cpp	(revision 18456)
+++ /issm/trunk-jpl/src/m/contrib/gravity/vfsa.cpp	(revision 18457)
@@ -8,4 +8,5 @@
 #include <cassert>
 #include <gsl/gsl_multifit.h>
+#include <pthread.h>
 
 using namespace std;
@@ -126,5 +127,5 @@
 void vec2grid(Matrix *V,Matrix *V1,Matrix *V2,int nx, int ny);
 void msplit( Matrix *m, Matrix *m1,Matrix *m2,double dlevel);
-void plouff(Matrix *g,Matrix *Pobs,Matrix *Pp,Matrix * mesh,Matrix *rho,int dx,int dy, int dn,int m,int n,int l);
+void* plouffT(void* vpthread_handle);
 void vec2gridsimple(Matrix *V,Matrix *V1,int nx, int ny);
 void reshape(Matrix* V,Matrix* V1,int nx,int ny);
@@ -135,4 +136,27 @@
 void newmodelgen(Matrix* m0,Matrix* m1,Matrix* bathy,Matrix* icethick,int mx,int my,double T,double ptval,double mmax,double mmax2,double ctr,double sd);
 double coolshed(double T0,int k,double c,double D);
+void   LaunchThread(void* function(void*), void* usr,int num_threads);
+/*}}}*/
+
+/*Multithreading structures {{{*/
+typedef struct{
+	void* usr;
+	int   my_thread;
+	int   num_threads;
+} pthread_handle;
+
+typedef struct{
+	Matrix *g;
+	Matrix *Pobs;
+	Matrix *Pp;
+	Matrix *mesh;
+	Matrix *rho;
+	int dx;
+	int dy;
+	int dn;
+	int m;
+	int n;
+	int l;
+} AppStruct;
 /*}}}*/
 
@@ -391,4 +415,5 @@
 			}
 		}
+	cout<<totaliter<<endl;
 	}
 
@@ -400,6 +425,6 @@
 		savefile2 << m_min->GetValue(i,0)<<"   "<< m_min->GetValue(i,1)<<"   "<< m_min->GetValue(i,2)<<endl;
 	}
-	savefile1.close;
-	savefile2.close;
+	savefile1.close();
+	savefile2.close();
 			/*}}}*/
 	return 0;
@@ -484,77 +509,4 @@
 	m2->AddNumber(dlevel);
 }/*}}}*/
-void plouff(Matrix *g,Matrix *Pobs,Matrix *Pp,Matrix * mesh,Matrix *rho,int dx,int dy, int dn,int m,int n,int l){/*{{{*/
-		double gg=6.673e-11;
-		int si,sj,id,s;
-		double R,Q,P;
-		Matrix *xp= new Matrix(1,2);
-		Matrix *yp= new Matrix(1,2);
-		Matrix *xpp= new Matrix(1,2);
-		Matrix *ypp= new Matrix(1,2);
-		Matrix *U= new Matrix(l,4);
-		Matrix *U1=new Matrix(1,4);
-		Matrix *U2=new Matrix(1,4);
-		Matrix *gl= new Matrix(1,l-1);
-		bool test=true;
-		for(int c=0;c<n;c++){
-			g->SetValue(c,0,0);
-			for(int a=0;a<m;a++){
-				test=true;
-				xp->SetValue(0,0,Pp->GetValue(a,0)-Pobs->GetValue(c,0));
-				xp->SetValue(0,1,Pp->GetValue(a,0)-Pobs->GetValue(c,0)+dx);
-				if(xp->GetValue(0,0)<0 && xp->GetValue(0,0)<xp->GetValue(0,1) && xp->GetValue(0,0)*xp->GetValue(0,1)>=0){
-					xpp->SetValue(0,0,xp->GetValue(0,1));
-					xpp->SetValue(0,1,xp->GetValue(0,0));
-					xp->MatrixAbs(xpp);
-				}
-				yp->SetValue(0,0,Pp->GetValue(a,1)-Pobs->GetValue(c,1));
-				yp->SetValue(0,1,Pp->GetValue(a,1)-Pobs->GetValue(c,1)+dy);
-				if(yp->GetValue(0,0)<0 && yp->GetValue(0,0)<yp->GetValue(0,1) && yp->GetValue(0,0)*yp->GetValue(0,1)>=0){
-					ypp->SetValue(0,0,yp->GetValue(0,1));
-					ypp->SetValue(0,1,yp->GetValue(0,0));
-					yp->MatrixAbs(ypp);
-				}
-				P=sqrt(xp->GetValue(0,0)*xp->GetValue(0,0)+yp->GetValue(0,0)*yp->GetValue(0,0));
-				if(P>dn){
-					test=false;
-					for(int i=0;i<l-1;i++){
-						gl->SetValue(0,i,0);
-					}
-				}
-				if(test==true){
-					si=1;
-					sj=1;
-					id=0;
-					for(int i=0;i<2;i++){
-						si*=-1;
-						for(int j=0;j<2;j++){
-							si*=-1;
-							s=si*sj;
-							for(int k=0;k<l;k++){
-								R=sqrt(xp->GetValue(0,i)*xp->GetValue(0,i)+yp->GetValue(0,j)*yp->GetValue(0,j)+mesh->GetValue(a,k)*mesh->GetValue(a,k));
-								Q=atan(xp->GetValue(0,i)*yp->GetValue(0,j)/(mesh->GetValue(a,k)*R));
-								U->SetValue(k,id,s*(mesh->GetValue(a,k)*Q-xp->GetValue(0,i)*log(R+yp->GetValue(0,j))-yp->GetValue(0,j)*log(R+xp->GetValue(0,i))));
-							}
-							id++;
-						}
-					}
-					for(int b=0;b<l-1;b++){
-						U->ExtractLine(U1,b);
-						U->ExtractLine(U2,b+1);
-						gl->SetValue(0,b,rho->GetValue(0,b)*(U1->MatrixInternSum()*(-1)+U2->MatrixInternSum()));
-					}
-				}
-				g->SetValue(c,0,g->GetValue(c,0)+gg*gl->MatrixInternSum());
-			}
-		}
-		delete xp;
-		delete yp;
-		delete xpp;
-		delete ypp;
-		delete gl;
-		delete U;
-		delete U1;
-		delete U2;
-	}/*}}}*/
 void vec2gridsimple(Matrix *V,Matrix *V1,int nx, int ny){/*{{{*/
 	for(int i=0;i<ny;i++){
@@ -585,6 +537,25 @@
 	double e=0;
 	msplit(m0,m1,m2,dlevel);
-   plouff(g1,Pobs,Pp,m1,rho1,dx,dy,dn,mx*my,nx*ny,4);
-   plouff(g2,Pobs,Pp,m2,rho2,dx,dy,dn,mx*my,nx*ny,3);
+	
+	/*Multithreaded core*/
+	int       num_threads = 4;
+	AppStruct usr;
+	usr.g = g1;
+	usr.Pobs = Pobs;
+	usr.Pp=Pp;
+	usr.mesh= m1;
+	usr.rho= rho1;
+	usr.dx=dx;
+	usr.dy=dy;
+	usr.dn=dn;
+	usr.m=mx*my;
+	usr.n=nx*ny;
+	usr.l=4;
+	LaunchThread(plouffT,(void*)&usr,num_threads);
+	usr.g = g2;
+	usr.mesh= m2;
+	usr.rho= rho2;
+	usr.l=3;
+	LaunchThread(plouffT,(void*)&usr,num_threads);
 	g->MatrixSum(g1,g2);
 	vec2gridsimple(g,gcalgr,nx,ny);
@@ -696,2 +667,128 @@
 	return T1;
 }/*}}}*/
+void* plouffT(void* vpthread_handle){/*{{{*/
+
+	/*recover this thread info*/
+	pthread_handle *handle = (pthread_handle*)vpthread_handle;
+	int my_thread   = handle->my_thread;
+	int num_threads = handle->num_threads;
+
+	/*Recover struct*/
+	AppStruct *usr = (AppStruct*)handle->usr;
+	Matrix *g = usr->g;
+	Matrix *Pobs = usr->Pobs;
+	Matrix *Pp = usr->Pp;
+	Matrix *mesh = usr->mesh;
+	Matrix *rho = usr->rho;
+	int dx =usr->dx;
+	int dy =usr->dy;
+	int dn =usr->dn;
+	int m =usr->m;
+	int n =usr->n;
+	int l =usr->l;
+
+	double gg=6.673e-11;
+	int si,sj,id,s;
+	double R,Q,P;
+	Matrix *xp= new Matrix(1,2);
+	Matrix *yp= new Matrix(1,2);
+	Matrix *xpp= new Matrix(1,2);
+	Matrix *ypp= new Matrix(1,2);
+	Matrix *U= new Matrix(l,4);
+	Matrix *U1=new Matrix(1,4);
+	Matrix *U2=new Matrix(1,4);
+	Matrix *gl= new Matrix(1,l-1);
+	bool test=true;
+	for(int c=my_thread;c<n;c+=num_threads){
+		g->SetValue(c,0,0);
+		for(int a=0;a<m;a++){
+			test=true;
+			xp->SetValue(0,0,Pp->GetValue(a,0)-Pobs->GetValue(c,0));
+			xp->SetValue(0,1,Pp->GetValue(a,0)-Pobs->GetValue(c,0)+dx);
+			if(xp->GetValue(0,0)<0 && xp->GetValue(0,0)<xp->GetValue(0,1) && xp->GetValue(0,0)*xp->GetValue(0,1)>=0){
+				xpp->SetValue(0,0,xp->GetValue(0,1));
+				xpp->SetValue(0,1,xp->GetValue(0,0));
+				xp->MatrixAbs(xpp);
+			}
+			yp->SetValue(0,0,Pp->GetValue(a,1)-Pobs->GetValue(c,1));
+			yp->SetValue(0,1,Pp->GetValue(a,1)-Pobs->GetValue(c,1)+dy);
+			if(yp->GetValue(0,0)<0 && yp->GetValue(0,0)<yp->GetValue(0,1) && yp->GetValue(0,0)*yp->GetValue(0,1)>=0){
+				ypp->SetValue(0,0,yp->GetValue(0,1));
+				ypp->SetValue(0,1,yp->GetValue(0,0));
+				yp->MatrixAbs(ypp);
+			}
+			P=sqrt(xp->GetValue(0,0)*xp->GetValue(0,0)+yp->GetValue(0,0)*yp->GetValue(0,0));
+			if(P>dn){
+				test=false;
+				for(int i=0;i<l-1;i++){
+					gl->SetValue(0,i,0);
+				}
+			}
+			if(test==true){
+				si=1;
+				sj=1;
+				id=0;
+				for(int i=0;i<2;i++){
+					si*=-1;
+					for(int j=0;j<2;j++){
+						si*=-1;
+						s=si*sj;
+						for(int k=0;k<l;k++){
+							R=sqrt(xp->GetValue(0,i)*xp->GetValue(0,i)+yp->GetValue(0,j)*yp->GetValue(0,j)+mesh->GetValue(a,k)*mesh->GetValue(a,k));
+							Q=atan(xp->GetValue(0,i)*yp->GetValue(0,j)/(mesh->GetValue(a,k)*R));
+							U->SetValue(k,id,s*(mesh->GetValue(a,k)*Q-xp->GetValue(0,i)*log(R+yp->GetValue(0,j))-yp->GetValue(0,j)*log(R+xp->GetValue(0,i))));
+						}
+						id++;
+					}
+				}
+				for(int b=0;b<l-1;b++){
+					U->ExtractLine(U1,b);
+					U->ExtractLine(U2,b+1);
+					gl->SetValue(0,b,rho->GetValue(0,b)*(U1->MatrixInternSum()*(-1)+U2->MatrixInternSum()));
+				}
+			}
+			g->SetValue(c,0,g->GetValue(c,0)+gg*gl->MatrixInternSum());
+		}
+	}
+	delete xp;
+	delete yp;
+	delete xpp;
+	delete ypp;
+	delete gl;
+	delete U;
+	delete U1;
+	delete U2;
+
+	return NULL;
+}/*}}}*/
+void LaunchThread(void* function(void*), void* usr,int num_threads){/*{{{*/
+
+	int i;
+	int            *status  = NULL;
+	pthread_t      *threads = NULL;
+	pthread_handle *handles = NULL;
+
+	/*dynamically allocate: */
+	threads=(pthread_t*)malloc(num_threads*sizeof(pthread_t));
+	handles=(pthread_handle*)malloc(num_threads*sizeof(pthread_handle));
+
+	for(i=0;i<num_threads;i++){
+		handles[i].usr=usr;
+		handles[i].my_thread  =i;
+		handles[i].num_threads=num_threads;
+	}
+	for(i=0;i<num_threads;i++){
+		if(pthread_create(threads+i,NULL,function,(void*)(handles+i))){
+			std::cerr<<"pthread_create error";
+		}
+	}
+	for(i=0;i<num_threads;i++){
+		if(pthread_join(threads[i],(void**)&status)){
+			std::cerr<<"pthread_join error";
+		}
+	}
+
+	/*Free ressources:*/
+	delete threads;
+	delete handles;
+}/*}}}*/