/[escript]/trunk/paso/src/PasoUtil.c
ViewVC logotype

Diff of /trunk/paso/src/PasoUtil.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3317 by gross, Tue Oct 19 00:28:00 2010 UTC revision 3318 by caltinay, Thu Oct 28 01:05:36 2010 UTC
# Line 59  bool_t Paso_Util_isAny(dim_t N,index_t* Line 59  bool_t Paso_Util_isAny(dim_t N,index_t*
59  index_t Paso_Util_cumsum(dim_t N,index_t* array) {  index_t Paso_Util_cumsum(dim_t N,index_t* array) {
60     index_t out=0,tmp;     index_t out=0,tmp;
61     dim_t i;     dim_t i;
62  if (omp_get_max_threads()>1) {  #ifdef _OPENMP
63       const int num_threads=omp_get_max_threads();
64       const int thread_num=omp_get_thread_num();
65    #else
66       const int num_threads=1;
67       const int thread_num=0;
68    #endif
69       if (num_threads>1) {
70        index_t *partial_sums=NULL,sum;        index_t *partial_sums=NULL,sum;
71        partial_sums=TMPMEMALLOC(omp_get_max_threads(),index_t);        partial_sums=TMPMEMALLOC(num_threads, index_t);
72        #pragma omp parallel private(sum,i)        #pragma omp parallel private(sum,i)
73        {        {
74          sum=0;          sum=0;
75          #pragma omp for schedule(static)          #pragma omp for schedule(static)
76          for (i=0;i<N;++i) sum+=array[i];          for (i=0;i<N;++i) sum+=array[i];
77    
78          partial_sums[omp_get_thread_num()]=sum;          partial_sums[thread_num]=sum;
79        }        }
80    
81        {        {
82            out=0;            out=0;
83            for (i=0;i<omp_get_max_threads();++i) {            for (i=0;i<num_threads;++i) {
84               tmp=out;               tmp=out;
85               out+=partial_sums[i];               out+=partial_sums[i];
86               partial_sums[i]=tmp;               partial_sums[i]=tmp;
# Line 82  if (omp_get_max_threads()>1) { Line 89  if (omp_get_max_threads()>1) {
89            
90        #pragma omp parallel private(sum,tmp,i)        #pragma omp parallel private(sum,tmp,i)
91        {        {
92          sum=partial_sums[omp_get_thread_num()];          sum=partial_sums[thread_num];
93          #pragma omp for schedule(static)          #pragma omp for schedule(static)
94          for (i=0;i<N;++i) {          for (i=0;i<N;++i) {
95            tmp=sum;            tmp=sum;
# Line 91  if (omp_get_max_threads()>1) { Line 98  if (omp_get_max_threads()>1) {
98          }          }
99        }        }
100        TMPMEMFREE(partial_sums);        TMPMEMFREE(partial_sums);
101  } else {     } else {
102        for (i=0;i<N;++i) {        for (i=0;i<N;++i) {
103           tmp=out;           tmp=out;
104           out+=array[i];           out+=array[i];
105           array[i]=tmp;           array[i]=tmp;
106        }        }
107  }     }
108     return out;     return out;
109  }  }
110    
# Line 105  index_t Paso_Util_cumsum_maskedTrue(dim_ Line 112  index_t Paso_Util_cumsum_maskedTrue(dim_
112     index_t out=0,tmp;     index_t out=0,tmp;
113     dim_t i;     dim_t i;
114     index_t *partial_sums=NULL,sum;     index_t *partial_sums=NULL,sum;
115      #ifdef _OPENMP
116       const int num_threads=omp_get_max_threads();
117  if (omp_get_max_threads()>1) {     const int thread_num=omp_get_thread_num();
118     partial_sums=TMPMEMALLOC(omp_get_max_threads(),index_t);  #else
119     #pragma omp parallel private(sum,i)     const int num_threads=1;
120     {     const int thread_num=0;
121        sum=0;  #endif
122        #pragma omp for schedule(static)  
123        for (i=0;i<N;++i) {     if (num_threads>1) {
124       if (mask[i]) {        partial_sums=TMPMEMALLOC(num_threads, index_t);
125          array[i] =1;        #pragma omp parallel private(sum,i)
126          sum++;        {
127       } else {           sum=0;
128          array[i] =0;           #pragma omp for schedule(static)
129       }           for (i=0;i<N;++i) {
130                if (mask[i]) {
131                   array[i] =1;
132                   sum++;
133                } else {
134                   array[i] =0;
135                }
136             }
137             partial_sums[thread_num]=sum;
138        }        }
139        partial_sums[omp_get_thread_num()]=sum;        
140     }        {
141               out=0;
142     {           for (i=0;i<num_threads;++i) {
143        out=0;              tmp=out;
144        for (i=0;i<omp_get_max_threads();++i) {              out+=partial_sums[i];
145       tmp=out;              partial_sums[i]=tmp;
146       out+=partial_sums[i];           }
147       partial_sums[i]=tmp;        }
148        }        
149     }        #pragma omp parallel private(sum,tmp,i)
150            {
151     #pragma omp parallel private(sum,tmp,i)           sum=partial_sums[thread_num];
152     {           #pragma omp for schedule(static)
153        sum=partial_sums[omp_get_thread_num()];           for (i=0;i<N;++i) {
154        #pragma omp for schedule(static)              if (mask[i]) {
155                   tmp=sum;
156                   sum+=array[i];
157                   array[i]=tmp;
158                } else {
159                   array[i]=-1;
160                }
161             }
162          }
163          TMPMEMFREE(partial_sums);
164       } else { /* num_threads=1 */
165        for (i=0;i<N;++i) {        for (i=0;i<N;++i) {
166       if (mask[i]) {           if (mask[i]) {
167          tmp=sum;              array[i]=out;
168          sum+=array[i];              out++;  
169          array[i]=tmp;           } else {
170       } else {              array[i]=-1;
171          array[i]=-1;           }
      }  
       }  
    }  
    TMPMEMFREE(partial_sums);  
 } else {  
    for (i=0;i<N;++i) {  
         
       if (mask[i]) {  
       array[i]=out;  
       out++;    
       } else {  
      array[i]=-1;  
172        }        }
173     }     }
 }  
174     return out;     return out;
175  }  }
176    
177  index_t Paso_Util_cumsum_maskedFalse(dim_t N,index_t* array, bool_t* mask) {  index_t Paso_Util_cumsum_maskedFalse(dim_t N,index_t* array, bool_t* mask) {
178     index_t out=0,tmp;     index_t out=0,tmp;
179     dim_t i;     dim_t i;
   
180     index_t *partial_sums=NULL,sum;     index_t *partial_sums=NULL,sum;
181    #ifdef _OPENMP
182       const int num_threads=omp_get_max_threads();
183       const int thread_num=omp_get_thread_num();
184    #else
185       const int num_threads=1;
186       const int thread_num=0;
187    #endif
188    
189         if (num_threads>1) {
190  if (omp_get_max_threads()>1) {        partial_sums=TMPMEMALLOC(num_threads,index_t);
191     partial_sums=TMPMEMALLOC(omp_get_max_threads(),index_t);        #pragma omp parallel private(sum,i)
192     #pragma omp parallel private(sum,i)        {
193     {           sum=0;
194        sum=0;           #pragma omp for schedule(static)
195        #pragma omp for schedule(static)           for (i=0;i<N;++i) {
196        for (i=0;i<N;++i) {              if (! mask[i]) {
197       if (! mask[i]) {                 array[i] =1;
198          array[i] =1;                 sum++;
199          sum++;              } else {
200       } else {                 array[i] =0;
201          array[i] =0;              }
202       }           }
203             partial_sums[thread_num]=sum;
204        }        }
205        partial_sums[omp_get_thread_num()]=sum;        
206     }        {
207               out=0;
208     {           for (i=0;i<num_threads;++i) {
209        out=0;              tmp=out;
210        for (i=0;i<omp_get_max_threads();++i) {              out+=partial_sums[i];
211       tmp=out;              partial_sums[i]=tmp;
212       out+=partial_sums[i];            }
213       partial_sums[i]=tmp;        }
214        }        
215     }        #pragma omp parallel private(sum,tmp,i)
216            {
217     #pragma omp parallel private(sum,tmp,i)           sum=partial_sums[thread_num];
218     {           #pragma omp for schedule(static)
219        sum=partial_sums[omp_get_thread_num()];           for (i=0;i<N;++i) {
220        #pragma omp for schedule(static)              if (! mask[i]) {
221                   tmp=sum;
222                   sum+=array[i];
223                   array[i]=tmp;
224                } else {
225                   array[i]=-1;
226                }
227             }
228          }
229          TMPMEMFREE(partial_sums);
230       } else {
231        for (i=0;i<N;++i) {        for (i=0;i<N;++i) {
232       if (! mask[i]) {           if (! mask[i]) {
233          tmp=sum;              array[i]=out;
234          sum+=array[i];              out++;  
235          array[i]=tmp;           } else {
236       } else {              array[i]=-1;
237          array[i]=-1;           }
      }  
       }  
    }  
    TMPMEMFREE(partial_sums);  
 } else {  
    for (i=0;i<N;++i) {  
       if (! mask[i]) {  
      array[i]=out;  
      out++;    
       } else {  
      array[i]=-1;  
238        }        }
239     }     }
 }  
