/[escript]/trunk/finley/src/Util.cpp
ViewVC logotype

Annotation of /trunk/finley/src/Util.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4789 - (hide annotations)
Fri Mar 21 00:30:32 2014 UTC (5 years, 3 months ago) by sshaw
File size: 12750 byte(s)
changing more awkward logic into more obvious (and correct, though i doubt a negative length ever happens) logic
1 jgs 82
2 jfenwick 3981 /*****************************************************************************
3 ksteube 1811 *
4 jfenwick 4657 * Copyright (c) 2003-2014 by University of Queensland
5 jfenwick 3981 * http://www.uq.edu.au
6 ksteube 1811 *
7     * Primary Business: Queensland, Australia
8     * Licensed under the Open Software License version 3.0
9     * http://www.opensource.org/licenses/osl-3.0.php
10     *
11 jfenwick 3981 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 jfenwick 4657 * Development 2012-2013 by School of Earth Sciences
13     * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 jfenwick 3981 *
15     *****************************************************************************/
16 jgs 82
17 ksteube 1811
18 caltinay 4431 /****************************************************************************
19 jgs 82
20 caltinay 4431 Some utility routines
21 jgs 82
22 caltinay 4431 *****************************************************************************/
23 jgs 82
24     #include "Finley.h"
25     #include "Util.h"
26 caltinay 4441 #include "esysUtils/index.h"
27 woo409 757
28 caltinay 4441 #include <algorithm> // std::sort
29 caltinay 4428 #include <limits>
30 jfenwick 3259
31 caltinay 4441 namespace finley {
32     namespace util {
33 jgs 82
34 caltinay 4496 /// comparison function for sortValueAndIndex
35 caltinay 4441 bool ValueAndIndexCompare(const std::pair<int,int> &i, const std::pair<int, int> &j)
36     {
37 sshaw 4672 // to ensure we have a strict ordering as required by std
38     if (i.first == j.first)
39     return i.second < j.second;
40     return i.first < j.first;
41 jgs 147 }
42    
43 caltinay 4441 /// orders a ValueAndIndexList by value
44     void sortValueAndIndex(ValueAndIndexList& array)
45     {
46     std::sort(array.begin(), array.end(), ValueAndIndexCompare);
47 jgs 82 }
48    
49 caltinay 4441 /// gathers values into vector out from vector in using index:
50     /// out(1:numData, 1:len) := in(1:numData, index(1:len))
51     void gather(int len, const int* index, int numData, const double* in, double* out)
52     {
53     for (int s=0; s<len; s++) {
54     for (int i=0; i<numData; i++) {
55     out[INDEX2(i,s,numData)] = in[INDEX2(i,index[s],numData)];
56     }
57 jgs 82 }
58     }
59    
60 caltinay 4441 /// adds a vector in into out using an index:
61     /// out(1:numData,index[p])+=in(1:numData,p) where
62     /// p={k=1...len, index[k]<upperBound}
63     void addScatter(const int len, const int* index, const int numData, const double* in, double* out, const int upperBound)
64     {
65     for (int s=0; s<len; s++) {
66     for (int i=0; i<numData; i++) {
67     if (index[s] < upperBound) {
68     out[INDEX2(i,index[s],numData)]+=in[INDEX2(i,s,numData)];
69     }
70     }
71     }
72 gross 798 }
73 bcumming 751
74 caltinay 4441 /// multiplies two matrices: A(1:A1,1:A2) := B(1:A1,1:B2)*C(1:B2,1:A2)
75 caltinay 4499 void smallMatMult(int A1, int A2, double* A, int B2,
76     const std::vector<double>& B,
77     const std::vector<double>& C)
78 caltinay 4441 {
79     for (int i=0; i<A1; i++) {
80     for (int j=0; j<A2; j++) {
81     double sum=0.;
82     for (int s=0; s<B2; s++)
83     sum+=B[INDEX2(i,s,A1)]*C[INDEX2(s,j,B2)];
84     A[INDEX2(i,j,A1)]=sum;
85     }
86 jgs 82 }
87     }
88 caltinay 4431
89 caltinay 4441 /// multiplies a set of matrices with a single matrix:
90     /// A(1:A1,1:A2,i)=B(1:A1,1:B2,i)*C(1:B2,1:A2) for i=1,len
91 caltinay 4499 void smallMatSetMult1(int len, int A1, int A2, double* A, int B2,
92     const std::vector<double>& B,
93     const std::vector<double>& C)
94 caltinay 4441 {
95     for (int q=0; q<len; q++) {
96     for (int i=0; i<A1; i++) {
97     for (int j=0; j<A2; j++) {
98     double sum=0.;
99     for (int s=0; s<B2; s++)
100     sum+=B[INDEX3(i,s,q,A1,B2)]*C[INDEX2(s,j,B2)];
101     A[INDEX3(i,j,q,A1,A2)]=sum;
102     }
103     }
104 gross 2748 }
105     }
106 jgs 82
107 caltinay 4441 /// inverts the set of dim x dim matrices A(:,:,1:len) with dim=1,2,3
108     /// the inverse and determinant are returned.
109     void invertSmallMat(int len, int dim, const double* A, double *invA, double* det)
110     {
111     switch(dim) {
112     case 1:
113     for (int q=0; q<len; q++) {
114     const double D=A[q];
115     if (ABS(D) > 0) {
116     det[q]=D;
117     invA[q]=1./D;
118     } else {
119 caltinay 4496 setError(ZERO_DIVISION_ERROR, "InvertSmallMat: Non-regular matrix");
120 caltinay 4441 break;
121     }
122 jgs 82 }
123 caltinay 4441 break;
124 jgs 82
125 caltinay 4441 case 2:
126     for (int q=0; q<len; q++) {
127     const double A11=A[INDEX3(0,0,q,2,2)];
128     const double A12=A[INDEX3(0,1,q,2,2)];
129     const double A21=A[INDEX3(1,0,q,2,2)];
130     const double A22=A[INDEX3(1,1,q,2,2)];
131 jgs 82
132 caltinay 4441 const double D = A11*A22-A12*A21;
133     if (ABS(D) > 0) {
134     det[q]=D;
135     invA[INDEX3(0,0,q,2,2)]= A22/D;
136     invA[INDEX3(1,0,q,2,2)]=-A21/D;
137     invA[INDEX3(0,1,q,2,2)]=-A12/D;
138     invA[INDEX3(1,1,q,2,2)]= A11/D;
139     } else {
140 caltinay 4496 setError(ZERO_DIVISION_ERROR, "InvertSmallMat: Non-regular matrix");
141 caltinay 4441 break;
142     }
143 jgs 82 }
144 caltinay 4441 break;
145 jgs 82
146 caltinay 4441 case 3:
147     for (int q=0; q<len; q++) {
148     const double A11=A[INDEX3(0,0,q,3,3)];
149     const double A21=A[INDEX3(1,0,q,3,3)];
150     const double A31=A[INDEX3(2,0,q,3,3)];
151     const double A12=A[INDEX3(0,1,q,3,3)];
152     const double A22=A[INDEX3(1,1,q,3,3)];
153     const double A32=A[INDEX3(2,1,q,3,3)];
154     const double A13=A[INDEX3(0,2,q,3,3)];
155     const double A23=A[INDEX3(1,2,q,3,3)];
156     const double A33=A[INDEX3(2,2,q,3,3)];
157 jgs 82
158 caltinay 4441 const double D = A11*(A22*A33-A23*A32) + A12*(A31*A23-A21*A33) + A13*(A21*A32-A31*A22);
159     if (ABS(D) > 0) {
160     det[q]=D;
161     invA[INDEX3(0,0,q,3,3)]=(A22*A33-A23*A32)/D;
162     invA[INDEX3(1,0,q,3,3)]=(A31*A23-A21*A33)/D;
163     invA[INDEX3(2,0,q,3,3)]=(A21*A32-A31*A22)/D;
164     invA[INDEX3(0,1,q,3,3)]=(A13*A32-A12*A33)/D;
165     invA[INDEX3(1,1,q,3,3)]=(A11*A33-A31*A13)/D;
166     invA[INDEX3(2,1,q,3,3)]=(A12*A31-A11*A32)/D;
167     invA[INDEX3(0,2,q,3,3)]=(A12*A23-A13*A22)/D;
168     invA[INDEX3(1,2,q,3,3)]=(A13*A21-A11*A23)/D;
169     invA[INDEX3(2,2,q,3,3)]=(A11*A22-A12*A21)/D;
170     } else {
171 caltinay 4496 setError(ZERO_DIVISION_ERROR, "InvertSmallMat: Non-regular matrix");
172 caltinay 4441 break;
173     }
174 jgs 82 }
175 caltinay 4441 break;
176 jgs 82
177 caltinay 4441 default:
178 caltinay 4496 setError(VALUE_ERROR, "InvertSmallMat: dim must be <=3");
179 caltinay 4441 break;
180     }
181 jgs 82 }
182    
183 caltinay 4441 /// returns the normalized vector Normal[dim,len] orthogonal to A(:,0,q) and
184     /// A(:,1,q) in the case of dim=3, or the vector A(:,0,q) in the case of dim=2
185     void normalVector(int len, int dim, int dim1, const double* A, double* Normal)
186     {
187     int q;
188     double A11,A12,CO_A13,A21,A22,CO_A23,A31,A32,CO_A33,length,invlength;
189 jgs 82
190 caltinay 4441 switch(dim) {
191     case 1:
192     for (q=0;q<len;q++) Normal[q]=1;
193     break;
194     case 2:
195     for (q=0;q<len;q++) {
196     A11=A[INDEX3(0,0,q,2,dim1)];
197     A21=A[INDEX3(1,0,q,2,dim1)];
198     length = sqrt(A11*A11+A21*A21);
199 sshaw 4789 if (length <= 0) {
200 caltinay 4496 setError(ZERO_DIVISION_ERROR, __FILE__ ": area equals zero.");
201 caltinay 4441 return;
202     } else {
203     invlength=1./length;
204     Normal[INDEX2(0,q,2)]=A21*invlength;
205     Normal[INDEX2(1,q,2)]=-A11*invlength;
206     }
207 jgs 82 }
208 caltinay 4441 break;
209     case 3:
210     for (q=0;q<len;q++) {
211     A11=A[INDEX3(0,0,q,3,dim1)];
212     A21=A[INDEX3(1,0,q,3,dim1)];
213     A31=A[INDEX3(2,0,q,3,dim1)];
214     A12=A[INDEX3(0,1,q,3,dim1)];
215     A22=A[INDEX3(1,1,q,3,dim1)];
216     A32=A[INDEX3(2,1,q,3,dim1)];
217     CO_A13=A21*A32-A31*A22;
218     CO_A23=A31*A12-A11*A32;
219     CO_A33=A11*A22-A21*A12;
220     length=sqrt(CO_A13*CO_A13+CO_A23*CO_A23+CO_A33*CO_A33);
221 sshaw 4789 if (length <= 0) {
222 caltinay 4496 setError(ZERO_DIVISION_ERROR, __FILE__ ": area equals zero.");
223 caltinay 4441 return;
224     } else {
225     invlength=1./length;
226     Normal[INDEX2(0,q,3)]=CO_A13*invlength;
227     Normal[INDEX2(1,q,3)]=CO_A23*invlength;
228     Normal[INDEX2(2,q,3)]=CO_A33*invlength;
229     }
230     }
231     break;
232     }
233 jgs 82 }
234    
235 caltinay 4431 /// calculates the minimum value from a dim X N integer array
236 caltinay 4441 int getMinInt(int dim, int N, const int* values)
237 caltinay 4431 {
238     int out = std::numeric_limits<int>::max();
239 caltinay 4441 if (values && dim*N > 0) {
240 caltinay 4431 out=values[0];
241     #pragma omp parallel
242     {
243     int out_local=out;
244     #pragma omp for
245     for (int j=0; j<N; j++) {
246     for (int i=0; i<dim; i++)
247     out_local=std::min(out_local, values[INDEX2(i,j,dim)]);
248     }
249     #pragma omp critical
250     out=std::min(out_local, out);
251     }
252     }
253     return out;
254 jgs 82 }
255 caltinay 4431
256     /// calculates the maximum value from a dim X N integer array
257 caltinay 4441 int getMaxInt(int dim, int N, const int* values)
258 caltinay 4431 {
259     int out = std::numeric_limits<int>::min();
260 caltinay 4441 if (values && dim*N > 0) {
261 caltinay 4431 out=values[0];
262     #pragma omp parallel
263     {
264     int out_local=out;
265     #pragma omp for
266     for (int j=0; j<N; j++) {
267     for (int i=0; i<dim; i++)
268     out_local=std::max(out_local, values[INDEX2(i,j,dim)]);
269     }
270     #pragma omp critical
271     out=std::max(out_local, out);
272     }
273     }
274     return out;
275 jgs 82 }
276 caltinay 4428
277 caltinay 4441 std::pair<int,int> getMinMaxInt(int dim, int N, const int* values)
278 caltinay 4428 {
279     int vmin = std::numeric_limits<int>::max();
280     int vmax = std::numeric_limits<int>::min();
281     if (values && dim*N > 0) {
282     vmin = vmax = values[0];
283     #pragma omp parallel
284     {
285     int vmin_local=vmin;
286     int vmax_local=vmax;
287     #pragma omp for
288     for (int j=0; j<N; j++) {
289     for (int i=0; i<dim; i++) {
290     vmin_local=std::min(vmin_local, values[INDEX2(i,j,dim)]);
291     vmax_local=std::max(vmax_local, values[INDEX2(i,j,dim)]);
292     }
293     }
294     #pragma omp critical
295 caltinay 4431 {
296     vmin=std::min(vmin_local, vmin);
297     vmax=std::max(vmax_local, vmax);
298     }
299 caltinay 4428 }
300     }
301     return std::pair<int,int>(vmin,vmax);
302     }
303 jgs 82
304 caltinay 4492 /// calculates the minimum and maximum value from an integer array of length N
305     /// disregarding the value 'ignore'
306     std::pair<int,int> getFlaggedMinMaxInt(int N, const int* values, int ignore)
307 caltinay 4441 {
308 caltinay 4492 int vmin = std::numeric_limits<int>::max();
309     int vmax = std::numeric_limits<int>::min();
310     if (values && N > 0) {
311     vmin = vmax = values[0];
312 caltinay 4441 #pragma omp parallel
313     {
314 caltinay 4492 int vmin_local=vmin;
315     int vmax_local=vmax;
316 caltinay 4441 #pragma omp for
317 caltinay 4492 for (int i=0; i<N; i++) {
318     if (values[i] != ignore) {
319     vmin_local=std::min(vmin_local, values[i]);
320     vmax_local=std::max(vmax_local, values[i]);
321 caltinay 4441 }
322     }
323 caltinay 4431 #pragma omp critical
324 caltinay 4492 {
325     vmin=std::min(vmin_local, vmin);
326     vmax=std::max(vmax_local, vmax);
327 caltinay 4441 }
328 jgs 82 }
329 caltinay 4441 }
330 caltinay 4492 return std::pair<int,int>(vmin,vmax);
331 jgs 82 }
332    
333 caltinay 4441 /// determines the indices of the positive entries in mask returning the
334     /// length of index.
335 caltinay 4496 std::vector<int> packMask(const std::vector<short>& mask)
336 caltinay 4441 {
337 caltinay 4496 std::vector<int> index;
338     for (int k=0; k<mask.size(); k++) {
339 caltinay 4441 if (mask[k] >= 0) {
340 caltinay 4496 index.push_back(k);
341 jgs 113 }
342 caltinay 4441 }
343 caltinay 4496 return index;
344 jgs 113 }
345 caltinay 3639
346 caltinay 4441 void setValuesInUse(const int *values, const int numValues,
347     std::vector<int>& valuesInUse, Esys_MPIInfo* mpiinfo)
348 gross 1716 {
349 caltinay 4428 int lastFoundValue=INDEX_T_MIN;
350     bool allFound=false;
351 jgs 82
352 caltinay 4428 valuesInUse.clear();
353    
354     while (!allFound) {
355 caltinay 4431 // find smallest value bigger than lastFoundValue
356 caltinay 4428 int minFoundValue = INDEX_T_MAX;
357     #pragma omp parallel
358 gross 1716 {
359 caltinay 4428 int local_minFoundValue=minFoundValue;
360     #pragma omp for
361     for (int i=0; i<numValues; i++) {
362     const int val=values[i];
363     if ((val>lastFoundValue) && (val<local_minFoundValue))
364     local_minFoundValue=val;
365 gross 1716 }
366 caltinay 4428 #pragma omp critical
367 gross 2425 {
368 caltinay 4428 if (local_minFoundValue<minFoundValue)
369     minFoundValue=local_minFoundValue;
370 gross 2425 }
371 caltinay 4428 }
372     #ifdef ESYS_MPI
373     int local_minFoundValue=minFoundValue;
374     MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, mpiinfo->comm);
375     #endif
376 gross 2425
377 caltinay 4431 // if we found a new value we need to add this to valuesInUse
378     if (minFoundValue < INDEX_T_MAX) {
379     valuesInUse.push_back(minFoundValue);
380     lastFoundValue=minFoundValue;
381     } else {
382     allFound=true;
383     }
384 caltinay 4428 }
385 gross 1716 }
386    
387 caltinay 4441 } // namespace util
388     } // namespace finley
389 gross 1716

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/lapack2681/finley/src/Util.cpp:2682-2741 /branches/pasowrap/finley/src/Util.cpp:3661-3674 /branches/py3_attempt2/finley/src/Util.cpp:3871-3891 /branches/restext/finley/src/Util.cpp:2610-2624 /branches/ripleygmg_from_3668/finley/src/Util.cpp:3669-3791 /branches/stage3.0/finley/src/Util.cpp:2569-2590 /branches/symbolic_from_3470/finley/src/Util.cpp:3471-3974 /release/3.0/finley/src/Util.cpp:2591-2601 /trunk/finley/src/Util.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26