Changeset 2549

Timestamp:

10/27/09 15:28:33 (15 years ago)

Author:

Eric.Larour

Message:

Added multi-threading to InterpFromGridToMeshx.

Location:

issm/trunk/src/c

Files:

• InterpFromGridToMeshx/InterpFromGridToMeshx.cpp (modified) (4 diffs)
• InterpFromGridToMeshx/InterpFromGridToMeshx.h (modified) (1 diff)
• InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp (added)
• InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dxt.cpp (modified) (2 diffs)
• Makefile.am (modified) (4 diffs)
• shared/Threads/PartitionRange.cpp (added)
• shared/Threads/issm_threads.h (modified) (1 diff)

issm/trunk/src/c/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp

@@ -2 +2 @@
  * \brief  "c" core code for interpolating values from a structured grid.
  */ 
+
+#ifdef HAVE_CONFIG_H
+        #include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
 
 #include "./InterpFromGridToMeshx.h"
@@ -9 +16 @@
 #define __FUNCT__ "InterpFromGridToMeshx"
 
-int findindices(int* pm,int* pn,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid);
 
 int InterpFromGridToMeshx( Vec* pdata_mesh,double* x_in, int x_rows, double* y_in, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods,double default_value) {
@@ -20 +26 @@
         double* x=NULL;
         double* y=NULL;
-        int i,m,n;
-        double x_grid,y_grid;
-        double xi,eta;
-        double G1,G2,G3,G4,data_value;
-        double area_1,area_2,area_3;
-        double area;
+        double  x_grid,y_grid;
+        int     i;
+
+        /*threading: */
+        InterpFromGridToMeshxThreadStruct gate;
+        int num=1;
+
+        #ifdef _MULTITHREADING_
+        num=_NUMTHREADS_;
+        #endif
 
         /*Some checks on arguments: */
@@ -58 +68 @@
         }
 
-        /*Linear (triangle) interpolation: */
-        for ( i=MPI_Lowerrow(nods); i<MPI_Upperrow(nods); i++) {
+        /*initialize thread parameters: */
+        gate.x_mesh=x_mesh;
+        gate.y_mesh=y_mesh;
+        gate.x_rows=x_rows;
+        gate.y_rows=y_rows;
+        gate.x=x;
+        gate.y=y;
+        gate.nods=nods;
+        gate.data_mesh=data_mesh;
+        gate.data=data;
+        gate.M=M;
+        gate.N=N;
 
-                x_grid=*(x_mesh+i);
-                y_grid=*(y_mesh+i);
-
-                /*Find indices m and n into y and x, for which  y(m)<=y_grids<=y(m+1) and x(n)<=x_grid<=x(n+1)*/
-                if(findindices(&n,&m,x,x_rows, y,y_rows, x_grid,y_grid)){
-                        
-                        /*Get area*/
-                        area=(x[n+1]-x[n])*(y[m+1]-y[m]);
-
-                        /*is it the upper right triangle?*/
-                        /*2'     3'
-                         *+-----+
-                         *1\    |
-                         *| \   |
-                         *|  \  |
-                         *|   \ |
-                         *|    \|
-                         *2----3+1' */
-                        if ((x_grid-x[n])/(x[n+1]-x[n])<(y_grid-y[m])/(y[m+1]-y[m])){
-
-                                /*Then find values of data at each summit*/
-                                G1=*(data+m*N+n);
-                                G2=*(data+(m+1)*N+n+1);
-                                G3=*(data+(m+1)*N+n);
-
-                                /*Get first area coordinate*/
-                                area_1=((y[m+1]-y_grid)*(x[n+1]-x[n]))/area;
-                                /*Get second area coordinate*/
-                                area_2=((x_grid-x[n])*(y[m+1]-y[m]))/area;
-                                /*Get third area coordinate = 1-area1-area2*/
-                                area_3=1-area_1-area_2;
-
-                                /*interpolate*/
-                                data_value=area_1*G1+area_2*G2+area_3*G3;
-                        }
-                        else {
-
-                                /*Then find values of data at each summit*/
-                                G1=*(data+(m+1)*N+n+1);
-                                G2=*(data+m*N+n);
-                                G3=*(data+m*N+n+1);
-
-                                /*Get first area coordinate*/
-                                area_1=((y_grid-y[m])*(x[n+1]-x[n]))/area;
-                                /*Get second area coordinate*/
-                                area_2=((x[n+1]-x_grid)*(y[m+1]-y[m]))/area;
-                                /*Get third area coordinate = 1-area1-area2*/
-                                area_3=1-area_1-area_2;
-
-                                /*interpolate*/
-                                data_value=area_1*G1+area_2*G2+area_3*G3;
-                        }
-
-                        /*Treat NANs*/
-                        if (isnan(data_value)){
-                                data_value=default_value;
-                        }
-                }
-                else{
-                        data_value=default_value;
-                }
-                VecSetValue(data_mesh,i,data_value,INSERT_VALUES);
-        }
+        /*launch the thread manager with InterpFromGridToMeshxt as a core: */
+        LaunchThread(InterpFromGridToMeshxt,(void*)&gate,num);
 
         /*Assign output pointers:*/
         *pdata_mesh=data_mesh;
 }