240    
241     return out;     return out;
242  }  }
# Line 231  index_t Paso_Util_arg_max(dim_t n, dim_t Line 248  index_t Paso_Util_arg_max(dim_t n, dim_t
248     index_t argmax=-1;     index_t argmax=-1;
249     index_t lmax=-1;     index_t lmax=-1;
250     index_t li=-1;     index_t li=-1;
251    #ifdef _OPENMP
252       const int num_threads=omp_get_max_threads();
253    #else
254       const int num_threads=1;
255    #endif
256        
257     if (n>0) {     if (n>0) {
         
258        max=lambda[0];        max=lambda[0];
259        argmax=0;        argmax=0;
260        if (omp_get_max_threads()>1) {        if (num_threads>1) {
261       #pragma omp parallel private(i,lmax,li)       #pragma omp parallel private(i,lmax,li)
262       {       {
263          lmax=max;          lmax=max;
# Line 276  index_t Paso_Util_arg_max(dim_t n, dim_t Line 297  index_t Paso_Util_arg_max(dim_t n, dim_t
297  void Paso_zeroes(const dim_t n, double* x)  void Paso_zeroes(const dim_t n, double* x)
298  {  {
299     dim_t i,local_n,rest,n_start,n_end,q;     dim_t i,local_n,rest,n_start,n_end,q;
300    #ifdef _OPENMP
301     const int num_threads=omp_get_max_threads();     const int num_threads=omp_get_max_threads();
302    #else
303       const int num_threads=1;
304    #endif
305    
306     #pragma omp parallel for private(i,local_n,rest,n_start,n_end,q)     #pragma omp parallel for private(i,local_n,rest,n_start,n_end,q)
307     for (i=0;i<num_threads;++i) {     for (i=0;i<num_threads;++i) {

Legend:
Removed from v.3317  
changed lines
  Added in v.3318

  ViewVC Help
Powered by ViewVC 1.1.26