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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (show annotations)
Mon Jan 20 03:37:18 2020 UTC (4 months, 1 week ago) by uqaeller
File size: 13036 byte(s)
Updated the copyright header.


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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/4.0fordebian/finley/src/Util.cpp:5567-5588 /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 /branches/trilinos_from_5897/finley/src/Util.cpp:5898-6118 /release/3.0/finley/src/Util.cpp:2591-2601 /release/4.0/finley/src/Util.cpp:5380-5406 /trunk/finley/src/Util.cpp:4257-4344

  ViewVC Help
Powered by ViewVC 1.1.26