/[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 3282 by jfenwick, Mon Oct 11 01:48:14 2010 UTC revision 3283 by gross, Mon Oct 18 22:39:28 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     #ifdef _OPENMP  if (omp_get_max_threads()>1) {
63        index_t *partial_sums=NULL,sum;        index_t *partial_sums=NULL,sum;
64        partial_sums=TMPMEMALLOC(omp_get_max_threads(),index_t);        partial_sums=TMPMEMALLOC(omp_get_max_threads(),index_t);
65        #pragma omp parallel private(sum,i)        #pragma omp parallel private(sum,i)
# Line 71  index_t Paso_Util_cumsum(dim_t N,index_t Line 71  index_t Paso_Util_cumsum(dim_t N,index_t
71          partial_sums[omp_get_thread_num()]=sum;          partial_sums[omp_get_thread_num()]=sum;
72        }        }
73    
74          {        {
75            out=0;            out=0;
76            for (i=0;i<omp_get_max_threads();++i) {            for (i=0;i<omp_get_max_threads();++i) {
77               tmp=out;               tmp=out;
78               out+=partial_sums[i];               out+=partial_sums[i];
79               partial_sums[i]=tmp;               partial_sums[i]=tmp;
80             }             }
81          }        }
82            
83        #pragma omp parallel private(sum,tmp,i)        #pragma omp parallel private(sum,tmp,i)
84        {        {
# Line 91  index_t Paso_Util_cumsum(dim_t N,index_t Line 91  index_t Paso_Util_cumsum(dim_t N,index_t
91          }          }
92        }        }
93        TMPMEMFREE(partial_sums);        TMPMEMFREE(partial_sums);
94     #else  } else {
95        for (i=0;i<N;++i) {        for (i=0;i<N;++i) {
96           tmp=out;           tmp=out;
97           out+=array[i];           out+=array[i];
98           array[i]=tmp;           array[i]=tmp;
99        }        }
100     #endif  }
101       return out;
102    }
103    
104    index_t Paso_Util_cumsum_maskedTrue(dim_t N,index_t* array, bool_t* mask) {
105       index_t out=0,tmp;
106       dim_t i;
107       index_t *partial_sums=NULL,sum;
108      
109    
110    if (omp_get_max_threads()>1) {
111       partial_sums=TMPMEMALLOC(omp_get_max_threads(),index_t);
112       #pragma omp parallel private(sum,i)
113       {
114          sum=0;
115          #pragma omp for schedule(static)
116          for (i=0;i<N;++i) {
117         if (mask[i]) {
118            array[i] =1;
119            sum++;
120         } else {
121            array[i] =0;
122         }
123          }
124          partial_sums[omp_get_thread_num()]=sum;
125       }
126      
127       {
128          out=0;
129          for (i=0;i<omp_get_max_threads();++i) {
130         tmp=out;
131         out+=partial_sums[i];
132         partial_sums[i]=tmp;
133          }
134       }
135      
136       #pragma omp parallel private(sum,tmp,i)
137       {
138          sum=partial_sums[omp_get_thread_num()];
139          #pragma omp for schedule(static)
140          for (i=0;i<N;++i) {
141         if (mask[i]) {
142            tmp=sum;
143            sum+=array[i];
144            array[i]=tmp;
145         } else {
146            array[i]=-1;
147         }
148          }
149       }
150       TMPMEMFREE(partial_sums);
151    } else {
152       for (i=0;i<N;++i) {
153          
154          if (mask[i]) {
155          array[i]=out;
156          out++;  
157          } else {
158         array[i]=-1;
159          }
160       }
161    }
162       return out;
163    }
164    
165    index_t Paso_Util_cumsum_maskedFalse(dim_t N,index_t* array, bool_t* mask) {
166       index_t out=0,tmp;
167       dim_t i;
168    
169       index_t *partial_sums=NULL,sum;
170    
171      
172    if (omp_get_max_threads()>1) {
173       partial_sums=TMPMEMALLOC(omp_get_max_threads(),index_t);
174       #pragma omp parallel private(sum,i)
175       {
176          sum=0;
177          #pragma omp for schedule(static)
178          for (i=0;i<N;++i) {
179         if (! mask[i]) {
180            array[i] =1;
181            sum++;
182         } else {
183            array[i] =0;
184         }
185          }
186          partial_sums[omp_get_thread_num()]=sum;
187       }
188      
189       {
190          out=0;
191          for (i=0;i<omp_get_max_threads();++i) {
192         tmp=out;
193         out+=partial_sums[i];
194         partial_sums[i]=tmp;
195          }
196       }
197      
198       #pragma omp parallel private(sum,tmp,i)
199       {
200          sum=partial_sums[omp_get_thread_num()];
201          #pragma omp for schedule(static)
202          for (i=0;i<N;++i) {
203         if (! mask[i]) {
204            tmp=sum;
205            sum+=array[i];
206            array[i]=tmp;
207         } else {
208            array[i]=-1;
209         }
210          }
211       }
212       TMPMEMFREE(partial_sums);
213    } else {
214       for (i=0;i<N;++i) {
215          if (! mask[i]) {
216         array[i]=out;
217         out++;  
218          } else {
219         array[i]=-1;
220          }
221       }
222    }
223    
224     return out;     return out;
225  }  }
226    
227    /* return the index  to the largest entry in lambda takes is maximum */
228    index_t Paso_Util_arg_max(dim_t n, dim_t* lambda) {
229       dim_t i;
230       index_t max=-1;
231       index_t argmax=-1;
232       index_t lmax=-1;
233       index_t li=-1;
234      
235       if (n>0) {
236          
237          max=lambda[0];
238          argmax=0;
239          if (omp_get_max_threads()>1) {
240         #pragma omp parallel private(i,lmax,li)
241         {
242            lmax=max;
243            li=lmax;
244            #pragma omp for schedule(static)
245            for (i=0;i<n;++i) {
246               if(lmax<lambda[i]){
247              lmax=lambda[i];
248              li=i;
249               }
250            }
251            #pragma omp critical
252            {
253               if (max<=lmax) {
254              if(max==lmax && argmax>li)
255              {
256                 argmax=li;
257              }
258              if (max < lmax)
259              {
260                 max=lmax;
261                 argmax=li;
262              }
263               }
264            }
265         }
266          } else {
267         for (i=0;i<n;++i) {
268            if(max<lambda[i]){
269               max=lambda[i];
270               argmax=i;
271            }
272         }
273          }
274       }
275       return argmax;
276    }
277    
278  /*  /*
279   * set x to zero:   * set x to zero:
280   */   */
Legend:
Removed from v.3282  
changed lines
  Added in v.3283
  ViewVC Help
Powered by ViewVC 1.1.26