/[escript]/branches/clazy/escriptcore/src/DataVectorOps.h
ViewVC logotype

Contents of /branches/clazy/escriptcore/src/DataVectorOps.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6517 - (show annotations)
Mon Mar 6 06:35:31 2017 UTC (2 years, 1 month ago) by jfenwick
File MIME type: text/plain
File size: 57682 byte(s)
Work towards complex lazy.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 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 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17 #ifndef __ESCRIPT_DATAMATHS_H__
18 #define __ESCRIPT_DATAMATHS_H__
19
20 #include "DataAbstract.h"
21 #include "DataException.h"
22 #include "ArrayOps.h"
23 #include "LapackInverseHelper.h"
24 #include "DataTagged.h"
25 #include <complex>
26 /**
27 \file DataVectorOps.h
28 \brief Describes binary operations performed on DataVector.
29
30
31 For operations on DataReady see BinaryDataReadyOp.h.
32 For operations on double* see ArrayOps.h.
33 */
34
35
36 namespace escript
37 {
38
39 /**
40 In order to properly identify the datapoints, in most cases, the vector, shape and offset of the point must all be supplied.
41 Note that vector in this context refers to a data vector storing datapoints not a mathematical vector. (However, datapoints within the data vector could represent scalars, vectors, matricies, ...).
42 */
43
44
45 /**
46 \brief
47 Perform a matrix multiply of the given views.
48
49 NB: Only multiplies together the two given datapoints,
50 would need to call this over all data-points to multiply the entire
51 Data objects involved.
52
53 \param left,right - vectors containing the datapoints
54 \param leftShape,rightShape - shapes of datapoints in the vectors
55 \param leftOffset,rightOffset - beginnings of datapoints in the vectors
56 \param result - Vector to store the resulting datapoint in
57 \param resultShape - expected shape of the resulting datapoint
58 */
59 ESCRIPT_DLL_API
60 void
61 matMult(const DataTypes::RealVectorType& left,
62 const DataTypes::ShapeType& leftShape,
63 DataTypes::RealVectorType::size_type leftOffset,
64 const DataTypes::RealVectorType& right,
65 const DataTypes::ShapeType& rightShape,
66 DataTypes::RealVectorType::size_type rightOffset,
67 DataTypes::RealVectorType& result,
68 const DataTypes::ShapeType& resultShape);
69 // Hmmmm why is there no offset for the result??
70
71
72
73
74 /**
75 \brief
76 Determine the shape of the result array for a matrix multiplication
77 of the given views.
78
79 \param left,right - shapes of the left and right matricies
80 \return the shape of the matrix which would result from multiplying left and right
81 */
82 ESCRIPT_DLL_API
83 DataTypes::ShapeType
84 determineResultShape(const DataTypes::ShapeType& left,
85 const DataTypes::ShapeType& right);
86
87
88 /**
89 \brief
90 computes a symmetric matrix from your square matrix A: (A + transpose(A)) / 2
91
92 \param in - vector containing the matrix A
93 \param inShape - shape of the matrix A
94 \param inOffset - the beginning of A within the vector in
95 \param ev - vector to store the output matrix
96 \param evShape - expected shape of the output matrix
97 \param evOffset - starting location for storing ev in vector ev
98 */
99 template<typename VEC>
100 inline
101 void
102 symmetric(const VEC& in,
103 const DataTypes::ShapeType& inShape,
104 typename VEC::size_type inOffset,
105 VEC& ev,
106 const DataTypes::ShapeType& evShape,
107 typename VEC::size_type evOffset)
108 {
109 if (DataTypes::getRank(inShape) == 2) {
110 int i0, i1;
111 int s0=inShape[0];
112 int s1=inShape[1];
113 for (i0=0; i0<s0; i0++) {
114 for (i1=0; i1<s1; i1++) {
115 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)] + in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)]) / 2.0;
116 }
117 }
118 }
119 else if (DataTypes::getRank(inShape) == 4) {
120 int i0, i1, i2, i3;
121 int s0=inShape[0];
122 int s1=inShape[1];
123 int s2=inShape[2];
124 int s3=inShape[3];
125 for (i0=0; i0<s0; i0++) {
126 for (i1=0; i1<s1; i1++) {
127 for (i2=0; i2<s2; i2++) {
128 for (i3=0; i3<s3; i3++) {
129 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i3)] + in[inOffset+DataTypes::getRelIndex(inShape,i2,i3,i0,i1)]) / 2.0;
130 }
131 }
132 }
133 }
134 }
135 }
136
137 /**
138 \brief
139 computes a antisymmetric matrix from your square matrix A: (A - transpose(A)) / 2
140
141 \param in - vector containing the matrix A
142 \param inShape - shape of the matrix A
143 \param inOffset - the beginning of A within the vector in
144 \param ev - vector to store the output matrix
145 \param evShape - expected shape of the output matrix
146 \param evOffset - starting location for storing ev in vector ev
147 */
148 template<typename VEC>
149 inline
150 void
151 antisymmetric(const VEC& in,
152 const DataTypes::ShapeType& inShape,
153 typename VEC::size_type inOffset,
154 VEC& ev,
155 const DataTypes::ShapeType& evShape,
156 typename VEC::size_type evOffset)
157 {
158 if (DataTypes::getRank(inShape) == 2) {
159 int i0, i1;
160 int s0=inShape[0];
161 int s1=inShape[1];
162 for (i0=0; i0<s0; i0++) {
163 for (i1=0; i1<s1; i1++) {
164 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)] - in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)]) / 2.0;
165 }
166 }
167 }
168 else if (DataTypes::getRank(inShape) == 4) {
169 int i0, i1, i2, i3;
170 int s0=inShape[0];
171 int s1=inShape[1];
172 int s2=inShape[2];
173 int s3=inShape[3];
174 for (i0=0; i0<s0; i0++) {
175 for (i1=0; i1<s1; i1++) {
176 for (i2=0; i2<s2; i2++) {
177 for (i3=0; i3<s3; i3++) {
178 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = (in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i3)] - in[inOffset+DataTypes::getRelIndex(inShape,i2,i3,i0,i1)]) / 2.0;
179 }
180 }
181 }
182 }
183 }
184 }
185
186
187
188 /**
189 \brief
190 computes an hermitian matrix from your square matrix A: (A + adjoint(A)) / 2
191
192 \param in - vector containing the matrix A
193 \param inShape - shape of the matrix A
194 \param inOffset - the beginning of A within the vector in
195 \param ev - vector to store the output matrix
196 \param evShape - expected shape of the output matrix
197 \param evOffset - starting location for storing ev in vector ev
198 */
199 void
200 hermitian(const DataTypes::CplxVectorType& in,
201 const DataTypes::ShapeType& inShape,
202 DataTypes::CplxVectorType::size_type inOffset,
203 DataTypes::CplxVectorType& ev,
204 const DataTypes::ShapeType& evShape,
205 DataTypes::CplxVectorType::size_type evOffset);
206
207 /**
208 \brief
209 computes a antihermitian matrix from your square matrix A: (A - adjoint(A)) / 2
210
211 \param in - vector containing the matrix A
212 \param inShape - shape of the matrix A
213 \param inOffset - the beginning of A within the vector in
214 \param ev - vector to store the output matrix
215 \param evShape - expected shape of the output matrix
216 \param evOffset - starting location for storing ev in vector ev
217 */
218 void
219 antihermitian(const DataTypes::CplxVectorType& in,
220 const DataTypes::ShapeType& inShape,
221 typename DataTypes::CplxVectorType::size_type inOffset,
222 DataTypes::CplxVectorType& ev,
223 const DataTypes::ShapeType& evShape,
224 typename DataTypes::CplxVectorType::size_type evOffset);
225
226 /**
227 \brief
228 computes the trace of a matrix
229
230 \param in - vector containing the input matrix
231 \param inShape - shape of the input matrix
232 \param inOffset - the beginning of the input matrix within the vector "in"
233 \param ev - vector to store the output matrix
234 \param evShape - expected shape of the output matrix
235 \param evOffset - starting location for storing the output matrix in vector ev
236 \param axis_offset
237 */
238 template <class VEC>
239 inline
240 void
241 trace(const VEC& in,
242 const DataTypes::ShapeType& inShape,
243 typename VEC::size_type inOffset,
244 VEC& ev,
245 const DataTypes::ShapeType& evShape,
246 typename VEC::size_type evOffset,
247 int axis_offset)
248 {
249 for (int j=0;j<DataTypes::noValues(evShape);++j)
250 {
251 ev[evOffset+j]=0;
252 }
253 if (DataTypes::getRank(inShape) == 2) {
254 int s0=inShape[0]; // Python wrapper limits to square matrix
255 int i;
256 for (i=0; i<s0; i++) {
257 ev[evOffset/*+DataTypes::getRelIndex(evShape)*/] += in[inOffset+DataTypes::getRelIndex(inShape,i,i)];
258 }
259 }
260 else if (DataTypes::getRank(inShape) == 3) {
261 if (axis_offset==0) {
262 int s0=inShape[0];
263 int s2=inShape[2];
264 int i0, i2;
265 for (i0=0; i0<s0; i0++) {
266 for (i2=0; i2<s2; i2++) {
267 ev[evOffset+DataTypes::getRelIndex(evShape,i2)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i0,i2)];
268 }
269 }
270 }
271 else if (axis_offset==1) {
272 int s0=inShape[0];
273 int s1=inShape[1];
274 int i0, i1;
275 for (i0=0; i0<s0; i0++) {
276 for (i1=0; i1<s1; i1++) {
277 ev[evOffset+DataTypes::getRelIndex(evShape,i0)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i1)];
278 }
279 }
280 }
281 }
282 else if (DataTypes::getRank(inShape) == 4) {
283 if (axis_offset==0) {
284 int s0=inShape[0];
285 int s2=inShape[2];
286 int s3=inShape[3];
287 int i0, i2, i3;
288 for (i0=0; i0<s0; i0++) {
289 for (i2=0; i2<s2; i2++) {
290 for (i3=0; i3<s3; i3++) {
291 ev[evOffset+DataTypes::getRelIndex(evShape,i2,i3)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i0,i2,i3)];
292 }
293 }
294 }
295 }
296 else if (axis_offset==1) {
297 int s0=inShape[0];
298 int s1=inShape[1];
299 int s3=inShape[3];
300 int i0, i1, i3;
301 for (i0=0; i0<s0; i0++) {
302 for (i1=0; i1<s1; i1++) {
303 for (i3=0; i3<s3; i3++) {
304 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i3)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i1,i3)];
305 }
306 }
307 }
308 }
309 else if (axis_offset==2) {
310 int s0=inShape[0];
311 int s1=inShape[1];
312 int s2=inShape[2];
313 int i0, i1, i2;
314 for (i0=0; i0<s0; i0++) {
315 for (i1=0; i1<s1; i1++) {
316 for (i2=0; i2<s2; i2++) {
317 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] += in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i2)];
318 }
319 }
320 }
321 }
322 }
323 }
324
325
326 /**
327 \brief
328 Transpose each data point of this Data object around the given axis.
329
330 \param in - vector containing the input matrix
331 \param inShape - shape of the input matrix
332 \param inOffset - the beginning of the input matrix within the vector "in"
333 \param ev - vector to store the output matrix
334 \param evShape - expected shape of the output matrix
335 \param evOffset - starting location for storing the output matrix in vector ev
336 \param axis_offset
337 */
338 ESCRIPT_DLL_API
339 template <class VEC>
340 inline
341 void
342 transpose(const VEC& in,
343 const DataTypes::ShapeType& inShape,
344 typename VEC::size_type inOffset,
345 VEC& ev,
346 const DataTypes::ShapeType& evShape,
347 typename VEC::size_type evOffset,
348 int axis_offset)
349 {
350 int inRank=DataTypes::getRank(inShape);
351 if ( inRank== 4) {
352 int s0=evShape[0];
353 int s1=evShape[1];
354 int s2=evShape[2];
355 int s3=evShape[3];
356 int i0, i1, i2, i3;
357 if (axis_offset==1) {
358 for (i0=0; i0<s0; i0++) {
359 for (i1=0; i1<s1; i1++) {
360 for (i2=0; i2<s2; i2++) {
361 for (i3=0; i3<s3; i3++) {
362 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i3,i0,i1,i2)];
363 }
364 }
365 }
366 }
367 }
368 else if (axis_offset==2) {
369 for (i0=0; i0<s0; i0++) {
370 for (i1=0; i1<s1; i1++) {
371 for (i2=0; i2<s2; i2++) {
372 for (i3=0; i3<s3; i3++) {
373 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i3,i0,i1)];
374 }
375 }
376 }
377 }
378 }
379 else if (axis_offset==3) {
380 for (i0=0; i0<s0; i0++) {
381 for (i1=0; i1<s1; i1++) {
382 for (i2=0; i2<s2; i2++) {
383 for (i3=0; i3<s3; i3++) {
384 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i2,i3,i0)];
385 }
386 }
387 }
388 }
389 }
390 else {
391 for (i0=0; i0<s0; i0++) {
392 for (i1=0; i1<s1; i1++) {
393 for (i2=0; i2<s2; i2++) {
394 for (i3=0; i3<s3; i3++) {
395 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2,i3)];
396 }
397 }
398 }
399 }
400 }
401 }
402 else if (inRank == 3) {
403 int s0=evShape[0];
404 int s1=evShape[1];
405 int s2=evShape[2];
406 int i0, i1, i2;
407 if (axis_offset==1) {
408 for (i0=0; i0<s0; i0++) {
409 for (i1=0; i1<s1; i1++) {
410 for (i2=0; i2<s2; i2++) {
411 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i0,i1)];
412 }
413 }
414 }
415 }
416 else if (axis_offset==2) {
417 for (i0=0; i0<s0; i0++) {
418 for (i1=0; i1<s1; i1++) {
419 for (i2=0; i2<s2; i2++) {
420 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i2,i0)];
421 }
422 }
423 }
424 }
425 else {
426 // Copy the matrix unchanged
427 for (i0=0; i0<s0; i0++) {
428 for (i1=0; i1<s1; i1++) {
429 for (i2=0; i2<s2; i2++) {
430 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i2)];
431 }
432 }
433 }
434 }
435 }
436 else if (inRank == 2) {
437 int s0=evShape[0];
438 int s1=evShape[1];
439 int i0, i1;
440 if (axis_offset==1) {
441 for (i0=0; i0<s0; i0++) {
442 for (i1=0; i1<s1; i1++) {
443 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)];
444 }
445 }
446 }
447 else {
448 for (i0=0; i0<s0; i0++) {
449 for (i1=0; i1<s1; i1++) {
450 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1)];
451 }
452 }
453 }
454 }
455 else if (inRank == 1) {
456 int s0=evShape[0];
457 int i0;
458 for (i0=0; i0<s0; i0++) {
459 ev[evOffset+DataTypes::getRelIndex(evShape,i0)] = in[inOffset+DataTypes::getRelIndex(inShape,i0)];
460 }
461 }
462 else if (inRank == 0) {
463 ev[evOffset/*+DataTypes::getRelIndex(evShape,)*/] = in[inOffset/*+DataTypes::getRelIndex(inShape,)*/];
464 }
465 else {
466 throw DataException("Error - DataArrayView::transpose can only be calculated for rank 0, 1, 2, 3 or 4 objects.");
467 }
468 }
469
470 /**
471 \brief
472 swaps the components axis0 and axis1.
473
474 \param in - vector containing the input matrix
475 \param inShape - shape of the input matrix
476 \param inOffset - the beginning of the input matrix within the vector "in"
477 \param ev - vector to store the output matrix
478 \param evShape - expected shape of the output matrix
479 \param evOffset - starting location for storing the output matrix in vector ev
480 \param axis0 - axis index
481 \param axis1 - axis index
482 */
483 ESCRIPT_DLL_API
484 template <class VEC>
485 inline
486 void
487 swapaxes(const VEC& in,
488 const DataTypes::ShapeType& inShape,
489 typename VEC::size_type inOffset,
490 VEC& ev,
491 const DataTypes::ShapeType& evShape,
492 typename VEC::size_type evOffset,
493 int axis0,
494 int axis1)
495 {
496 int inRank=DataTypes::getRank(inShape);
497 if (inRank == 4) {
498 int s0=evShape[0];
499 int s1=evShape[1];
500 int s2=evShape[2];
501 int s3=evShape[3];
502 int i0, i1, i2, i3;
503 if (axis0==0) {
504 if (axis1==1) {
505 for (i0=0; i0<s0; i0++) {
506 for (i1=0; i1<s1; i1++) {
507 for (i2=0; i2<s2; i2++) {
508 for (i3=0; i3<s3; i3++) {
509 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0,i2,i3)];
510 }
511 }
512 }
513 }
514 } else if (axis1==2) {
515 for (i0=0; i0<s0; i0++) {
516 for (i1=0; i1<s1; i1++) {
517 for (i2=0; i2<s2; i2++) {
518 for (i3=0; i3<s3; i3++) {
519 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i1,i0,i3)];
520 }
521 }
522 }
523 }
524
525 } else if (axis1==3) {
526 for (i0=0; i0<s0; i0++) {
527 for (i1=0; i1<s1; i1++) {
528 for (i2=0; i2<s2; i2++) {
529 for (i3=0; i3<s3; i3++) {
530 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i3,i1,i2,i0)];
531 }
532 }
533 }
534 }
535 }
536 } else if (axis0==1) {
537 if (axis1==2) {
538 for (i0=0; i0<s0; i0++) {
539 for (i1=0; i1<s1; i1++) {
540 for (i2=0; i2<s2; i2++) {
541 for (i3=0; i3<s3; i3++) {
542 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i2,i1,i3)];
543 }
544 }
545 }
546 }
547 } else if (axis1==3) {
548 for (i0=0; i0<s0; i0++) {
549 for (i1=0; i1<s1; i1++) {
550 for (i2=0; i2<s2; i2++) {
551 for (i3=0; i3<s3; i3++) {
552 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i3,i2,i1)];
553 }
554 }
555 }
556 }
557 }
558 } else if (axis0==2) {
559 if (axis1==3) {
560 for (i0=0; i0<s0; i0++) {
561 for (i1=0; i1<s1; i1++) {
562 for (i2=0; i2<s2; i2++) {
563 for (i3=0; i3<s3; i3++) {
564 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2,i3)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i1,i3,i2)];
565 }
566 }
567 }
568 }
569 }
570 }
571
572 } else if ( inRank == 3) {
573 int s0=evShape[0];
574 int s1=evShape[1];
575 int s2=evShape[2];
576 int i0, i1, i2;
577 if (axis0==0) {
578 if (axis1==1) {
579 for (i0=0; i0<s0; i0++) {
580 for (i1=0; i1<s1; i1++) {
581 for (i2=0; i2<s2; i2++) {
582 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0,i2)];
583 }
584 }
585 }
586 } else if (axis1==2) {
587 for (i0=0; i0<s0; i0++) {
588 for (i1=0; i1<s1; i1++) {
589 for (i2=0; i2<s2; i2++) {
590 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i2,i1,i0)];
591 }
592 }
593 }
594 }
595 } else if (axis0==1) {
596 if (axis1==2) {
597 for (i0=0; i0<s0; i0++) {
598 for (i1=0; i1<s1; i1++) {
599 for (i2=0; i2<s2; i2++) {
600 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1,i2)] = in[inOffset+DataTypes::getRelIndex(inShape,i0,i2,i1)];
601 }
602 }
603 }
604 }
605 }
606 } else if ( inRank == 2) {
607 int s0=evShape[0];
608 int s1=evShape[1];
609 int i0, i1;
610 if (axis0==0) {
611 if (axis1==1) {
612 for (i0=0; i0<s0; i0++) {
613 for (i1=0; i1<s1; i1++) {
614 ev[evOffset+DataTypes::getRelIndex(evShape,i0,i1)] = in[inOffset+DataTypes::getRelIndex(inShape,i1,i0)];
615 }
616 }
617 }
618 }
619 } else {
620 throw DataException("Error - DataArrayView::swapaxes can only be calculated for rank 2, 3 or 4 objects.");
621 }
622 }
623
624 /**
625 \brief
626 solves a local eigenvalue problem
627
628 \param in - vector containing the input matrix
629 \param inShape - shape of the input matrix
630 \param inOffset - the beginning of the input matrix within the vector "in"
631 \param ev - vector to store the eigenvalues
632 \param evShape - expected shape of the eigenvalues
633 \param evOffset - starting location for storing the eigenvalues in vector ev
634 */
635 ESCRIPT_DLL_API
636 inline
637 void
638 eigenvalues(const DataTypes::RealVectorType& in,
639 const DataTypes::ShapeType& inShape,
640 typename DataTypes::RealVectorType::size_type inOffset,
641 DataTypes::RealVectorType& ev,
642 const DataTypes::ShapeType& evShape,
643 typename DataTypes::RealVectorType::size_type evOffset)
644 {
645 typename DataTypes::RealVectorType::ElementType in00,in10,in20,in01,in11,in21,in02,in12,in22;
646 typename DataTypes::RealVectorType::ElementType ev0,ev1,ev2;
647 int s=inShape[0];
648 if (s==1) {
649 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
650 eigenvalues1(in00,&ev0);
651 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
652
653 } else if (s==2) {
654 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
655 in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
656 in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
657 in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
658 eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
659 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
660 ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
661
662 } else if (s==3) {
663 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
664 in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
665 in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
666 in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
667 in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
668 in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
669 in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
670 in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
671 in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
672 eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
673 &ev0,&ev1,&ev2);
674 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
675 ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
676 ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
677
678 }
679 }
680
681 inline
682 void
683 eigenvalues(const DataTypes::CplxVectorType& in,
684 const DataTypes::ShapeType& inShape,
685 typename DataTypes::CplxVectorType::size_type inOffset,
686 DataTypes::CplxVectorType& ev,
687 const DataTypes::ShapeType& evShape,
688 typename DataTypes::CplxVectorType::size_type evOffset)
689 {
690 typename DataTypes::CplxVectorType::ElementType in00,in10,in20,in01,in11,in21,in02,in12,in22;
691 typename DataTypes::CplxVectorType::ElementType ev0,ev1,ev2;
692 int s=inShape[0];
693 if (s==1) {
694 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
695 eigenvalues1(in00,&ev0);
696 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
697
698 } else if (s==2) {
699 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
700 in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
701 in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
702 in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
703 eigenvalues2(in00,(in01+in10)/2.,in11,&ev0,&ev1);
704 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
705 ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
706
707 } else if (s==3) {
708 // this doesn't work yet
709 // in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
710 // in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
711 // in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
712 // in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
713 // in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
714 // in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
715 // in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
716 // in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
717 // in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
718 // eigenvalues3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
719 // &ev0,&ev1,&ev2);
720 // ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
721 // ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
722 // ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
723
724 }
725 }
726
727
728 /**
729 \brief
730 solves a local eigenvalue problem
731
732 \param in - vector containing the input matrix
733 \param inShape - shape of the input matrix
734 \param inOffset - the beginning of the input matrix within the vector "in"
735 \param ev - vector to store the eigenvalues
736 \param evShape - expected shape of the eigenvalues
737 \param evOffset - starting location for storing the eigenvalues in ev
738 \param V - vector to store the eigenvectors
739 \param VShape - expected shape of the eigenvectors
740 \param VOffset - starting location for storing the eigenvectors in V
741 \param tol - Input - eigenvalues with relative difference tol are treated as equal
742 */
743 ESCRIPT_DLL_API
744 inline
745 void
746 eigenvalues_and_eigenvectors(const DataTypes::RealVectorType& in, const DataTypes::ShapeType& inShape,
747 DataTypes::RealVectorType::size_type inOffset,
748 DataTypes::RealVectorType& ev, const DataTypes::ShapeType& evShape,
749 DataTypes::RealVectorType::size_type evOffset,
750 DataTypes::RealVectorType& V, const DataTypes::ShapeType& VShape,
751 DataTypes::RealVectorType::size_type VOffset,
752 const double tol=1.e-13)
753 {
754 double in00,in10,in20,in01,in11,in21,in02,in12,in22;
755 double V00,V10,V20,V01,V11,V21,V02,V12,V22;
756 double ev0,ev1,ev2;
757 int s=inShape[0];
758 if (s==1) {
759 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
760 eigenvalues_and_eigenvectors1(in00,&ev0,&V00,tol);
761 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
762 V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
763 } else if (s==2) {
764 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
765 in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
766 in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
767 in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
768 eigenvalues_and_eigenvectors2(in00,(in01+in10)/2.,in11,
769 &ev0,&ev1,&V00,&V10,&V01,&V11,tol);
770 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
771 ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
772 V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
773 V[inOffset+DataTypes::getRelIndex(VShape,1,0)]=V10;
774 V[inOffset+DataTypes::getRelIndex(VShape,0,1)]=V01;
775 V[inOffset+DataTypes::getRelIndex(VShape,1,1)]=V11;
776 } else if (s==3) {
777 in00=in[inOffset+DataTypes::getRelIndex(inShape,0,0)];
778 in10=in[inOffset+DataTypes::getRelIndex(inShape,1,0)];
779 in20=in[inOffset+DataTypes::getRelIndex(inShape,2,0)];
780 in01=in[inOffset+DataTypes::getRelIndex(inShape,0,1)];
781 in11=in[inOffset+DataTypes::getRelIndex(inShape,1,1)];
782 in21=in[inOffset+DataTypes::getRelIndex(inShape,2,1)];
783 in02=in[inOffset+DataTypes::getRelIndex(inShape,0,2)];
784 in12=in[inOffset+DataTypes::getRelIndex(inShape,1,2)];
785 in22=in[inOffset+DataTypes::getRelIndex(inShape,2,2)];
786 eigenvalues_and_eigenvectors3(in00,(in01+in10)/2.,(in02+in20)/2.,in11,(in21+in12)/2.,in22,
787 &ev0,&ev1,&ev2,
788 &V00,&V10,&V20,&V01,&V11,&V21,&V02,&V12,&V22,tol);
789 ev[evOffset+DataTypes::getRelIndex(evShape,0)]=ev0;
790 ev[evOffset+DataTypes::getRelIndex(evShape,1)]=ev1;
791 ev[evOffset+DataTypes::getRelIndex(evShape,2)]=ev2;
792 V[inOffset+DataTypes::getRelIndex(VShape,0,0)]=V00;
793 V[inOffset+DataTypes::getRelIndex(VShape,1,0)]=V10;
794 V[inOffset+DataTypes::getRelIndex(VShape,2,0)]=V20;
795 V[inOffset+DataTypes::getRelIndex(VShape,0,1)]=V01;
796 V[inOffset+DataTypes::getRelIndex(VShape,1,1)]=V11;
797 V[inOffset+DataTypes::getRelIndex(VShape,2,1)]=V21;
798 V[inOffset+DataTypes::getRelIndex(VShape,0,2)]=V02;
799 V[inOffset+DataTypes::getRelIndex(VShape,1,2)]=V12;
800 V[inOffset+DataTypes::getRelIndex(VShape,2,2)]=V22;
801
802 }
803 }
804
805
806 /**
807 Inline function definitions.
808 */
809
810 template <class VEC>
811 inline
812 bool
813 checkOffset(const VEC& data,
814 const DataTypes::ShapeType& shape,
815 typename VEC::size_type offset)
816 {
817 return (data.size() >= (offset+DataTypes::noValues(shape)));
818 }
819
820 /**
821 * This assumes that all data involved have the same points per sample and same shape
822 */
823 template <class ResVEC, class LVEC, class RSCALAR>
824 void
825 binaryOpVectorRightScalar(ResVEC& res, // where result is to be stored
826 typename ResVEC::size_type resOffset, // offset in the result vector to start storing results
827 const typename ResVEC::size_type samplesToProcess, // number of samples to be updated in the result
828 const typename ResVEC::size_type sampleSize, // number of values in each sample
829 const LVEC& left, // LHS of calculation
830 typename LVEC::size_type leftOffset, // where to start reading LHS values
831 const RSCALAR* right, // RHS of the calculation
832 const bool rightreset, // true if RHS is providing a single sample of 1 value only
833 escript::ES_optype operation, // operation to perform
834 bool singleleftsample) // set to false for normal operation
835 {
836 size_t substep=(rightreset?0:1);
837 switch (operation)
838 {
839 case ADD:
840 {
841 #pragma omp parallel for
842 for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
843 {
844 typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
845 const RSCALAR* rpos=right+(rightreset?0:i*substep);
846
847 for (typename ResVEC::size_type j=0;j<sampleSize;++j)
848 {
849 res[i*sampleSize+resOffset+j]=left[leftbase+j]+*rpos;
850 }
851 }
852 }
853 break;
854 case POW:
855 {
856 #pragma omp parallel for
857 for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
858 {
859 typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
860 const RSCALAR* rpos=right+(rightreset?0:i*substep);
861
862 for (typename ResVEC::size_type j=0;j<sampleSize;++j)
863 {
864 res[i*sampleSize+resOffset+j]=pow(left[leftbase+j],*rpos);
865 }
866 }
867 }
868 break;
869 case SUB:
870 {
871 #pragma omp parallel for
872 for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
873 {
874 typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
875 const RSCALAR* rpos=right+(rightreset?0:i*substep);
876
877 for (typename ResVEC::size_type j=0;j<sampleSize;++j)
878 {
879 res[i*sampleSize+resOffset+j]=left[leftbase+j]-*rpos;
880 }
881 }
882 }
883 break;
884 case MUL:
885 {
886 #pragma omp parallel for
887 for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
888 {
889 typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
890 const RSCALAR* rpos=right+(rightreset?0:i*substep);
891
892 for (typename ResVEC::size_type j=0;j<sampleSize;++j)
893 {
894 res[i*sampleSize+resOffset+j]=left[leftbase+j] * *rpos;
895 }
896 }
897 }
898 break;
899 case DIV:
900 {
901 #pragma omp parallel for
902 for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
903 {
904 typename LVEC::size_type leftbase=leftOffset+(singleleftsample?0:i*sampleSize);
905 const RSCALAR* rpos=right+(rightreset?0:i*substep);
906
907 for (typename ResVEC::size_type j=0;j<sampleSize;++j)
908 {
909 res[i*sampleSize+resOffset+j]=left[leftbase+j]/ *rpos;
910 }
911 }
912 }
913 break;
914 default:
915 throw DataException("Unsupported binary operation");
916 }
917 }
918
919 template<>
920 void
921 binaryOpVectorRightScalar(DataTypes::RealVectorType& res, // where result is to be stored
922 typename DataTypes::RealVectorType::size_type resOffset, // offset in the result vector to start storing results
923 const typename DataTypes::RealVectorType::size_type samplesToProcess, // number of samples to be updated in the result
924 const typename DataTypes::RealVectorType::size_type sampleSize, // number of values in each sample
925 const DataTypes::RealVectorType& left, // LHS of calculation
926 typename DataTypes::RealVectorType::size_type leftOffset, // where to start reading LHS values
927 const DataTypes::real_t* right, // RHS of the calculation
928 const bool rightreset, // true if RHS is providing a single sample of 1 value only
929 escript::ES_optype operation, // operation to perform
930 bool singleleftsample);
931
932 /**
933 * This assumes that all data involved have the same points per sample and same shape
934 */
935 template <class ResVEC, class LSCALAR, class RVEC>
936 void
937 binaryOpVectorLeftScalar(ResVEC& res, // where result is to be stored
938 typename ResVEC::size_type resOffset, // offset in the result vector to start storing results
939 const typename ResVEC::size_type samplesToProcess, // number of samples to be updated in the result
940 const typename ResVEC::size_type sampleSize, // number of values in each sample
941 const LSCALAR* left, // LHS of calculation
942 const bool leftreset, // true if LHS is providing a single sample of 1 value only
943 const RVEC& right, // RHS of the calculation
944 typename RVEC::size_type rightOffset, // where to start reading RHS values
945 escript::ES_optype operation, // operation to perform
946 bool singlerightsample) // right consists of a single sample
947 {
948 size_t substep=(leftreset?0:1);
949 switch (operation)
950 {
951 case ADD:
952 {
953 #pragma omp parallel for
954 for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
955 {
956 typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
957 const LSCALAR* lpos=left+(leftreset?0:i*substep);
958 for (typename ResVEC::size_type j=0;j<sampleSize;++j)
959 {
960 res[i*sampleSize+resOffset+j]=*lpos+right[rightbase+j];
961 }
962 }
963 }
964 break;
965 case POW:
966 {
967 #pragma omp parallel for
968 for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
969 {
970 typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
971 const LSCALAR* lpos=left+(leftreset?0:i*substep);
972 for (typename ResVEC::size_type j=0;j<sampleSize;++j)
973 {
974 res[i*sampleSize+resOffset+j]=pow(*lpos,right[rightbase+j]);
975 }
976 }
977 }
978 break;
979 case SUB:
980 {
981 #pragma omp parallel for
982 for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
983 {
984 typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
985 const LSCALAR* lpos=left+(leftreset?0:i*substep);
986 for (typename ResVEC::size_type j=0;j<sampleSize;++j)
987 {
988 res[i*sampleSize+resOffset+j]=*lpos-right[rightbase+j];
989 }
990 }
991 }
992 break;
993 case MUL:
994 {
995 #pragma omp parallel for
996 for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
997 {
998 typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
999 const LSCALAR* lpos=left+(leftreset?0:i*substep);
1000 for (typename ResVEC::size_type j=0;j<sampleSize;++j)
1001 {
1002 res[i*sampleSize+resOffset+j]=*lpos*right[rightbase+j];
1003 }
1004 }
1005 }
1006 break;
1007 case DIV:
1008 {
1009 #pragma omp parallel for
1010 for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1011 {
1012 typename RVEC::size_type rightbase=rightOffset+(singlerightsample?0:i*sampleSize);
1013 const LSCALAR* lpos=left+(leftreset?0:i*substep);
1014 for (typename ResVEC::size_type j=0;j<sampleSize;++j)
1015 {
1016 res[i*sampleSize+resOffset+j]=*lpos/right[rightbase+j];
1017 }
1018 }
1019 }
1020 break;
1021 default:
1022 throw DataException("Unsupported binary operation");
1023 }
1024 }
1025
1026 template <>
1027 void
1028 binaryOpVectorLeftScalar(DataTypes::RealVectorType& res, // where result is to be stored
1029 typename DataTypes::RealVectorType::size_type resOffset, // offset in the result vector to start storing results
1030 const typename DataTypes::RealVectorType::size_type samplesToProcess, // number of samples to be updated in the result
1031 const typename DataTypes::RealVectorType::size_type sampleSize, // number of values in each sample
1032 const DataTypes::real_t* left, // LHS of calculation
1033 const bool leftreset, // true if LHS is providing a single sample of 1 value only
1034 const DataTypes::RealVectorType& right, // RHS of the calculation
1035 typename DataTypes::RealVectorType::size_type rightOffset, // where to start reading RHS values
1036 escript::ES_optype operation, // operation to perform
1037 bool singlerightsample); // right consists of a single sample
1038
1039 /**
1040 * This assumes that all data involved have the same points per sample and same shape
1041 */
1042 template <class ResVEC, class LVEC, class RVEC>
1043 void
1044 binaryOpVector(ResVEC& res, // where result is to be stored
1045 typename ResVEC::size_type resOffset, // offset in the result vector to start storing results
1046 const typename ResVEC::size_type samplesToProcess, // number of samples to be updated in the result
1047 const typename ResVEC::size_type sampleSize, // number of values in each sample
1048 const LVEC& left, // LHS of calculation
1049 typename LVEC::size_type leftOffset, // where to start reading LHS values
1050 const bool leftreset, // Is LHS only supplying a single sample instead of a bunch of them
1051 const RVEC& right, // RHS of the calculation
1052 typename RVEC::size_type rightOffset, // where to start reading RHS values
1053 const bool rightreset, // Is RHS only supplying a single sample instead of a bunch of them
1054 escript::ES_optype operation) // operation to perform
1055 {
1056 switch (operation)
1057 {
1058 case ADD:
1059 {
1060 #pragma omp parallel for
1061 for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1062 {
1063 typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
1064 typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
1065 for (typename ResVEC::size_type j=0;j<sampleSize;++j)
1066 {
1067 res[i*sampleSize+resOffset+j]=left[leftbase+j]+right[rightbase+j];
1068 }
1069 }
1070 }
1071 break;
1072 case POW:
1073 {
1074 #pragma omp parallel for
1075 for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1076 {
1077 typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
1078 typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
1079 for (typename ResVEC::size_type j=0;j<sampleSize;++j)
1080 {
1081 res[i*sampleSize+resOffset+j]=pow(left[leftbase+j],right[rightbase+j]);
1082 }
1083 }
1084 }
1085 break;
1086 case SUB:
1087 {
1088 #pragma omp parallel for
1089 for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1090 {
1091 typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
1092 typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
1093 for (typename ResVEC::size_type j=0;j<sampleSize;++j)
1094 {
1095 res[i*sampleSize+resOffset+j]=left[leftbase+j]-right[rightbase+j];
1096 }
1097 }
1098 }
1099 break;
1100 case MUL:
1101 {
1102 #pragma omp parallel for
1103 for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1104 {
1105 typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
1106 typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
1107 for (typename ResVEC::size_type j=0;j<sampleSize;++j)
1108 {
1109 res[i*sampleSize+resOffset+j]=left[leftbase+j]*right[rightbase+j];
1110 }
1111 }
1112 }
1113 break;
1114 case DIV:
1115 {
1116 #pragma omp parallel for
1117 for (typename ResVEC::size_type i=0;i<samplesToProcess;++i)
1118 {
1119 typename LVEC::size_type leftbase=leftOffset+(leftreset?0:i*sampleSize);
1120 typename RVEC::size_type rightbase=rightOffset+(rightreset?0:i*sampleSize);
1121 for (typename ResVEC::size_type j=0;j<sampleSize;++j)
1122 {
1123 res[i*sampleSize+resOffset+j]=left[leftbase+j]/right[rightbase+j];
1124 }
1125 }
1126 }
1127 break;
1128 default:
1129 throw DataException("Unsupported binary operation");
1130 }
1131 }
1132
1133 template <>
1134 void
1135 binaryOpVector(DataTypes::RealVectorType& res, // where result is to be stored
1136 typename DataTypes::RealVectorType::size_type resOffset, // offset in the result vector to start storing results
1137 const typename DataTypes::RealVectorType::size_type samplesToProcess, // number of samples to be updated in the result
1138 const typename DataTypes::RealVectorType::size_type sampleSize, // number of values in each sample
1139 const DataTypes::RealVectorType& left, // LHS of calculation
1140 typename DataTypes::RealVectorType::size_type leftOffset, // where to start reading LHS values
1141 const bool leftreset, // Is LHS only supplying a single sample instead of a bunch of them
1142 const DataTypes::RealVectorType& right, // RHS of the calculation
1143 typename DataTypes::RealVectorType::size_type rightOffset, // where to start reading RHS values
1144 const bool rightreset, // Is RHS only supplying a single sample instead of a bunch of them
1145 escript::ES_optype operation); // operation to perform
1146
1147 #define OPVECLAZYBODY(X) \
1148 for (size_t j=0;j<onumsteps;++j)\
1149 {\
1150 for (size_t i=0;i<numsteps;++i,res+=resultStep) \
1151 { \
1152 for (size_t s=0; s<chunksize; ++s)\
1153 {\
1154 res[s] = X;\
1155 }\
1156 /* tensor_binary_operation< TYPE >(chunksize, &((*left)[lroffset]), &((*right)[rroffset]), resultp, X);*/ \
1157 lroffset+=leftstep; \
1158 rroffset+=rightstep; \
1159 }\
1160 lroffset+=oleftstep;\
1161 rroffset+=orightstep;\
1162 }
1163
1164 /**
1165 * This assumes that all data involved have the same points per sample and same shape
1166 * This version is to be called from within DataLazy.
1167 * It does not have openmp around loops because it will be evaluating individual samples
1168 * (Which will be done within an enclosing openmp region.
1169 */
1170 template <class ResELT, class LELT, class RELT>
1171 void
1172 binaryOpVectorLazyArithmeticHelper(ResELT* res,
1173 const LELT* left,
1174 const RELT* right,
1175 const size_t chunksize,
1176 const size_t onumsteps,
1177 const size_t numsteps,
1178 const size_t resultStep,
1179 const size_t leftstep,
1180 const size_t rightstep,
1181 const size_t oleftstep,
1182 const size_t orightstep,
1183 size_t lroffset,
1184 size_t rroffset,
1185 escript::ES_optype operation) // operation to perform
1186 {
1187 switch (operation)
1188 {
1189 case ADD:
1190 OPVECLAZYBODY((left[lroffset+s]+right[rroffset+s]));
1191 break;
1192 case POW:
1193 OPVECLAZYBODY(pow(left[lroffset+s],right[rroffset+s]))
1194 break;
1195 case SUB:
1196 OPVECLAZYBODY(left[lroffset+s]-right[rroffset+s])
1197 break;
1198 case MUL:
1199 OPVECLAZYBODY(left[lroffset+s]*right[rroffset+s])
1200 break;
1201 case DIV:
1202 OPVECLAZYBODY(left[lroffset+s]/right[rroffset+s])
1203 break;
1204 default:
1205 ESYS_ASSERT(false, "Invalid operation. This should never happen!");
1206 // I can't throw here because this will be called inside a parallel section
1207 }
1208 }
1209
1210
1211 /**
1212 * This assumes that all data involved have the same points per sample and same shape
1213 * This version is to be called from within DataLazy.
1214 * It does not have openmp around loops because it will be evaluating individual samples
1215 * (Which will be done within an enclosing openmp region.
1216 */
1217 template <class ResELT, class LELT, class RELT>
1218 void
1219 binaryOpVectorLazyRelationalHelper(ResELT* res,
1220 const LELT* left,
1221 const RELT* right,
1222 const size_t chunksize,
1223 const size_t onumsteps,
1224 const size_t numsteps,
1225 const size_t resultStep,
1226 const size_t leftstep,
1227 const size_t rightstep,
1228 const size_t oleftstep,
1229 const size_t orightstep,
1230 size_t lroffset,
1231 size_t rroffset,
1232 escript::ES_optype operation) // operation to perform
1233 {
1234 switch (operation)
1235 {
1236 case LESS:
1237 OPVECLAZYBODY(left[lroffset+s]<right[rroffset+s])
1238 break;
1239 case GREATER:
1240 OPVECLAZYBODY(left[lroffset+s]>right[rroffset+s])
1241 break;
1242 case GREATER_EQUAL:
1243 OPVECLAZYBODY(left[lroffset+s]>=right[rroffset+s])
1244 break;
1245 case LESS_EQUAL:
1246 OPVECLAZYBODY(left[lroffset+s]<=right[rroffset+s])
1247 break;
1248 default:
1249 ESYS_ASSERT(false, "Invalid operation. This should never happen!");
1250 // I can't throw here because this will be called inside a parallel section
1251 }
1252 }
1253
1254
1255 /**
1256 * This assumes that all data involved have the same points per sample and same shape
1257 */
1258 /* trying to make a single version for all Tagged+Expanded interactions */
1259 template <class ResVEC, class LVEC, class RVEC>
1260 void
1261 binaryOpVectorTagged(ResVEC& res, // where result is to be stored
1262 const typename ResVEC::size_type samplesToProcess, // number of samples to be updated in the result
1263 const typename ResVEC::size_type DPPSample, // number of datapoints per sample
1264 const typename ResVEC::size_type DPSize, // datapoint size
1265 const LVEC& left, // LHS of calculation
1266 bool leftscalar,
1267 const RVEC& right, // RHS of the calculation
1268 bool rightscalar,
1269 bool lefttagged, // true if left object is the tagged one
1270 const DataTagged& tagsource, // where to get tag offsets from
1271 escript::ES_optype operation) // operation to perform
1272 {
1273 typename ResVEC::size_type lstep=leftscalar?1:DPSize;
1274 typename ResVEC::size_type rstep=rightscalar?1:DPSize;
1275 typename ResVEC::size_type limit=samplesToProcess*DPPSample;
1276 switch (operation)
1277 {
1278 case ADD:
1279 {
1280 #pragma omp parallel for
1281 for (typename ResVEC::size_type i=0;i<limit;++i)
1282 {
1283 typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
1284 typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
1285
1286 for (typename ResVEC::size_type j=0;j<DPSize;++j)
1287 {
1288 res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]+right[rightbase+j*(!rightscalar)];
1289 }
1290 }
1291 }
1292 break;
1293 case POW:
1294 {
1295 #pragma omp parallel for
1296 for (typename ResVEC::size_type i=0;i<limit;++i)
1297 {
1298 typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
1299 typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
1300
1301 for (typename ResVEC::size_type j=0;j<DPSize;++j)
1302 {
1303 res[i*DPSize+j]=pow(left[leftbase+j*(!leftscalar)],right[rightbase+j*(!rightscalar)]);
1304 }
1305 }
1306 }
1307 break;
1308 case SUB:
1309 {
1310 #pragma omp parallel for
1311 for (typename ResVEC::size_type i=0;i<limit;++i)
1312 {
1313 typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
1314 typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
1315
1316 for (typename ResVEC::size_type j=0;j<DPSize;++j)
1317 {
1318 res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]-right[rightbase+j*(!rightscalar)];
1319 }
1320 }
1321 }
1322 break;
1323 case MUL:
1324 {
1325 #pragma omp parallel for
1326 for (typename ResVEC::size_type i=0;i<limit;++i)
1327 {
1328 typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
1329 typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
1330
1331 for (typename ResVEC::size_type j=0;j<DPSize;++j)
1332 {
1333 res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]*right[rightbase+j*(!rightscalar)];
1334 }
1335 }
1336 }
1337 break;
1338 case DIV:
1339 {
1340 #pragma omp parallel for
1341 for (typename ResVEC::size_type i=0;i<limit;++i)
1342 {
1343 typename LVEC::size_type leftbase=(lefttagged?tagsource.getPointOffset(i/DPPSample,0):i*lstep); // only one of these
1344 typename RVEC::size_type rightbase=(lefttagged?i*rstep:tagsource.getPointOffset(i/DPPSample,0)); // will apply
1345
1346 for (typename ResVEC::size_type j=0;j<DPSize;++j)
1347 {
1348 res[i*DPSize+j]=left[leftbase+j*(!leftscalar)]/right[rightbase+j*(!rightscalar)];
1349 }
1350 }
1351 }
1352 break;
1353 default:
1354 throw DataException("Unsupported binary operation");
1355 }
1356 }
1357
1358 template<>
1359 void
1360 binaryOpVectorTagged(DataTypes::RealVectorType& res, // where result is to be stored
1361 const typename DataTypes::RealVectorType::size_type samplesToProcess, // number of samples to be updated in the result
1362 const typename DataTypes::RealVectorType::size_type DPPSample, // number of datapoints per sample
1363 const typename DataTypes::RealVectorType::size_type DPSize, // datapoint size
1364 const DataTypes::RealVectorType& left, // LHS of calculation
1365 const bool leftscalar,
1366 const DataTypes::RealVectorType& right, // RHS of the calculation
1367 const bool rightscalar,
1368 const bool lefttagged, // true if left object is the tagged one
1369 const DataTagged& tagsource, // where to get tag offsets from
1370 escript::ES_optype operation);
1371
1372
1373
1374
1375 /**
1376 \brief
1377 Perform the given data point reduction operation on the data point
1378 specified by the given offset into the view. Reduces all elements of
1379 the data point using the given operation, returning the result as a
1380 scalar. Operation must be a pointer to a function.
1381
1382 Called by escript::algorithm.
1383
1384 \param left - vector containing the datapoint
1385 \param shape - shape of datapoints in the vector
1386 \param offset - beginning of datapoint in the vector
1387 \param operation - Input -
1388 Operation to apply. Must be a pointer to a function.
1389 \param initial_value
1390 */
1391 template <class BinaryFunction>
1392 inline
1393 DataTypes::real_t
1394 reductionOpVector(const DataTypes::RealVectorType& left,
1395 const DataTypes::ShapeType& leftShape,
1396 DataTypes::RealVectorType::size_type offset,
1397 BinaryFunction operation,
1398 DataTypes::real_t initial_value)
1399 {
1400 ESYS_ASSERT((left.size()>0)&&checkOffset(left,leftShape,offset),
1401 "Couldn't perform reductionOp due to insufficient storage.");
1402 DataTypes::real_t current_value=initial_value;
1403 for (DataTypes::RealVectorType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
1404 current_value=operation(current_value,left[offset+i]);
1405 }
1406 return current_value;
1407 }
1408
1409 template <class BinaryFunction>
1410 inline
1411 DataTypes::real_t
1412 reductionOpVector(const DataTypes::CplxVectorType& left,
1413 const DataTypes::ShapeType& leftShape,
1414 DataTypes::CplxVectorType::size_type offset,
1415 BinaryFunction operation,
1416 DataTypes::real_t initial_value)
1417 {
1418 ESYS_ASSERT((left.size()>0)&&checkOffset(left,leftShape,offset),
1419 "Couldn't perform reductionOp due to insufficient storage.");
1420 DataTypes::real_t current_value=initial_value;
1421 for (DataTypes::RealVectorType::size_type i=0;i<DataTypes::noValues(leftShape);i++) {
1422 current_value=operation(current_value,left[offset+i]);
1423 }
1424 return current_value;
1425 }
1426
1427
1428 /**
1429 \brief
1430 computes the inverses of square (up to 3x3) matricies
1431
1432 \param in - vector containing the input matricies
1433 \param inShape - shape of the input matricies
1434 \param inOffset - the beginning of the input matricies within the vector "in"
1435 \param out - vector to store the inverses
1436 \param outShape - expected shape of the inverses
1437 \param outOffset - starting location for storing the inverses in out
1438 \param count - number of matricies to invert
1439 \param helper - associated working storage
1440
1441 \exception DataException if input and output are not the correct shape or if any of the matricies are not invertible.
1442 \return 0 on success, on failure the return value should be passed to matrixInverseError(int err).
1443 */
1444 int
1445 matrix_inverse(const DataTypes::RealVectorType& in,
1446 const DataTypes::ShapeType& inShape,
1447 DataTypes::RealVectorType::size_type inOffset,
1448 DataTypes::RealVectorType& out,
1449 const DataTypes::ShapeType& outShape,
1450 DataTypes::RealVectorType::size_type outOffset,
1451 int count,
1452 LapackInverseHelper& helper);
1453
1454 /**
1455 \brief
1456 throws an appropriate exception based on failure of matrix_inverse.
1457
1458 \param err - error code returned from matrix_inverse
1459 \warning do not call in a parallel region since it throws.
1460 */
1461 void
1462 matrixInverseError(int err);
1463
1464 /**
1465 \brief returns true if the vector contains NaN
1466
1467 */
1468 inline
1469 bool
1470 vectorHasNaN(const DataTypes::RealVectorType& in, DataTypes::RealVectorType::size_type inOffset, size_t count)
1471 {
1472 for (size_t z=inOffset;z<inOffset+count;++z)
1473 {
1474 if (nancheck(in[z]))
1475 {
1476 return true;
1477 }
1478 }
1479 return false;
1480 }
1481
1482 } // end namespace escript
1483
1484 #endif // __ESCRIPT_DATAMATHS_H__
1485

  ViewVC Help
Powered by ViewVC 1.1.26