-
-int findindices(int* pm,int* pn,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid){
-
-        int foundx=0;
-        int foundy=0;
-        int i;
-        int m=-1;
-        int n=-1;
-
-        for (i=0;i<x_rows-1;i++){
-                if ( (*(x+i)<=xgrid) && (xgrid<*(x+i+1)) ){
-                        m=i;
-                        foundx= 1;
-                        break;
-                }
-        }
-        if(*(x+x_rows-1)==xgrid){
-                m=x_rows-2;
-                foundx=1;
-        }
-        
-        for (i=0;i<y_rows-1;i++){
-                if ( (*(y+i)<=ygrid) && (ygrid<*(y+i+1)) ){
-                        n=i;
-                        foundy= 1;
-                        break;
-                }
-        }
-        if(*(y+y_rows-1)==ygrid){
-                m=y_rows-2;
-                foundy=1;
-        }
-
-        /*Assign output pointers:*/
-        *pm=m;
-        *pn=n;
-        return foundx*foundy;
-}

issm/trunk/src/c/InterpFromGridToMeshx/InterpFromGridToMeshx.h

@@ -8 +8 @@
 #include "../toolkits/toolkits.h"
 
+
+
+/*threading: */
+typedef struct{
+
+        double* x_mesh;
+        double* y_mesh;
+        int     x_rows;
+        int     y_rows;
+        double* x;
+        double* y;
+        int     nods;
+        double  default_value;
+        Vec     data_mesh;
+        double* data;
+        int     M;
+        int     N;
+
+} InterpFromGridToMeshxThreadStruct;
+
+int findindices(int* pm,int* pn,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid);
+
 int InterpFromGridToMeshx( Vec* pdata_mesh,double* x, int x_rows, double* y, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods, double default_value);
+
+void* InterpFromGridToMeshxt(void* vInterpFromGridToMeshxThreadStruct);
+
 
 #endif /* _INTERPFROMGRIDTOMESHX_H */

issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dxt.cpp

@@ -45 +45 @@
         double area_1,area_2,area_3;
         double data_value;
-        int    step;
 
         /*recover handle and gate: */
@@ -72 +71 @@
 
         /*distribute elements across threads :*/
-        step=(int)floor((double)nels_data/(double)num_threads);
-        for(i=0;i<(my_thread+1);i++){
-                i0=i*step;
-                if(i==(num_threads-1))i1=nels_data;
-                else i1=i0+step;
-        }
+        PartitionRange(&i0,&i1,nels_data,num_threads,my_thread);
         
         for (i=i0;i<i1;i++){

issm/trunk/src/c/Makefile.am

@@ -89 +89 @@
                                         ./shared/Threads/issm_threads.h\
                                         ./shared/Threads/LaunchThread.cpp\
+                                        ./shared/Threads/PartitionRange.cpp\
                                         ./shared/Matlab/matlabshared.h\
                                         ./shared/Matlab/PrintfFunction.cpp\
@@ -249 +250 @@
                                         ./InterpFromGridToMeshx/InterpFromGridToMeshx.cpp\
                                         ./InterpFromGridToMeshx/InterpFromGridToMeshx.h\
+                                        ./InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp\
                                         ./InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp\
                                         ./InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dxt.cpp\
@@ -386 +388 @@
                                         ./shared/Threads/issm_threads.h\
                                         ./shared/Threads/LaunchThread.cpp\
+                                        ./shared/Threads/PartitionRange.cpp\
                                         ./shared/Alloc/alloc.h\
                                         ./shared/Alloc/alloc.cpp\
@@ -545 +548 @@
                                         ./InterpFromGridToMeshx/InterpFromGridToMeshx.cpp\
                                         ./InterpFromGridToMeshx/InterpFromGridToMeshx.h\
+                                        ./InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp\
                                         ./InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp\
                                         ./InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dxt.cpp\

issm/trunk/src/c/shared/Threads/issm_threads.h

@@ -18 +18 @@
  * or just serially if not requested: */
 void LaunchThread(void* function(void*), void* gate,int num_threads);
+void PartitionRange(int* pi0,int* pi1, int num_el,int num_threads,int my_thread);
 
 #endif //ifndef _ISSM_THREADS_H_