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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4789 - (show annotations)
Fri Mar 21 00:30:32 2014 UTC (5 years, 2 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
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2014 by University of Queensland
5 * http://www.uq.edu.au
6 *
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 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17
18 /****************************************************************************
19
20 Some utility routines
21
22 *****************************************************************************/
23
24 #include "Finley.h"
25 #include "Util.h"
26 #include "esysUtils/index.h"
27
28 #include <algorithm> // std::sort
29 #include <limits>
30
31 namespace finley {
32 namespace util {
33
34 /// comparison function for sortValueAndIndex
35 bool ValueAndIndexCompare(const std::pair<int,int> &i, const std::pair<int, int> &j)
36 {
37 // 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 }
42
43 /// orders a ValueAndIndexList by value
44 void sortValueAndIndex(ValueAndIndexList& array)
45 {
46 std::sort(array.begin(), array.end(), ValueAndIndexCompare);
47 }
48
49 /// 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 }
58 }
59
60 /// 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 }
73
74 /// multiplies two matrices: A(1:A1,1:A2) := B(1:A1,1:B2)*C(1:B2,1:A2)
75 void smallMatMult(int A1, int A2, double* A, int B2,
76 const std::vector<double>& B,
77 const std::vector<double>& C)
78 {
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 }
87 }
88
89 /// 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 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 {
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 }
105 }
106
107 /// 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 setError(ZERO_DIVISION_ERROR, "InvertSmallMat: Non-regular matrix");
120 break;
121 }
122 }
123 break;
124
125 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
132 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 setError(ZERO_DIVISION_ERROR, "InvertSmallMat: Non-regular matrix");
141 break;
142 }
143 }
144 break;
145
146 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
158 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 setError(ZERO_DIVISION_ERROR, "InvertSmallMat: Non-regular matrix");
172 break;
173 }
174 }
175 break;
176
177 default:
178 setError(VALUE_ERROR, "InvertSmallMat: dim must be <=3");
179 break;
180 }
181 }
182
183 /// 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
190 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 if (length <= 0) {
200 setError(ZERO_DIVISION_ERROR, __FILE__ ": area equals zero.");
201 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 }
208 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 if (length <= 0) {
222 setError(ZERO_DIVISION_ERROR, __FILE__ ": area equals zero.");
223 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 }
234
235 /// calculates the minimum value from a dim X N integer array
236 int getMinInt(int dim, int N, const int* values)
237 {
238 int out = std::numeric_limits<int>::max();
239 if (values && dim*N > 0) {
240 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 }
255
256 /// calculates the maximum value from a dim X N integer array
257 int getMaxInt(int dim, int N, const int* values)
258 {
259 int out = std::numeric_limits<int>::min();
260 if (values && dim*N > 0) {
261 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 }
276
277 std::pair<int,int> getMinMaxInt(int dim, int N, const int* values)
278 {
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 {
296 vmin=std::min(vmin_local, vmin);
297 vmax=std::max(vmax_local, vmax);
298 }
299 }
300 }
301 return std::pair<int,int>(vmin,vmax);
302 }
303
304 /// 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 {
308 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 #pragma omp parallel
313 {
314 int vmin_local=vmin;
315 int vmax_local=vmax;
316 #pragma omp for
317 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 }
322 }
323 #pragma omp critical
324 {
325 vmin=std::min(vmin_local, vmin);
326 vmax=std::max(vmax_local, vmax);
327 }
328 }
329 }
330 return std::pair<int,int>(vmin,vmax);
331 }
332
333 /// determines the indices of the positive entries in mask returning the
334 /// length of index.
335 std::vector<int> packMask(const std::vector<short>& mask)
336 {
337 std::vector<int> index;
338 for (int k=0; k<mask.size(); k++) {
339 if (mask[k] >= 0) {
340 index.push_back(k);
341 }
342 }
343 return index;
344 }
345
346 void setValuesInUse(const int *values, const int numValues,
347 std::vector<int>& valuesInUse, Esys_MPIInfo* mpiinfo)
348 {
349 int lastFoundValue=INDEX_T_MIN;
350 bool allFound=false;
351
352 valuesInUse.clear();
353
354 while (!allFound) {
355 // find smallest value bigger than lastFoundValue
356 int minFoundValue = INDEX_T_MAX;
357 #pragma omp parallel
358 {
359 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 }
366 #pragma omp critical
367 {
368 if (local_minFoundValue<minFoundValue)
369 minFoundValue=local_minFoundValue;
370 }
371 }
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
377 // 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 }
385 }
386
387 } // namespace util
388 } // namespace finley
389

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