Changeset 18457

Timestamp:

08/21/14 16:06:57 (11 years ago)

Author:

abuzzi

Message:

MODIF:final version of atsu cpy multithread

File:

• issm/trunk-jpl/src/m/contrib/gravity/vfsa.cpp (modified) (8 diffs)

issm/trunk-jpl/src/m/contrib/gravity/vfsa.cpp

@@ -8 +8 @@
 #include <cassert>
 #include <gsl/gsl_multifit.h>
+#include <pthread.h>
 
 using namespace std;
@@ -126 +127 @@
 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 +136 @@
 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 +415 @@
                         }
                 }
+        cout<<totaliter<<endl;
         }
 
@@ -400 +425 @@
                 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 +509 @@
         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 +537 @@
         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 +667 @@
         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;
+}/*}}}*/