Intrepid
Intrepid_RealSpaceToolsDef.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38// Denis Ridzal (dridzal@sandia.gov), or
39// Kara Peterson (kjpeter@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
50namespace Intrepid {
51
52
53
54template<class Scalar>
55void RealSpaceTools<Scalar>::absval(Scalar* absArray, const Scalar* inArray, const int size) {
56 for (size_t i=0; i<size; i++) {
57 absArray[i] = std::abs(inArray[i]);
58 }
59}
60
61
62
63template<class Scalar>
64void RealSpaceTools<Scalar>::absval(Scalar* inoutAbsArray, const int size) {
65 for (size_t i=0; i<size; i++) {
66 inoutAbsArray[i] = std::abs(inoutAbsArray[i]);
67 }
68}
69
70
71
72template<class Scalar>
73template<class ArrayAbs, class ArrayIn>
74void RealSpaceTools<Scalar>::absval(ArrayAbs & absArray, const ArrayIn & inArray) {
75#ifdef HAVE_INTREPID_DEBUG
76 TEUCHOS_TEST_FOR_EXCEPTION( ( getrank(inArray) != getrank(absArray) ),
77 std::invalid_argument,
78 ">>> ERROR (RealSpaceTools::absval): Array arguments must have identical ranks!");
79 for (size_t i=0; i<getrank(inArray); i++) {
80 TEUCHOS_TEST_FOR_EXCEPTION( ( static_cast<size_t>(inArray.dimension(i)) != static_cast<size_t>(absArray.dimension(i)) ),
81 std::invalid_argument,
82 ">>> ERROR (RealSpaceTools::absval): Dimensions of array arguments do not agree!");
83 }
84#endif
85
86 ArrayWrapper<Scalar,ArrayAbs, Rank<ArrayAbs >::value, false>absArrayWrap(absArray);
87 ArrayWrapper<Scalar,ArrayIn, Rank<ArrayIn >::value, true>inArrayWrap(inArray);
88
89 int inArrayRank=getrank(inArray);
90
91
92 if(inArrayRank==5){
93 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++)
94 for (size_t j=0; j<static_cast<size_t>(static_cast<size_t>(inArray.dimension(1))); j++)
95 for (size_t k=0; k<static_cast<size_t>(static_cast<size_t>(inArray.dimension(2))); k++)
96 for (size_t l=0; l<static_cast<size_t>(static_cast<size_t>(inArray.dimension(3))); l++)
97 for (size_t m=0; m<static_cast<size_t>(static_cast<size_t>(inArray.dimension(4))); m++){
98 absArrayWrap(i,j,k,l,m) = std::abs(inArrayWrap(i,j,k,l,m));
99 }
100 }else if(inArrayRank==4){
101 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++)
102 for (size_t j=0; j<static_cast<size_t>(static_cast<size_t>(inArray.dimension(1))); j++)
103 for (size_t k=0; k<static_cast<size_t>(static_cast<size_t>(inArray.dimension(2))); k++)
104 for (size_t l=0; l<static_cast<size_t>(static_cast<size_t>(inArray.dimension(3))); l++){
105 absArrayWrap(i,j,k,l) = std::abs(inArrayWrap(i,j,k,l));
106 }
107 }else if(inArrayRank==3){
108 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++)
109 for (size_t j=0; j<static_cast<size_t>(static_cast<size_t>(inArray.dimension(1))); j++)
110 for (size_t k=0; k<static_cast<size_t>(static_cast<size_t>(inArray.dimension(2))); k++){
111 absArrayWrap(i,j,k) = std::abs(inArrayWrap(i,j,k));
112 }
113 }else if(inArrayRank==2){
114 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++)
115 for (size_t j=0; j<static_cast<size_t>(static_cast<size_t>(inArray.dimension(1))); j++){
116 absArrayWrap(i,j) = std::abs(inArrayWrap(i,j));
117 }
118 }else if(inArrayRank==1){
119 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++){
120 absArrayWrap(i) = std::abs(inArrayWrap(i));
121
122 }
123 }
124}
125
126
127
128template<class Scalar>
129template<class ArrayInOut>
130void RealSpaceTools<Scalar>::absval(ArrayInOut & inoutAbsArray) {
131 for (size_t i=0; i<(size_t)inoutAbsArray.size(); i++) {
132 inoutAbsArray[i] = std::abs(inoutAbsArray[i]);
133 }
134}
135
136
137
138template<class Scalar>
139Scalar RealSpaceTools<Scalar>::vectorNorm(const Scalar* inVec, const size_t dim, const ENorm normType) {
140 Scalar temp = (Scalar)0;
141 switch(normType) {
142 case NORM_TWO:
143 for(size_t i = 0; i < dim; i++){
144 temp += inVec[i]*inVec[i];
145 }
146 temp = std::sqrt(temp);
147 break;
148 case NORM_INF:
149 temp = std::abs(inVec[0]);
150 for(size_t i = 1; i < dim; i++){
151 Scalar absData = std::abs(inVec[i]);
152 if (temp < absData) temp = absData;
153 }
154 break;
155 case NORM_ONE:
156 for(size_t i = 0; i < dim; i++){
157 temp += std::abs(inVec[i]);
158 }
159 break;
160 default:
161 TEUCHOS_TEST_FOR_EXCEPTION( ( (normType != NORM_TWO) && (normType != NORM_INF) && (normType != NORM_ONE) ),
162 std::invalid_argument,
163 ">>> ERROR (RealSpaceTools::vectorNorm): Invalid argument normType.");
164 }
165 return temp;
166}
167
168template<class Scalar>
169template<class ArrayIn>
170Scalar RealSpaceTools<Scalar>::vectorNorm(const ArrayIn & inVec, const ENorm normType) {
171#ifdef HAVE_INTREPID_DEBUG
172 TEUCHOS_TEST_FOR_EXCEPTION( ( !(getrank(inVec) >= 1 && getrank(inVec) <= 5) ),
173 std::invalid_argument,
174 ">>> ERROR (RealSpaceTools::vectorNorm): Vector argument must have rank 1!");
175#endif
176 ArrayWrapper<Scalar,ArrayIn, Rank<ArrayIn >::value, true>inVecWrap(inVec);
177 int inVecRank=getrank(inVecWrap);
178 Scalar temp = (Scalar)0;
179 switch(normType) {
180 case NORM_TWO:{
181 if(inVecRank==5){
182 for (size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
183 for (size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
184 for (size_t k=0; k<static_cast<size_t>(inVec.dimension(2)); k++)
185 for (size_t l=0; l<static_cast<size_t>(inVec.dimension(3)); l++)
186 for (size_t m=0; m<static_cast<size_t>(inVec.dimension(4)); m++)
187 temp += inVecWrap(i,j,k,l,m)*inVecWrap(i,j,k,l,m);
188 }else if(inVecRank==4){
189 for (size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
190 for (size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
191 for (size_t k=0; k<static_cast<size_t>(inVec.dimension(2)); k++)
192 for (size_t l=0; l<static_cast<size_t>(inVec.dimension(3)); l++)
193 temp += inVecWrap(i,j,k,l)*inVecWrap(i,j,k,l);
194 }else if(inVecRank==3){
195 for (size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
196 for (size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
197 for (size_t k=0; k<static_cast<size_t>(inVec.dimension(2)); k++)
198 temp += inVecWrap(i,j,k)*inVecWrap(i,j,k);
199 }else if(inVecRank==2){
200 for (size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
201 for (size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
202 temp += inVecWrap(i,j)*inVecWrap(i,j);
203 }else if(inVecRank==1){
204 for (size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
205 temp += inVecWrap(i)*inVecWrap(i);
206 }
207 temp = std::sqrt(temp);
208 }
209 break;
210 case NORM_INF:{
211
212 if(inVecRank==5){
213 temp = std::abs(inVecWrap(0,0,0,0,0));
214 for (size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
215 for (size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
216 for (size_t k=0; k<static_cast<size_t>(inVec.dimension(2)); k++)
217 for (size_t l=0; l<static_cast<size_t>(inVec.dimension(3)); l++)
218 for (size_t m=1; m<static_cast<size_t>(inVec.dimension(4)); m++){
219 Scalar absData = std::abs(inVecWrap(i,j,k,l,m));
220 if (temp < absData) temp = absData;
221 }
222 }else if(inVecRank==4){
223 temp = std::abs(inVecWrap(0,0,0,0));
224 for (size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
225 for (size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
226 for (size_t k=0; k<static_cast<size_t>(inVec.dimension(2)); k++)
227 for (size_t l=1; l<static_cast<size_t>(inVec.dimension(3)); l++){
228 Scalar absData = std::abs(inVecWrap(i,j,k,l));
229 if (temp < absData) temp = absData;
230 }
231 }else if(inVecRank==3){
232 temp = std::abs(inVecWrap(0,0,0));
233 for (size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
234 for (size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
235 for (size_t k=0; k<static_cast<size_t>(inVec.dimension(2)); k++){
236 Scalar absData = std::abs(inVecWrap(i,j,k));
237 if (temp < absData) temp = absData;
238 }
239 }else if(inVecRank==2){
240 temp = std::abs(inVecWrap(0,0));
241 for (size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
242 for (size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++){
243 Scalar absData = std::abs(inVecWrap(i,j));
244 if (temp < absData) temp = absData;
245 }
246 }else if(inVecRank==1){
247 temp = std::abs(inVecWrap(0));
248 for (size_t i=1; i<static_cast<size_t>(inVec.dimension(0)); i++){
249 Scalar absData = std::abs(inVecWrap(i));
250 if (temp < absData) temp = absData;
251 }
252 }
253}
254 break;
255 case NORM_ONE:{
256 if(inVecRank==5){
257 for (size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
258 for (size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
259 for (size_t k=0; k<static_cast<size_t>(inVec.dimension(2)); k++)
260 for (size_t l=0; l<static_cast<size_t>(inVec.dimension(3)); l++)
261 for (size_t m=0; m<static_cast<size_t>(inVec.dimension(4)); m++){
262 temp += std::abs(inVecWrap(i,j,k,l,m));
263 }
264 }else if(inVecRank==4){
265 for (size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
266 for (size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
267 for (size_t k=0; k<static_cast<size_t>(inVec.dimension(2)); k++)
268 for (size_t l=0; l<static_cast<size_t>(inVec.dimension(3)); l++){
269 temp += std::abs(inVecWrap(i,j,k,l));
270 }
271 }else if(inVecRank==3){
272 for (size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
273 for (size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++)
274 for (size_t k=0; k<static_cast<size_t>(inVec.dimension(2)); k++){
275 temp += std::abs(inVecWrap(i,j,k));
276 }
277 }else if(inVecRank==2){
278 for (size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++)
279 for (size_t j=0; j<static_cast<size_t>(inVec.dimension(1)); j++){
280 temp += std::abs(inVecWrap(i,j));
281 }
282 }else if(inVecRank==1){
283 for (size_t i=0; i<static_cast<size_t>(inVec.dimension(0)); i++){
284 temp += std::abs(inVecWrap(i));
285 }
286 }
287}
288 break;
289 default:
290 TEUCHOS_TEST_FOR_EXCEPTION( ( (normType != NORM_TWO) && (normType != NORM_INF) && (normType != NORM_ONE) ),
291 std::invalid_argument,
292 ">>> ERROR (RealSpaceTools::vectorNorm): Invalid argument normType.");
293 }
294 return temp;
295}
296/*
297template<class Scalar>
298template<class ArrayIn>
299Scalar RealSpaceTools<Scalar>::vectorNorm(const ArrayIn & inVec, const ENorm normType) {
300
301#ifdef HAVE_INTREPID_DEBUG
302 TEUCHOS_TEST_FOR_EXCEPTION( ( inVec.rank() != 1 ),
303 std::invalid_argument,
304 ">>> ERROR (RealSpaceTools::vectorNorm): Vector argument must have rank 1!");
305#endif
306
307 int size = inVec.size();
308
309 Scalar temp = (Scalar)0;
310 switch(normType) {
311 case NORM_TWO:
312 for(size_t i = 0; i < size; i++){
313 temp += inVec[i]*inVec[i];
314 }
315 temp = std::sqrt(temp);
316 break;
317 case NORM_INF:
318 temp = std::abs(inVec[0]);
319 for(size_t i = 1; i < size; i++){
320 Scalar absData = std::abs(inVec[i]);
321 if (temp < absData) temp = absData;
322 }
323 break;
324 case NORM_ONE:
325 for(size_t i = 0; i < size; i++){
326 temp += std::abs(inVec[i]);
327 }
328 break;
329 default:
330 TEUCHOS_TEST_FOR_EXCEPTION( ( (normType != NORM_TWO) && (normType != NORM_INF) && (normType != NORM_ONE) ),
331 std::invalid_argument,
332 ">>> ERROR (RealSpaceTools::vectorNorm): Invalid argument normType.");
333 }
334 return temp;
335}*/
336template<class Scalar>
337template<class ArrayNorm, class ArrayIn>
338void RealSpaceTools<Scalar>::vectorNorm(ArrayNorm & normArray, const ArrayIn & inVecs, const ENorm normType) {
339
340 ArrayWrapper<Scalar,ArrayNorm, Rank<ArrayNorm >::value, false>normArrayWrap(normArray);
341 ArrayWrapper<Scalar,ArrayIn, Rank<ArrayIn >::value, true>inVecsWrap(inVecs);
342
343 size_t arrayRank = getrank(inVecs);
344#ifdef HAVE_INTREPID_DEBUG
345 TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != getrank(normArray)+1 ),
346 std::invalid_argument,
347 ">>> ERROR (RealSpaceTools::vectorNorm): Ranks of norm and vector array arguments are incompatible!");
348 TEUCHOS_TEST_FOR_EXCEPTION( ( (arrayRank < 2) || (arrayRank > 3) ),
349 std::invalid_argument,
350 ">>> ERROR (RealSpaceTools::vectorNorm): Rank of vector array must be 2 or 3!");
351 for (size_t i=0; i<arrayRank-1; i++) {
352 TEUCHOS_TEST_FOR_EXCEPTION( ( inVecs.dimension(i) != normArray.dimension(i) ),
353 std::invalid_argument,
354 ">>> ERROR (RealSpaceTools::vectorNorm): Dimensions of norm and vector arguments do not agree!");
355 }
356#endif
357
358 size_t dim_i0 = 1; // first index dimension (e.g. cell index)
359 size_t dim_i1 = 1; // second index dimension (e.g. point index)
360 size_t dim = static_cast<size_t>(inVecs.dimension(arrayRank-1)); // spatial dimension
361
362 // determine i0 and i1 dimensions
363 switch(arrayRank) {
364 case 3:
365 dim_i0 = static_cast<size_t>(inVecs.dimension(0));
366 dim_i1 = static_cast<size_t>(inVecs.dimension(1));
367 switch(normType) {
368 case NORM_TWO: {
369 for (size_t i0=0; i0<dim_i0; i0++) {
370 for (size_t i1=0; i1<dim_i1; i1++) {
371 Scalar temp = (Scalar)0;
372 for(size_t i = 0; i < dim; i++){
373 temp += inVecsWrap(i0,i1,i)*inVecsWrap(i0,i1,i);
374 }
375 normArrayWrap(i0,i1) = std::sqrt(temp);
376 }
377 }
378 break;
379 } // case NORM_TWO
380
381 case NORM_INF: {
382 for (size_t i0=0; i0<dim_i0; i0++) {
383 for (size_t i1=0; i1<dim_i1; i1++) {
384 Scalar temp = (Scalar)0;
385 temp = std::abs(inVecsWrap(i0,i1,0));
386 for(size_t i = 1; i < dim; i++){
387 Scalar absData = std::abs(inVecsWrap(i0,i1,i));
388 if (temp < absData) temp = absData;
389 }
390 normArrayWrap(i0,i1) = temp;
391 }
392 }
393 break;
394 } // case NORM_INF
395
396 case NORM_ONE: {
397 for (size_t i0=0; i0<dim_i0; i0++) {
398 for (size_t i1=0; i1<dim_i1; i1++) {
399 Scalar temp = (Scalar)0;
400 for(size_t i = 0; i < dim; i++){
401 temp += std::abs(inVecsWrap(i0,i1,i));
402 }
403 normArrayWrap(i0,i1) = temp;
404 }
405 }
406 break;
407 } // case NORM_ONE
408
409 default:
410 TEUCHOS_TEST_FOR_EXCEPTION( ( (normType != NORM_TWO) && (normType != NORM_INF) && (normType != NORM_ONE) ),
411 std::invalid_argument,
412 ">>> ERROR (RealSpaceTools::vectorNorm): Invalid argument normType.");
413 }
414
415
416
417 break;
418 case 2:
419 dim_i1 = static_cast<size_t>(inVecs.dimension(0));
420 switch(normType) {
421 case NORM_TWO: {
422
423 for (size_t i1=0; i1<dim_i1; i1++) {
424 Scalar temp = (Scalar)0;
425 for(size_t i = 0; i < dim; i++){
426 temp += inVecsWrap(i1,i)*inVecsWrap(i1,i);
427 }
428 normArrayWrap(i1) = std::sqrt(temp);
429 }
430
431 break;
432 } // case NORM_TWO
433
434 case NORM_INF: {
435 for (size_t i1=0; i1<dim_i1; i1++) {
436 Scalar temp = (Scalar)0;
437 temp = std::abs(inVecsWrap(i1,0));
438 for(size_t i = 1; i < dim; i++){
439 Scalar absData = std::abs(inVecsWrap(i1,i));
440 if (temp < absData) temp = absData;
441 }
442 normArrayWrap(i1) = temp;
443 }
444 break;
445 } // case NORM_INF
446
447 case NORM_ONE: {
448 for (size_t i1=0; i1<dim_i1; i1++) {
449 Scalar temp = (Scalar)0;
450 for(size_t i = 0; i < dim; i++){
451 temp += std::abs(inVecsWrap(i1,i));
452 }
453 normArrayWrap(i1) = temp;
454 }
455 break;
456 } // case NORM_ONE
457
458 default:
459 TEUCHOS_TEST_FOR_EXCEPTION( ( (normType != NORM_TWO) && (normType != NORM_INF) && (normType != NORM_ONE) ),
460 std::invalid_argument,
461 ">>> ERROR (RealSpaceTools::vectorNorm): Invalid argument normType.");
462 }
463
464
465
466 break;
467 }
468
469
470}
471/*
472
473template<class Scalar>
474template<class ArrayNorm, class ArrayIn>
475void RealSpaceTools<Scalar>::vectorNorm(ArrayNorm & normArray, const ArrayIn & inVecs, const ENorm normType) {
476
477 int arrayRank = inVecs.rank();
478
479#ifdef HAVE_INTREPID_DEBUG
480 TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != normArray.rank()+1 ),
481 std::invalid_argument,
482 ">>> ERROR (RealSpaceTools::vectorNorm): Ranks of norm and vector array arguments are incompatible!");
483 TEUCHOS_TEST_FOR_EXCEPTION( ( (arrayRank < 2) || (arrayRank > 3) ),
484 std::invalid_argument,
485 ">>> ERROR (RealSpaceTools::vectorNorm): Rank of vector array must be 2 or 3!");
486 for (int i=0; i<arrayRank-1; i++) {
487 TEUCHOS_TEST_FOR_EXCEPTION( ( inVecs.dimension(i) != normArray.dimension(i) ),
488 std::invalid_argument,
489 ">>> ERROR (RealSpaceTools::vectorNorm): Dimensions of norm and vector arguments do not agree!");
490 }
491#endif
492
493 size_t dim_i0 = 1; // first index dimension (e.g. cell index)
494 size_t dim_i1 = 1; // second index dimension (e.g. point index)
495 size_t dim = inVecs.dimension(arrayRank-1); // spatial dimension
496
497 // determine i0 and i1 dimensions
498 switch(arrayRank) {
499 case 3:
500 dim_i0 = inVecs.dimension(0);
501 dim_i1 = inVecs.dimension(1);
502 break;
503 case 2:
504 dim_i1 = inVecs.dimension(0);
505 break;
506 }
507
508 switch(normType) {
509 case NORM_TWO: {
510 int offset_i0, offset, normOffset;
511 for (size_t i0=0; i0<dim_i0; i0++) {
512 offset_i0 = i0*dim_i1;
513 for (size_t i1=0; i1<dim_i1; i1++) {
514 offset = offset_i0 + i1;
515 normOffset = offset;
516 offset *= dim;
517 Scalar temp = (Scalar)0;
518 for(size_t i = 0; i < dim; i++){
519 temp += inVecs[offset+i]*inVecs[offset+i];
520 }
521 normArray[normOffset] = std::sqrt(temp);
522 }
523 }
524 break;
525 } // case NORM_TWO
526
527 case NORM_INF: {
528 int offset_i0, offset, normOffset;
529 for (size_t i0=0; i0<dim_i0; i0++) {
530 offset_i0 = i0*dim_i1;
531 for (size_t i1=0; i1<dim_i1; i1++) {
532 offset = offset_i0 + i1;
533 normOffset = offset;
534 offset *= dim;
535 Scalar temp = (Scalar)0;
536 temp = std::abs(inVecs[offset]);
537 for(size_t i = 1; i < dim; i++){
538 Scalar absData = std::abs(inVecs[offset+i]);
539 if (temp < absData) temp = absData;
540 }
541 normArray[normOffset] = temp;
542 }
543 }
544 break;
545 } // case NORM_INF
546
547 case NORM_ONE: {
548 int offset_i0, offset, normOffset;
549 for (size_t i0=0; i0<dim_i0; i0++) {
550 offset_i0 = i0*dim_i1;
551 for (size_t i1=0; i1<dim_i1; i1++) {
552 offset = offset_i0 + i1;
553 normOffset = offset;
554 offset *= dim;
555 Scalar temp = (Scalar)0;
556 for(size_t i = 0; i < dim; i++){
557 temp += std::abs(inVecs[offset+i]);
558 }
559 normArray[normOffset] = temp;
560 }
561 }
562 break;
563 } // case NORM_ONE
564
565 default:
566 TEUCHOS_TEST_FOR_EXCEPTION( ( (normType != NORM_TWO) && (normType != NORM_INF) && (normType != NORM_ONE) ),
567 std::invalid_argument,
568 ">>> ERROR (RealSpaceTools::vectorNorm): Invalid argument normType.");
569 }
570}
571
572
573*/
574
575template<class Scalar>
576void RealSpaceTools<Scalar>::transpose(Scalar* transposeMat, const Scalar* inMat, const size_t dim) {
577 for(size_t i=0; i < dim; i++){
578 transposeMat[i*dim+i]=inMat[i*dim+i]; // Set diagonal elements
579 for(size_t j=i+1; j < dim; j++){
580 transposeMat[i*dim+j]=inMat[j*dim+i]; // Set off-diagonal elements
581 transposeMat[j*dim+i]=inMat[i*dim+j];
582 }
583 }
584}
585
586
587
588template<class Scalar>
589template<class ArrayTranspose, class ArrayIn>
590void RealSpaceTools<Scalar>::transpose(ArrayTranspose & transposeMats, const ArrayIn & inMats) {
591 size_t arrayRank = getrank(inMats);
592#ifdef HAVE_INTREPID_DEBUG
593 TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != getrank(transposeMats) ),
594 std::invalid_argument,
595 ">>> ERROR (RealSpaceTools::transpose): Matrix array arguments do not have identical ranks!");
596 TEUCHOS_TEST_FOR_EXCEPTION( ( (arrayRank < 2) || (arrayRank > 4) ),
597 std::invalid_argument,
598 ">>> ERROR (RealSpaceTools::transpose): Rank of matrix array must be 2, 3, or 4!");
599 for (size_t i=0; i<arrayRank; i++) {
600 TEUCHOS_TEST_FOR_EXCEPTION( ( static_cast<size_t>(inMats.dimension(i)) != static_cast<size_t>(transposeMats.dimension(i)) ),
601 std::invalid_argument,
602 ">>> ERROR (RealSpaceTools::transpose): Dimensions of matrix arguments do not agree!");
603 }
604 TEUCHOS_TEST_FOR_EXCEPTION( ( static_cast<size_t>(inMats.dimension(arrayRank-2)) != static_cast<size_t>(inMats.dimension(arrayRank-1)) ),
605 std::invalid_argument,
606 ">>> ERROR (RealSpaceTools::transpose): Matrices are not square!");
607#endif
608 size_t dim_i0 = 1; // first index dimension (e.g. cell index)
609 size_t dim_i1 = 1; // second index dimension (e.g. point index)
610 size_t dim = static_cast<size_t>(inMats.dimension(arrayRank-2)); // spatial dimension
611
612
613 ArrayWrapper<Scalar,ArrayTranspose,Rank<ArrayTranspose>::value,false>transposeArr(transposeMats);
614 ArrayWrapper<Scalar,ArrayIn,Rank<ArrayIn>::value,true>inputArr(inMats);
615 // determine i0 and i1 dimensions
616 switch(arrayRank) {
617 case 4:
618 dim_i0 = static_cast<size_t>(inMats.dimension(0));
619 dim_i1 = static_cast<size_t>(inMats.dimension(1));
620
621 for (size_t i0=0; i0<dim_i0; i0++) {
622 for (size_t i1=0; i1<dim_i1; i1++) {
623 for(size_t i=0; i < dim; i++){
624 transposeArr(i0,i1,i,i)=inputArr(i0,i1,i,i);
625 for(size_t j=i+1; j < dim; j++){
626 transposeArr(i0,i1,i,j)=inputArr(i0,i1,j,i);
627 transposeArr(i0,i1,j,i)=inputArr(i0,i1,i,j);
628 }
629 }
630
631 } // i1
632 } // i0
633 break;
634 case 3:
635 dim_i1 = static_cast<size_t>(inMats.dimension(0));
636 for (size_t i1=0; i1<dim_i1; i1++) {
637 for(size_t i=0; i < dim; i++){
638 transposeArr(i1,i,i)=inputArr(i1,i,i);
639 for(size_t j=i+1; j < dim; j++){
640 transposeArr(i1,i,j)=inputArr(i1,j,i);
641 transposeArr(i1,j,i)=inputArr(i1,i,j);
642 }
643 } // i1
644 }
645 break;
646 }
647
648
649
650}
651
652template<class Scalar>
653void RealSpaceTools<Scalar>::inverse(Scalar* inverseMat, const Scalar* inMat, const size_t dim) {
654
655 switch(dim) {
656 case 3: {
657 size_t i, j, rowID = 0, colID = 0;
658 int rowperm[3]={0,1,2};
659 int colperm[3]={0,1,2}; // Complete pivoting
660 Scalar emax(0);
661
662 for(i=0; i < 3; ++i){
663 for(j=0; j < 3; ++j){
664 if( std::abs( inMat[i*3+j] ) > emax){
665 rowID = i; colID = j; emax = std::abs( inMat[i*3+j] );
666 }
667 }
668 }
669#ifdef HAVE_INTREPID_DEBUG
670 TEUCHOS_TEST_FOR_EXCEPTION( ( emax == (Scalar)0 ),
671 std::invalid_argument,
672 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
673#endif
674 if( rowID ){
675 rowperm[0] = rowID;
676 rowperm[rowID] = 0;
677 }
678 if( colID ){
679 colperm[0] = colID;
680 colperm[colID] = 0;
681 }
682 Scalar B[3][3], S[2][2], Bi[3][3]; // B=rowperm inMat colperm, S=Schur complement(Boo)
683 for(i=0; i < 3; ++i){
684 for(j=0; j < 3; ++j){
685 B[i][j] = inMat[rowperm[i]*3+colperm[j]];
686 }
687 }
688 B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot
689 for(i=0; i < 2; ++i){
690 for(j=0; j < 2; ++j){
691 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y'
692 }
693 }
694 Scalar detS = S[0][0]*S[1][1]- S[0][1]*S[1][0], Si[2][2];
695#ifdef HAVE_INTREPID_DEBUG
696 TEUCHOS_TEST_FOR_EXCEPTION( ( detS == (Scalar)0 ),
697 std::invalid_argument,
698 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
699#endif
700
701 Si[0][0] = S[1][1]/detS; Si[0][1] = -S[0][1]/detS;
702 Si[1][0] = -S[1][0]/detS; Si[1][1] = S[0][0]/detS;
703
704 for(j=0; j<2;j++)
705 Bi[0][j+1] = -( B[0][1]*Si[0][j] + B[0][2]* Si[1][j])/B[0][0];
706 for(i=0; i<2;i++)
707 Bi[i+1][0] = -(Si[i][0]*B[1][0] + Si[i][1]*B[2][0]);
708
709 Bi[0][0] = ((Scalar)1/B[0][0])-Bi[0][1]*B[1][0]-Bi[0][2]*B[2][0];
710 Bi[1][1] = Si[0][0];
711 Bi[1][2] = Si[0][1];
712 Bi[2][1] = Si[1][0];
713 Bi[2][2] = Si[1][1];
714 for(i=0; i < 3; ++i){
715 for(j=0; j < 3; ++j){
716 inverseMat[i*3+j] = Bi[colperm[i]][rowperm[j]]; // set inverse
717 }
718 }
719 break;
720 } // case 3
721
722 case 2: {
723
724 Scalar determinant = inMat[0]*inMat[3]-inMat[1]*inMat[2];;
725#ifdef HAVE_INTREPID_DEBUG
726 TEUCHOS_TEST_FOR_EXCEPTION( ( (inMat[0]==(Scalar)0) && (inMat[1]==(Scalar)0) &&
727 (inMat[2]==(Scalar)0) && (inMat[3]==(Scalar)0) ),
728 std::invalid_argument,
729 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
730 TEUCHOS_TEST_FOR_EXCEPTION( ( determinant == (Scalar)0 ),
731 std::invalid_argument,
732 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
733#endif
734 inverseMat[0] = inMat[3] / determinant;
735 inverseMat[1] = - inMat[1] / determinant;
736 //
737 inverseMat[2] = - inMat[2] / determinant;
738 inverseMat[3] = inMat[0] / determinant;
739 break;
740 } // case 2
741
742 case 1: {
743#ifdef HAVE_INTREPID_DEBUG
744 TEUCHOS_TEST_FOR_EXCEPTION( ( inMat[0] == (Scalar)0 ),
745 std::invalid_argument,
746 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
747#endif
748 inverseMat[0] = (Scalar)1 / inMat[0];
749 break;
750 } // case 1
751
752 } // switch (dim)
753}
754
755template<class Scalar>
756template<class ArrayInverse, class ArrayIn>
757void RealSpaceTools<Scalar>::inverse(ArrayInverse & inverseMats, const ArrayIn & inMats) {
758
759 ArrayWrapper<Scalar,ArrayInverse, Rank<ArrayInverse >::value, false>inverseMatsWrap(inverseMats);
760 ArrayWrapper<Scalar,ArrayIn, Rank<ArrayIn >::value, true>inMatsWrap(inMats);
761
762 size_t arrayRank = getrank(inMats);
763
764#ifdef HAVE_INTREPID_DEBUG
765 TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != getrank(inverseMats) ),
766 std::invalid_argument,
767 ">>> ERROR (RealSpaceTools::inverse): Matrix array arguments do not have identical ranks!");
768 TEUCHOS_TEST_FOR_EXCEPTION( ( (arrayRank < 2) || (arrayRank > 4) ),
769 std::invalid_argument,
770 ">>> ERROR (RealSpaceTools::inverse): Rank of matrix array must be 2, 3, or 4!");
771 for (size_t i=0; i<arrayRank; i++) {
772 TEUCHOS_TEST_FOR_EXCEPTION( ( static_cast<size_t>(inMats.dimension(i)) != static_cast<size_t>(inverseMats.dimension(i)) ),
773 std::invalid_argument,
774 ">>> ERROR (RealSpaceTools::inverse): Dimensions of matrix arguments do not agree!");
775 }
776 TEUCHOS_TEST_FOR_EXCEPTION( ( static_cast<size_t>(inMats.dimension(arrayRank-2)) != static_cast<size_t>(inMats.dimension(arrayRank-1)) ),
777 std::invalid_argument,
778 ">>> ERROR (RealSpaceTools::inverse): Matrices are not square!");
779 TEUCHOS_TEST_FOR_EXCEPTION( ( (static_cast<size_t>(inMats.dimension(arrayRank-2)) < 1) || (static_cast<size_t>(inMats.dimension(arrayRank-2)) > 3) ),
780 std::invalid_argument,
781 ">>> ERROR (RealSpaceTools::inverse): Spatial dimension must be 1, 2, or 3!");
782#endif
783
784 size_t dim_i0 = 1; // first index dimension (e.g. cell index)
785 size_t dim_i1 = 1; // second index dimension (e.g. point index)
786 size_t dim = static_cast<size_t>(inMats.dimension(arrayRank-2)); // spatial dimension
787
788 // determine i0 and i1 dimensions
789 switch(arrayRank) {
790 case 4:
791 dim_i0 = static_cast<size_t>(inMats.dimension(0));
792 dim_i1 = static_cast<size_t>(inMats.dimension(1));
793 switch(dim) {
794 case 3: {
795
796
797 for (size_t i0=0; i0<dim_i0; i0++) {
798
799 for (size_t i1=0; i1<dim_i1; i1++) {
800
801
802 size_t i, j, rowID = 0, colID = 0;
803 int rowperm[3]={0,1,2};
804 int colperm[3]={0,1,2}; // Complete pivoting
805 Scalar emax(0);
806
807 for(i=0; i < 3; ++i){
808 for(j=0; j < 3; ++j){
809 if( std::abs( inMatsWrap(i0,i1,i,j) ) > emax){
810 rowID = i; colID = j; emax = std::abs( inMatsWrap(i0,i1,i,j) );
811 }
812 }
813 }
814#ifdef HAVE_INTREPID_DEBUG
815#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
816 TEUCHOS_TEST_FOR_EXCEPTION( ( emax == (Scalar)0 ),
817 std::invalid_argument,
818 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
819#endif
820#endif
821 if( rowID ){
822 rowperm[0] = rowID;
823 rowperm[rowID] = 0;
824 }
825 if( colID ){
826 colperm[0] = colID;
827 colperm[colID] = 0;
828 }
829 Scalar B[3][3], S[2][2], Bi[3][3]; // B=rowperm inMat colperm, S=Schur complement(Boo)
830 for(i=0; i < 3; ++i){
831 for(j=0; j < 3; ++j){
832 B[i][j] = inMatsWrap(i0,i1,rowperm[i],colperm[j]);
833 }
834 }
835 B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot
836 for(i=0; i < 2; ++i){
837 for(j=0; j < 2; ++j){
838 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y'
839 }
840 }
841 Scalar detS = S[0][0]*S[1][1]- S[0][1]*S[1][0], Si[2][2];
842#ifdef HAVE_INTREPID_DEBUG
843#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
844 TEUCHOS_TEST_FOR_EXCEPTION( ( detS == (Scalar)0 ),
845 std::invalid_argument,
846 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
847#endif
848#endif
849
850 Si[0][0] = S[1][1]/detS; Si[0][1] = -S[0][1]/detS;
851 Si[1][0] = -S[1][0]/detS; Si[1][1] = S[0][0]/detS;
852
853 for(j=0; j<2;j++)
854 Bi[0][j+1] = -( B[0][1]*Si[0][j] + B[0][2]* Si[1][j])/B[0][0];
855 for(i=0; i<2;i++)
856 Bi[i+1][0] = -(Si[i][0]*B[1][0] + Si[i][1]*B[2][0]);
857
858 Bi[0][0] = ((Scalar)1/B[0][0])-Bi[0][1]*B[1][0]-Bi[0][2]*B[2][0];
859 Bi[1][1] = Si[0][0];
860 Bi[1][2] = Si[0][1];
861 Bi[2][1] = Si[1][0];
862 Bi[2][2] = Si[1][1];
863 for(i=0; i < 3; ++i){
864 for(j=0; j < 3; ++j){
865 inverseMatsWrap(i0,i1,i,j) = Bi[colperm[i]][rowperm[j]]; // set inverse
866 }
867 }
868 } // for i1
869 } // for i0
870 break;
871 } // case 3
872
873 case 2: {
874
875 for (size_t i0=0; i0<dim_i0; i0++) {
876
877 for (size_t i1=0; i1<dim_i1; i1++) {
878
879
880 Scalar determinant = inMatsWrap(i0,i1,0,0)*inMatsWrap(i0,i1,1,1)-inMatsWrap(i0,i1,0,1)*inMatsWrap(i0,i1,1,0);
881#ifdef HAVE_INTREPID_DEBUG
882#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
883 TEUCHOS_TEST_FOR_EXCEPTION( ( (inMatsWrap(i0,i1,0,0)==(Scalar)0) && (inMatsWrap(i0,i1,0,1)==(Scalar)0) &&
884 (inMatsWrap(i0,i1,1,0)==(Scalar)0) && (inMatsWrap(i0,i1,1,1)==(Scalar)0) ),
885 std::invalid_argument,
886 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
887 TEUCHOS_TEST_FOR_EXCEPTION( ( determinant == (Scalar)0 ),
888 std::invalid_argument,
889 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
890#endif
891#endif
892 inverseMatsWrap(i0,i1,0,0) = inMatsWrap(i0,i1,1,1) / determinant;
893 inverseMatsWrap(i0,i1,0,1) = - inMatsWrap(i0,i1,0,1) / determinant;
894 //
895 inverseMatsWrap(i0,i1,1,0) = - inMatsWrap(i0,i1,1,0) / determinant;
896 inverseMatsWrap(i0,i1,1,1) = inMatsWrap(i0,i1,0,0) / determinant;
897 } // for i1
898 } // for i0
899 break;
900 } // case 2
901
902 case 1: {
903 for (size_t i0=0; i0<dim_i0; i0++) {
904 for (size_t i1=0; i1<dim_i1; i1++) {
905#ifdef HAVE_INTREPID_DEBUG
906#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
907 TEUCHOS_TEST_FOR_EXCEPTION( ( inMatsWrap(i0,i1,0,0) == (Scalar)0 ),
908 std::invalid_argument,
909 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
910#endif
911#endif
912
913 inverseMatsWrap(i0,i1,0,0) = (Scalar)1 / inMatsWrap(i0,i1,0,0);
914 } // for i1
915 } // for i2
916
917
918 break;
919 } // case 1
920
921 } // switch (dim)
922 break;
923 case 3:
924 dim_i1 = static_cast<size_t>(inMats.dimension(0));
925 switch(dim) {
926 case 3: {
927
928 for (size_t i1=0; i1<dim_i1; i1++) {
929
930 size_t i, j, rowID = 0, colID = 0;
931 int rowperm[3]={0,1,2};
932 int colperm[3]={0,1,2}; // Complete pivoting
933 Scalar emax(0);
934
935 for(i=0; i < 3; ++i){
936 for(j=0; j < 3; ++j){
937 if( std::abs( inMatsWrap(i1,i,j) ) > emax){
938 rowID = i; colID = j; emax = std::abs( inMatsWrap(i1,i,j) );
939 }
940 }
941 }
942#ifdef HAVE_INTREPID_DEBUG
943#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
944 TEUCHOS_TEST_FOR_EXCEPTION( ( emax == (Scalar)0 ),
945 std::invalid_argument,
946 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
947#endif
948#endif
949
950 if( rowID ){
951 rowperm[0] = rowID;
952 rowperm[rowID] = 0;
953 }
954 if( colID ){
955 colperm[0] = colID;
956 colperm[colID] = 0;
957 }
958 Scalar B[3][3], S[2][2], Bi[3][3]; // B=rowperm inMat colperm, S=Schur complement(Boo)
959 for(i=0; i < 3; ++i){
960 for(j=0; j < 3; ++j){
961 B[i][j] = inMatsWrap(i1,rowperm[i],colperm[j]);
962 }
963 }
964 B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot
965 for(i=0; i < 2; ++i){
966 for(j=0; j < 2; ++j){
967 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y'
968 }
969 }
970 Scalar detS = S[0][0]*S[1][1]- S[0][1]*S[1][0], Si[2][2];
971#ifdef HAVE_INTREPID_DEBUG
972#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
973 TEUCHOS_TEST_FOR_EXCEPTION( ( detS == (Scalar)0 ),
974 std::invalid_argument,
975 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
976#endif
977#endif
978
979
980 Si[0][0] = S[1][1]/detS; Si[0][1] = -S[0][1]/detS;
981 Si[1][0] = -S[1][0]/detS; Si[1][1] = S[0][0]/detS;
982
983 for(j=0; j<2;j++)
984 Bi[0][j+1] = -( B[0][1]*Si[0][j] + B[0][2]* Si[1][j])/B[0][0];
985 for(i=0; i<2;i++)
986 Bi[i+1][0] = -(Si[i][0]*B[1][0] + Si[i][1]*B[2][0]);
987
988 Bi[0][0] = ((Scalar)1/B[0][0])-Bi[0][1]*B[1][0]-Bi[0][2]*B[2][0];
989 Bi[1][1] = Si[0][0];
990 Bi[1][2] = Si[0][1];
991 Bi[2][1] = Si[1][0];
992 Bi[2][2] = Si[1][1];
993 for(i=0; i < 3; ++i){
994 for(j=0; j < 3; ++j){
995 inverseMatsWrap(i1,i,j) = Bi[colperm[i]][rowperm[j]]; // set inverse
996 }
997 }
998 } // for i1
999
1000 break;
1001 } // case 3
1002
1003 case 2: {
1004
1005 for (size_t i1=0; i1<dim_i1; i1++) {
1006
1007
1008 Scalar determinant = inMatsWrap(i1,0,0)*inMatsWrap(i1,1,1)-inMatsWrap(i1,0,1)*inMatsWrap(i1,1,0);
1009#ifdef HAVE_INTREPID_DEBUG
1010#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
1011 TEUCHOS_TEST_FOR_EXCEPTION( ( (inMatsWrap(i1,0,0)==(Scalar)0) && (inMatsWrap(i1,0,1)==(Scalar)0) &&
1012 (inMatsWrap(i1,1,0)==(Scalar)0) && (inMatsWrap(i1,1,1)==(Scalar)0) ),
1013 std::invalid_argument,
1014 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
1015 TEUCHOS_TEST_FOR_EXCEPTION( ( determinant == (Scalar)0 ),
1016 std::invalid_argument,
1017 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
1018#endif
1019#endif
1020
1021 inverseMatsWrap(i1,0,0) = inMatsWrap(i1,1,1) / determinant;
1022 inverseMatsWrap(i1,0,1) = - inMatsWrap(i1,0,1) / determinant;
1023 //
1024 inverseMatsWrap(i1,1,0) = - inMatsWrap(i1,1,0) / determinant;
1025 inverseMatsWrap(i1,1,1) = inMatsWrap(i1,0,0) / determinant;
1026 } // for i1
1027
1028 break;
1029 } // case 2
1030
1031 case 1: {
1032
1033 for (size_t i1=0; i1<dim_i1; i1++) {
1034#ifdef HAVE_INTREPID_DEBUG
1035#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
1036 TEUCHOS_TEST_FOR_EXCEPTION( ( inMatsWrap(i1,0,0) == (Scalar)0 ),
1037 std::invalid_argument,
1038 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
1039#endif
1040#endif
1041
1042
1043 inverseMatsWrap(i1,0,0) = (Scalar)1 / inMatsWrap(i1,0,0);
1044 }
1045
1046
1047
1048 break;
1049 } // case 1
1050
1051 } // switch (dim)
1052 break;
1053 case 2:
1054 switch(dim) {
1055 case 3: {
1056 size_t i, j, rowID = 0, colID = 0;
1057 int rowperm[3]={0,1,2};
1058 int colperm[3]={0,1,2}; // Complete pivoting
1059 Scalar emax(0);
1060
1061 for(i=0; i < 3; ++i){
1062 for(j=0; j < 3; ++j){
1063 if( std::abs( inMatsWrap(i,j) ) > emax){
1064 rowID = i; colID = j; emax = std::abs( inMatsWrap(i,j) );
1065 }
1066 }
1067 }
1068#ifdef HAVE_INTREPID_DEBUG
1069#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
1070 TEUCHOS_TEST_FOR_EXCEPTION( ( emax == (Scalar)0 ),
1071 std::invalid_argument,
1072 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
1073#endif
1074#endif
1075
1076 if( rowID ){
1077 rowperm[0] = rowID;
1078 rowperm[rowID] = 0;
1079 }
1080 if( colID ){
1081 colperm[0] = colID;
1082 colperm[colID] = 0;
1083 }
1084 Scalar B[3][3], S[2][2], Bi[3][3]; // B=rowperm inMat colperm, S=Schur complement(Boo)
1085 for(i=0; i < 3; ++i){
1086 for(j=0; j < 3; ++j){
1087 B[i][j] = inMatsWrap(rowperm[i],colperm[j]);
1088 }
1089 }
1090 B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot
1091 for(i=0; i < 2; ++i){
1092 for(j=0; j < 2; ++j){
1093 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y'
1094 }
1095 }
1096 Scalar detS = S[0][0]*S[1][1]- S[0][1]*S[1][0], Si[2][2];
1097#ifdef HAVE_INTREPID_DEBUG
1098#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
1099 TEUCHOS_TEST_FOR_EXCEPTION( ( detS == (Scalar)0 ),
1100 std::invalid_argument,
1101 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
1102#endif
1103#endif
1104
1105 Si[0][0] = S[1][1]/detS; Si[0][1] = -S[0][1]/detS;
1106 Si[1][0] = -S[1][0]/detS; Si[1][1] = S[0][0]/detS;
1107
1108 for(j=0; j<2;j++)
1109 Bi[0][j+1] = -( B[0][1]*Si[0][j] + B[0][2]* Si[1][j])/B[0][0];
1110 for(i=0; i<2;i++)
1111 Bi[i+1][0] = -(Si[i][0]*B[1][0] + Si[i][1]*B[2][0]);
1112
1113 Bi[0][0] = ((Scalar)1/B[0][0])-Bi[0][1]*B[1][0]-Bi[0][2]*B[2][0];
1114 Bi[1][1] = Si[0][0];
1115 Bi[1][2] = Si[0][1];
1116 Bi[2][1] = Si[1][0];
1117 Bi[2][2] = Si[1][1];
1118 for(i=0; i < 3; ++i){
1119 for(j=0; j < 3; ++j){
1120 inverseMatsWrap(i,j) = Bi[colperm[i]][rowperm[j]]; // set inverse
1121 }
1122 }
1123
1124 break;
1125 } // case 3
1126
1127 case 2: {
1128 Scalar determinant = inMatsWrap(0,0)*inMatsWrap(1,1)-inMatsWrap(0,1)*inMatsWrap(1,0);
1129#ifdef HAVE_INTREPID_DEBUG
1130#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
1131 TEUCHOS_TEST_FOR_EXCEPTION( ( (inMatsWrap(0,0)==(Scalar)0) && (inMatsWrap(0,1)==(Scalar)0) &&
1132 (inMatsWrap(1,0)==(Scalar)0) && (inMatsWrap(1,1)==(Scalar)0) ),
1133 std::invalid_argument,
1134 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
1135 TEUCHOS_TEST_FOR_EXCEPTION( ( determinant == (Scalar)0 ),
1136 std::invalid_argument,
1137 ">>> ERROR (Matrix): Inverse of a singular matrix is undefined!");
1138#endif
1139#endif
1140 inverseMatsWrap(0,0) = inMatsWrap(1,1) / determinant;
1141 inverseMatsWrap(0,1) = - inMatsWrap(0,1) / determinant;
1142 //
1143 inverseMatsWrap(1,0) = - inMatsWrap(1,0) / determinant;
1144 inverseMatsWrap(1,1) = inMatsWrap(0,0) / determinant;
1145
1146 break;
1147 } // case 2
1148
1149 case 1: {
1150#ifdef HAVE_INTREPID_DEBUG
1151#ifdef HAVE_INTREPID_DEBUG_INF_CHECK
1152 TEUCHOS_TEST_FOR_EXCEPTION( ( inMatsWrap(0,0) == (Scalar)0 ),
1153 std::invalid_argument,
1154 ">>> ERROR (Matrix): Inverse of a zero matrix is undefined!");
1155#endif
1156#endif
1157 inverseMatsWrap(0,0) = (Scalar)1 / inMatsWrap(0,0);
1158 break;
1159 } // case 1
1160
1161 }
1162 break;
1163 }
1164
1165
1166}
1167
1168
1169
1170
1171template<class Scalar>
1172Scalar RealSpaceTools<Scalar>::det(const Scalar* inMat, const size_t dim) {
1173 Scalar determinant(0);
1174
1175 switch (dim) {
1176 case 3: {
1177 int i,j,rowID = 0;
1178 int colID = 0;
1179 int rowperm[3]={0,1,2};
1180 int colperm[3]={0,1,2}; // Complete pivoting
1181 Scalar emax(0);
1182
1183 for(i=0; i < 3; ++i){
1184 for(j=0; j < 3; ++j){
1185 if( std::abs( inMat[i*dim+j] ) > emax){
1186 rowID = i; colID = j; emax = std::abs( inMat[i*dim+j] );
1187 }
1188 }
1189 }
1190 if( emax > 0 ){
1191 if( rowID ){
1192 rowperm[0] = rowID;
1193 rowperm[rowID] = 0;
1194 }
1195 if( colID ){
1196 colperm[0] = colID;
1197 colperm[colID] = 0;
1198 }
1199 Scalar B[3][3], S[2][2]; // B=rowperm inMat colperm, S=Schur complement(Boo)
1200 for(i=0; i < 3; ++i){
1201 for(j=0; j < 3; ++j){
1202 B[i][j] = inMat[rowperm[i]*dim+colperm[j]];
1203 }
1204 }
1205 B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot
1206 for(i=0; i < 2; ++i){
1207 for(j=0; j < 2; ++j){
1208 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y'
1209 }
1210 }
1211 determinant = B[0][0] * (S[0][0] * S[1][1] - S[0][1] * S[1][0]); // det(B)
1212 if( rowID ) determinant = -determinant;
1213 if( colID ) determinant = -determinant;
1214 }
1215 break;
1216 } // case 3
1217
1218 case 2:
1219 determinant = inMat[0]*inMat[3]-
1220 inMat[1]*inMat[2];
1221 break;
1222
1223 case 1:
1224 determinant = inMat[0];
1225 break;
1226
1227 default:
1228 TEUCHOS_TEST_FOR_EXCEPTION( ( (dim != 1) && (dim != 2) && (dim != 3) ),
1229 std::invalid_argument,
1230 ">>> ERROR (Matrix): Invalid matrix dimension.");
1231 } // switch (dim)
1232
1233 return determinant;
1234}
1235
1236
1237
1238template<class Scalar>
1239template<class ArrayIn>
1240Scalar RealSpaceTools<Scalar>::det(const ArrayIn & inMat) {
1241
1242#ifdef HAVE_INTREPID_DEBUG
1243 TEUCHOS_TEST_FOR_EXCEPTION( (getrank(inMat) != 2),
1244 std::invalid_argument,
1245 ">>> ERROR (RealSpaceTools::det): Rank of matrix argument must be 2!");
1246 TEUCHOS_TEST_FOR_EXCEPTION( ( static_cast<size_t>(inMat.dimension(0)) != static_cast<size_t>(inMat.dimension(1)) ),
1247 std::invalid_argument,
1248 ">>> ERROR (RealSpaceTools::det): Matrix is not square!");
1249 TEUCHOS_TEST_FOR_EXCEPTION( ( (static_cast<size_t>(inMat.dimension(0)) < 1) || (static_cast<size_t>(inMat.dimension(0)) > 3) ),
1250 std::invalid_argument,
1251 ">>> ERROR (RealSpaceTools::det): Spatial dimension must be 1, 2, or 3!");
1252#endif
1253 ArrayWrapper<Scalar,ArrayIn, Rank<ArrayIn >::value, true>inMatWrap(inMat);
1254
1255 size_t dim = static_cast<size_t>(inMat.dimension(0));
1256 Scalar determinant(0);
1257
1258 switch (dim) {
1259 case 3: {
1260 int i,j,rowID = 0;
1261 int colID = 0;
1262 int rowperm[3]={0,1,2};
1263 int colperm[3]={0,1,2}; // Complete pivoting
1264 Scalar emax(0);
1265
1266 for(i=0; i < 3; ++i){
1267 for(j=0; j < 3; ++j){
1268 if( std::abs( inMatWrap(i,j) ) > emax){
1269 rowID = i; colID = j; emax = std::abs( inMatWrap(i,j) );
1270 }
1271 }
1272 }
1273 if( emax > 0 ){
1274 if( rowID ){
1275 rowperm[0] = rowID;
1276 rowperm[rowID] = 0;
1277 }
1278 if( colID ){
1279 colperm[0] = colID;
1280 colperm[colID] = 0;
1281 }
1282 Scalar B[3][3], S[2][2]; // B=rowperm inMat colperm, S=Schur complement(Boo)
1283 for(i=0; i < 3; ++i){
1284 for(j=0; j < 3; ++j){
1285 B[i][j] = inMatWrap(rowperm[i],colperm[j]);
1286 }
1287 }
1288 B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot
1289 for(i=0; i < 2; ++i){
1290 for(j=0; j < 2; ++j){
1291 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y'
1292 }
1293 }
1294 determinant = B[0][0] * (S[0][0] * S[1][1] - S[0][1] * S[1][0]); // det(B)
1295 if( rowID ) determinant = -determinant;
1296 if( colID ) determinant = -determinant;
1297 }
1298 break;
1299 } // case 3
1300
1301 case 2:
1302 determinant = inMatWrap(0,0)*inMatWrap(1,1)-
1303 inMatWrap(0,1)*inMatWrap(1,0);
1304 break;
1305
1306 case 1:
1307 determinant = inMatWrap(0,0);
1308 break;
1309
1310 default:
1311 TEUCHOS_TEST_FOR_EXCEPTION( ( (dim != 1) && (dim != 2) && (dim != 3) ),
1312 std::invalid_argument,
1313 ">>> ERROR (Matrix): Invalid matrix dimension.");
1314 } // switch (dim)
1315
1316 return determinant;
1317}
1318
1319template<class Scalar>
1320template<class ArrayDet, class ArrayIn>
1321void RealSpaceTools<Scalar>::det(ArrayDet & detArray, const ArrayIn & inMats) {
1322 ArrayWrapper<Scalar,ArrayDet, Rank<ArrayDet >::value, false>detArrayWrap(detArray);
1323 ArrayWrapper<Scalar,ArrayIn, Rank<ArrayIn >::value, true>inMatsWrap(inMats);
1324
1325 size_t matArrayRank = getrank(inMats);
1326#ifdef HAVE_INTREPID_DEBUG
1327 TEUCHOS_TEST_FOR_EXCEPTION( ( matArrayRank != getrank(detArray)+2 ),
1328 std::invalid_argument,
1329 ">>> ERROR (RealSpaceTools::det): Determinant and matrix array arguments do not have compatible ranks!");
1330 TEUCHOS_TEST_FOR_EXCEPTION( ( (matArrayRank < 3) || (matArrayRank > 4) ),
1331 std::invalid_argument,
1332 ">>> ERROR (RealSpaceTools::det): Rank of matrix array must be 3 or 4!");
1333 for (size_t i=0; i<matArrayRank-2; i++) {
1334 TEUCHOS_TEST_FOR_EXCEPTION( ( static_cast<size_t>(inMats.dimension(i)) != static_cast<size_t>(detArray.dimension(i)) ),
1335 std::invalid_argument,
1336 ">>> ERROR (RealSpaceTools::det): Dimensions of determinant and matrix array arguments do not agree!");
1337 }
1338 TEUCHOS_TEST_FOR_EXCEPTION( ( static_cast<size_t>(inMats.dimension(matArrayRank-2)) != static_cast<size_t>(inMats.dimension(matArrayRank-1)) ),
1339 std::invalid_argument,
1340 ">>> ERROR (RealSpaceTools::det): Matrices are not square!");
1341 TEUCHOS_TEST_FOR_EXCEPTION( ( (static_cast<size_t>(inMats.dimension(matArrayRank-2)) < 1) || (static_cast<size_t>(inMats.dimension(matArrayRank-2)) > 3) ),
1342 std::invalid_argument,
1343 ">>> ERROR (RealSpaceTools::det): Spatial dimension must be 1, 2, or 3!");
1344#endif
1345
1346 size_t dim_i0 = 1; // first index dimension (e.g. cell index)
1347 size_t dim_i1 = 1; // second index dimension (e.g. point index)
1348 size_t dim = inMats.dimension(matArrayRank-2); // spatial dimension
1349
1350 // determine i0 and i1 dimensions
1351 switch(matArrayRank) {
1352 case 4:
1353 dim_i0 = static_cast<size_t>(inMats.dimension(0));
1354 dim_i1 = static_cast<size_t>(inMats.dimension(1));
1355 switch(dim) {
1356 case 3: {
1357
1358
1359 for (size_t i0=0; i0<dim_i0; i0++) {
1360
1361 for (size_t i1=0; i1<dim_i1; i1++) {
1362
1363
1364 int i,j,rowID = 0;
1365 int colID = 0;
1366 int rowperm[3]={0,1,2};
1367 int colperm[3]={0,1,2}; // Complete pivoting
1368 Scalar emax(0), determinant(0);
1369
1370 for(i=0; i < 3; ++i){
1371 for(j=0; j < 3; ++j){
1372 if( std::abs( inMatsWrap(i0,i1,i,j) ) > emax){
1373 rowID = i; colID = j; emax = std::abs( inMatsWrap(i0,i1,i,j) );
1374 }
1375 }
1376 }
1377 if( emax > 0 ){
1378 if( rowID ){
1379 rowperm[0] = rowID;
1380 rowperm[rowID] = 0;
1381 }
1382 if( colID ){
1383 colperm[0] = colID;
1384 colperm[colID] = 0;
1385 }
1386 Scalar B[3][3], S[2][2]; // B=rowperm inMat colperm, S=Schur complement(Boo)
1387 for(i=0; i < 3; ++i){
1388 for(j=0; j < 3; ++j){
1389 B[i][j] = inMatsWrap(i0,i1,rowperm[i],colperm[j]);
1390 }
1391 }
1392 B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot
1393 for(i=0; i < 2; ++i){
1394 for(j=0; j < 2; ++j){
1395 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y'
1396 }
1397 }
1398 determinant = B[0][0] * (S[0][0] * S[1][1] - S[0][1] * S[1][0]); // det(B)
1399 if( rowID ) determinant = -determinant;
1400 if( colID ) determinant = -determinant;
1401 }
1402 detArrayWrap(i0,i1)= determinant;
1403 } // for i1
1404 } // for i0
1405 break;
1406 } // case 3
1407
1408 case 2: {
1409
1410 for (size_t i0=0; i0<dim_i0; i0++) {
1411
1412 for (size_t i1=0; i1<dim_i1; i1++) {
1413
1414 detArrayWrap(i0,i1) = inMatsWrap(i0,i1,0,0)*inMatsWrap(i0,i1,1,1)-inMatsWrap(i0,i1,0,1)*inMatsWrap(i0,i1,1,0);
1415 } // for i1
1416 } // for i0
1417 break;
1418 } // case 2
1419
1420 case 1: {
1421
1422 for (size_t i0=0; i0<dim_i0; i0++) {
1423
1424 for (size_t i1=0; i1<dim_i1; i1++) {
1425 detArrayWrap(i0,i1) = inMatsWrap(i0,i1,0,0);
1426 } // for i1
1427 } // for i2
1428 break;
1429 } // case 1
1430
1431 } // switch (dim)
1432 break;
1433 case 3:
1434 dim_i1 = static_cast<size_t>(inMats.dimension(0));
1435 switch(dim) {
1436 case 3: {
1437 for (size_t i1=0; i1<dim_i1; i1++) {
1438
1439 int i,j,rowID = 0;
1440 int colID = 0;
1441 int rowperm[3]={0,1,2};
1442 int colperm[3]={0,1,2}; // Complete pivoting
1443 Scalar emax(0), determinant(0);
1444
1445
1446 for(i=0; i < 3; ++i){
1447 for(j=0; j < 3; ++j){
1448 if( std::abs( inMatsWrap(i1,i,j) ) > emax){
1449 rowID = i; colID = j; emax = std::abs( inMatsWrap(i1,i,j) );
1450 }
1451 }
1452 }
1453 if( emax > 0 ){
1454 if( rowID ){
1455 rowperm[0] = rowID;
1456 rowperm[rowID] = 0;
1457 }
1458 if( colID ){
1459 colperm[0] = colID;
1460 colperm[colID] = 0;
1461 }
1462 Scalar B[3][3], S[2][2]; // B=rowperm inMat colperm, S=Schur complement(Boo)
1463 for(i=0; i < 3; ++i){
1464 for(j=0; j < 3; ++j){
1465 B[i][j] = inMatsWrap(i1,rowperm[i],colperm[j]);
1466 }
1467 }
1468 B[1][0] /= B[0][0]; B[2][0] /= B[0][0];// B(:,0)/=pivot
1469 for(i=0; i < 2; ++i){
1470 for(j=0; j < 2; ++j){
1471 S[i][j] = B[i+1][j+1] - B[i+1][0] * B[0][j+1]; // S = B -z*y'
1472 }
1473 }
1474 determinant = B[0][0] * (S[0][0] * S[1][1] - S[0][1] * S[1][0]); // det(B)
1475 if( rowID ) determinant = -determinant;
1476 if( colID ) determinant = -determinant;
1477 }
1478 detArrayWrap(i1) = determinant;
1479 }
1480 break;
1481 } // case 3
1482
1483 case 2: {
1484 for (size_t i1=0; i1<dim_i1; i1++) {
1485 detArrayWrap(i1) = inMatsWrap(i1,0,0)*inMatsWrap(i1,1,1)-inMatsWrap(i1,0,1)*inMatsWrap(i1,1,0);
1486 }
1487 break;
1488 } // case 2
1489
1490 case 1: {
1491 for (size_t i1=0; i1<dim_i1; i1++) {
1492 detArrayWrap(i1) = inMatsWrap(i1,0,0);
1493 }
1494 break;
1495 } // case 1
1496}
1497break;
1498
1499 } // switch (dim)
1500
1501 }
1502
1503
1504
1505template<class Scalar>
1506void RealSpaceTools<Scalar>::add(Scalar* sumArray, const Scalar* inArray1, const Scalar* inArray2, const int size) {
1507 for (size_t i=0; i<size; i++) {
1508 sumArray[i] = inArray1[i] + inArray2[i];
1509 }
1510}
1511
1512
1513
1514template<class Scalar>
1515void RealSpaceTools<Scalar>::add(Scalar* inoutSumArray, const Scalar* inArray, const int size) {
1516 for (size_t i=0; i<size; i++) {
1517 inoutSumArray[i] += inArray[i];
1518 }
1519}
1520
1521
1522
1523template<class Scalar>
1524template<class ArraySum, class ArrayIn1, class ArrayIn2>
1525void RealSpaceTools<Scalar>::add(ArraySum & sumArray, const ArrayIn1 & inArray1, const ArrayIn2 & inArray2) {
1526#ifdef HAVE_INTREPID_DEBUG
1527 TEUCHOS_TEST_FOR_EXCEPTION( ( (getrank(inArray1) != getrank(inArray2)) || (getrank(inArray1) != getrank(sumArray)) ),
1528 std::invalid_argument,
1529 ">>> ERROR (RealSpaceTools::add): Array arguments must have identical ranks!");
1530 for (size_t i=0; i<getrank(inArray1); i++) {
1531 TEUCHOS_TEST_FOR_EXCEPTION( ( (static_cast<size_t>(inArray1.dimension(i)) != static_cast<size_t>(inArray2.dimension(i))) || (static_cast<size_t>(inArray1.dimension(i)) != static_cast<size_t>(sumArray.dimension(i))) ),
1532 std::invalid_argument,
1533 ">>> ERROR (RealSpaceTools::add): Dimensions of array arguments do not agree!");
1534 }
1535#endif
1536
1537
1538
1539
1540 ArrayWrapper<Scalar,ArraySum, Rank<ArraySum >::value, false>sumArrayWrap(sumArray);
1541 ArrayWrapper<Scalar,ArrayIn1, Rank<ArrayIn1 >::value, true>inArray1Wrap(inArray1);
1542 ArrayWrapper<Scalar,ArrayIn2, Rank<ArrayIn2 >::value, true>inArray2Wrap(inArray2);
1543 int inArrayRank=getrank(inArray1);
1544
1545
1546 if(inArrayRank==5){
1547 for (size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++)
1548 for (size_t j=0; j<static_cast<size_t>(inArray1.dimension(1)); j++)
1549 for (size_t k=0; k<static_cast<size_t>(inArray1.dimension(2)); k++)
1550 for (size_t l=0; l<static_cast<size_t>(inArray1.dimension(3)); l++)
1551 for (size_t m=0; m<static_cast<size_t>(inArray1.dimension(4)); m++){
1552 sumArrayWrap(i,j,k,l,m) = inArray1Wrap(i,j,k,l,m)+inArray2Wrap(i,j,k,l,m);
1553 }
1554 }else if(inArrayRank==4){
1555 for (size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++)
1556 for (size_t j=0; j<static_cast<size_t>(inArray1.dimension(1)); j++)
1557 for (size_t k=0; k<static_cast<size_t>(inArray1.dimension(2)); k++)
1558 for (size_t l=0; l<static_cast<size_t>(inArray1.dimension(3)); l++){
1559 sumArrayWrap(i,j,k,l) = inArray1Wrap(i,j,k,l)+inArray2Wrap(i,j,k,l);
1560 }
1561 }else if(inArrayRank==3){
1562 for (size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++)
1563 for (size_t j=0; j<static_cast<size_t>(inArray1.dimension(1)); j++)
1564 for (size_t k=0; k<static_cast<size_t>(inArray1.dimension(2)); k++){
1565 sumArrayWrap(i,j,k) = inArray1Wrap(i,j,k)+inArray2Wrap(i,j,k);
1566 }
1567 }else if(inArrayRank==2){
1568 for (size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++)
1569 for (size_t j=0; j<static_cast<size_t>(inArray1.dimension(1)); j++){
1570 sumArrayWrap(i,j) = inArray1Wrap(i,j)+inArray2Wrap(i,j);
1571 }
1572 }else if(inArrayRank==1){
1573 for (size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++){
1574 sumArrayWrap(i) = inArray1Wrap(i)+inArray2Wrap(i);
1575
1576 }
1577 }
1578}
1579
1580
1581
1582template<class Scalar>
1583template<class ArraySum, class ArrayIn>
1584void RealSpaceTools<Scalar>::add(ArraySum & inoutSumArray, const ArrayIn & inArray) {
1585#ifdef HAVE_INTREPID_DEBUG
1586 TEUCHOS_TEST_FOR_EXCEPTION( ( getrank(inArray) != getrank(inoutSumArray) ),
1587 std::invalid_argument,
1588 ">>> ERROR (RealSpaceTools::add): Array arguments must have identical ranks!");
1589 for (size_t i=0; i<getrank(inArray); i++) {
1590 TEUCHOS_TEST_FOR_EXCEPTION( ( static_cast<size_t>(inArray.dimension(i)) != static_cast<size_t>(inoutSumArray.dimension(i)) ),
1591 std::invalid_argument,
1592 ">>> ERROR (RealSpaceTools::add): Dimensions of array arguments do not agree!");
1593 }
1594#endif
1595
1596 ArrayWrapper<Scalar,ArraySum, Rank<ArraySum >::value, false>inoutSumArrayWrap(inoutSumArray);
1597 ArrayWrapper<Scalar,ArrayIn, Rank<ArrayIn >::value, true>inArrayWrap(inArray);
1598 int inArrayRank=getrank(inArray);
1599
1600
1601 if(inArrayRank==5){
1602 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++)
1603 for (size_t j=0; j<static_cast<size_t>(static_cast<size_t>(inArray.dimension(1))); j++)
1604 for (size_t k=0; k<static_cast<size_t>(static_cast<size_t>(inArray.dimension(2))); k++)
1605 for (size_t l=0; l<static_cast<size_t>(static_cast<size_t>(inArray.dimension(3))); l++)
1606 for (size_t m=0; m<static_cast<size_t>(static_cast<size_t>(inArray.dimension(4))); m++){
1607 inoutSumArrayWrap(i,j,k,l,m) += inArrayWrap(i,j,k,l,m);
1608 }
1609 }else if(inArrayRank==4){
1610 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++)
1611 for (size_t j=0; j<static_cast<size_t>(static_cast<size_t>(inArray.dimension(1))); j++)
1612 for (size_t k=0; k<static_cast<size_t>(static_cast<size_t>(inArray.dimension(2))); k++)
1613 for (size_t l=0; l<static_cast<size_t>(static_cast<size_t>(inArray.dimension(3))); l++){
1614 inoutSumArrayWrap(i,j,k,l) += inArrayWrap(i,j,k,l);
1615 }
1616 }else if(inArrayRank==3){
1617 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++)
1618 for (size_t j=0; j<static_cast<size_t>(static_cast<size_t>(inArray.dimension(1))); j++)
1619 for (size_t k=0; k<static_cast<size_t>(inArray.dimension(2)); k++){
1620 inoutSumArrayWrap(i,j,k) += inArrayWrap(i,j,k);
1621 }
1622 }else if(inArrayRank==2){
1623 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++)
1624 for (size_t j=0; j<static_cast<size_t>(static_cast<size_t>(inArray.dimension(1))); j++){
1625 inoutSumArrayWrap(i,j) += inArrayWrap(i,j);
1626 }
1627 }else if(inArrayRank==1){
1628 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++){
1629 inoutSumArrayWrap(i) += inArrayWrap(i);
1630
1631 }
1632 }
1633}
1634
1635
1636
1637template<class Scalar>
1638void RealSpaceTools<Scalar>::subtract(Scalar* diffArray, const Scalar* inArray1, const Scalar* inArray2, const int size) {
1639 for (size_t i=0; i<size; i++) {
1640 diffArray[i] = inArray1[i] - inArray2[i];
1641 }
1642}
1643
1644
1645
1646template<class Scalar>
1647void RealSpaceTools<Scalar>::subtract(Scalar* inoutDiffArray, const Scalar* inArray, const int size) {
1648 for (int i=0; i<size; i++) {
1649 inoutDiffArray[i] -= inArray[i];
1650 }
1651}
1652
1653
1654
1655template<class Scalar>
1656template<class ArrayDiff, class ArrayIn1, class ArrayIn2>
1657void RealSpaceTools<Scalar>::subtract(ArrayDiff & diffArray, const ArrayIn1 & inArray1, const ArrayIn2 & inArray2) {
1658 ArrayWrapper<Scalar,ArrayDiff, Rank<ArrayDiff >::value, false>diffArrayWrap(diffArray);
1659 ArrayWrapper<Scalar,ArrayIn1, Rank<ArrayIn1 >::value, true>inArray1Wrap(inArray1);
1660 ArrayWrapper<Scalar,ArrayIn2, Rank<ArrayIn2 >::value, true>inArray2Wrap(inArray2);
1661 size_t inArray1Rank=getrank(inArray1);
1662#ifdef HAVE_INTREPID_DEBUG
1663 TEUCHOS_TEST_FOR_EXCEPTION( ( (getrank(inArray1) != getrank(inArray2)) || (getrank(inArray1) != getrank(diffArray)) ),
1664 std::invalid_argument,
1665 ">>> ERROR (RealSpaceTools::subtract): Array arguments must have identical ranks!");
1666 for (size_t i=0; i<getrank(inArray1); i++) {
1667 TEUCHOS_TEST_FOR_EXCEPTION( ( (static_cast<size_t>(inArray1.dimension(i)) != static_cast<size_t>(inArray2.dimension(i))) || (static_cast<size_t>(inArray1.dimension(i)) != static_cast<size_t>(diffArray.dimension(i))) ),
1668 std::invalid_argument,
1669 ">>> ERROR (RealSpaceTools::subtract): Dimensions of array arguments do not agree!");
1670 }
1671#endif
1672 if(inArray1Rank==5){
1673 for (size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++)
1674 for (size_t j=0; j<static_cast<size_t>(inArray1.dimension(1)); j++)
1675 for (size_t k=0; k<static_cast<size_t>(inArray1.dimension(2)); k++)
1676 for (size_t l=0; l<static_cast<size_t>(inArray1.dimension(3)); l++)
1677 for (size_t m=0; m<static_cast<size_t>(inArray1.dimension(4)); m++){
1678 diffArrayWrap(i,j,k,l,m) = inArray1Wrap(i,j,k,l,m)-inArray2Wrap(i,j,k,l,m);
1679 }
1680 }else if(inArray1Rank==4){
1681 for (size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++)
1682 for (size_t j=0; j<static_cast<size_t>(inArray1.dimension(1)); j++)
1683 for (size_t k=0; k<static_cast<size_t>(inArray1.dimension(2)); k++)
1684 for (size_t l=0; l<static_cast<size_t>(inArray1.dimension(3)); l++){
1685 diffArrayWrap(i,j,k,l) = inArray1Wrap(i,j,k,l)-inArray2Wrap(i,j,k,l);
1686 }
1687 }else if(inArray1Rank==3){
1688 for (size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++)
1689 for (size_t j=0; j<static_cast<size_t>(inArray1.dimension(1)); j++)
1690 for (size_t k=0; k<static_cast<size_t>(inArray1.dimension(2)); k++){
1691 diffArrayWrap(i,j,k) = inArray1Wrap(i,j,k)-inArray2Wrap(i,j,k);
1692 }
1693 }else if(inArray1Rank==2){
1694 for (size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++)
1695 for (size_t j=0; j<static_cast<size_t>(inArray1.dimension(1)); j++){
1696 diffArrayWrap(i,j) = inArray1Wrap(i,j)-inArray2Wrap(i,j);
1697 }
1698 }else if(inArray1Rank==1){
1699 for (size_t i=0; i<static_cast<size_t>(inArray1.dimension(0)); i++){
1700 diffArrayWrap(i) = inArray1Wrap(i)-inArray2Wrap(i);
1701
1702 }
1703 }
1704
1705}
1706
1707
1708template<class Scalar>
1709template<class ArrayDiff, class ArrayIn>
1710void RealSpaceTools<Scalar>::subtract(ArrayDiff & inoutDiffArray, const ArrayIn & inArray) {
1711 ArrayWrapper<Scalar,ArrayDiff, Rank<ArrayDiff >::value, false>inoutDiffArrayWrap(inoutDiffArray);
1712 ArrayWrapper<Scalar,ArrayIn, Rank<ArrayIn >::value, true>inArrayWrap(inArray);
1713 int inArrayRank=getrank(inArray);
1714#ifdef HAVE_INTREPID_DEBUG
1715 TEUCHOS_TEST_FOR_EXCEPTION( ( getrank(inArray) != getrank(inoutDiffArray) ),
1716 std::invalid_argument,
1717 ">>> ERROR (RealSpaceTools::subtract): Array arguments must have identical ranks!");
1718 for (size_t i=0; i<getrank(inArray); i++) {
1719 TEUCHOS_TEST_FOR_EXCEPTION( ( static_cast<size_t>(inArray.dimension(i)) != static_cast<size_t>(inoutDiffArray.dimension(i)) ),
1720 std::invalid_argument,
1721 ">>> ERROR (RealSpaceTools::subtract): Dimensions of array arguments do not agree!");
1722 }
1723#endif
1724
1725 if(inArrayRank==5){
1726 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++)
1727 for (size_t j=0; j<static_cast<size_t>(static_cast<size_t>(inArray.dimension(1))); j++)
1728 for (size_t k=0; k<static_cast<size_t>(static_cast<size_t>(inArray.dimension(2))); k++)
1729 for (size_t l=0; l<static_cast<size_t>(static_cast<size_t>(inArray.dimension(3))); l++)
1730 for (size_t m=0; m<static_cast<size_t>(static_cast<size_t>(inArray.dimension(4))); m++){
1731 inoutDiffArrayWrap(i,j,k,l,m) -= inArrayWrap(i,j,k,l,m);
1732 }
1733 }else if(inArrayRank==4){
1734 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++)
1735 for (size_t j=0; j<static_cast<size_t>(static_cast<size_t>(inArray.dimension(1))); j++)
1736 for (size_t k=0; k<static_cast<size_t>(static_cast<size_t>(inArray.dimension(2))); k++)
1737 for (size_t l=0; l<static_cast<size_t>(static_cast<size_t>(inArray.dimension(3))); l++){
1738 inoutDiffArrayWrap(i,j,k,l) -= inArrayWrap(i,j,k,l);
1739 }
1740 }else if(inArrayRank==3){
1741 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++)
1742 for (size_t j=0; j<static_cast<size_t>(static_cast<size_t>(inArray.dimension(1))); j++)
1743 for (size_t k=0; k<static_cast<size_t>(static_cast<size_t>(inArray.dimension(2))); k++){
1744 inoutDiffArrayWrap(i,j,k) -= inArrayWrap(i,j,k);
1745 }
1746 }else if(inArrayRank==2){
1747 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++)
1748 for (size_t j=0; j<static_cast<size_t>(static_cast<size_t>(inArray.dimension(1))); j++){
1749 inoutDiffArrayWrap(i,j) -= inArrayWrap(i,j);
1750 }
1751 }else if(inArrayRank==1){
1752 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++){
1753 inoutDiffArrayWrap(i) -= inArrayWrap(i);
1754
1755 }
1756 }
1757
1758}
1759
1760
1761template<class Scalar>
1762void RealSpaceTools<Scalar>::scale(Scalar* scaledArray, const Scalar* inArray, const int size, const Scalar scalar) {
1763 for (size_t i=0; i<size; i++) {
1764 scaledArray[i] = scalar*inArray[i];
1765 }
1766}
1767
1768
1769
1770template<class Scalar>
1771void RealSpaceTools<Scalar>::scale(Scalar* inoutScaledArray, const int size, const Scalar scalar) {
1772 for (size_t i=0; i<size; i++) {
1773 inoutScaledArray[i] *= scalar;
1774 }
1775}
1776
1777
1778
1779template<class Scalar>
1780template<class ArrayScaled, class ArrayIn>
1781void RealSpaceTools<Scalar>::scale(ArrayScaled & scaledArray, const ArrayIn & inArray, const Scalar scalar) {
1782#ifdef HAVE_INTREPID_DEBUG
1783 TEUCHOS_TEST_FOR_EXCEPTION( ( getrank(inArray) != getrank(scaledArray) ),
1784 std::invalid_argument,
1785 ">>> ERROR (RealSpaceTools::scale): Array arguments must have identical ranks!");
1786 for (size_t i=0; i<getrank(inArray); i++) {
1787 TEUCHOS_TEST_FOR_EXCEPTION( ( static_cast<size_t>(inArray.dimension(i)) != static_cast<size_t>(scaledArray.dimension(i)) ),
1788 std::invalid_argument,
1789 ">>> ERROR (RealSpaceTools::scale): Dimensions of array arguments do not agree!");
1790 }
1791#endif
1792
1793
1794 int inArrayRank=getrank(inArray);
1795 ArrayWrapper<Scalar,ArrayScaled, Rank<ArrayScaled >::value, false>scaledArrayWrap(scaledArray);
1796 ArrayWrapper<Scalar,ArrayIn, Rank<ArrayIn >::value, true>inArrayWrap(inArray);
1797 if(inArrayRank==5){
1798 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++)
1799 for (size_t j=0; j<static_cast<size_t>(static_cast<size_t>(inArray.dimension(1))); j++)
1800 for (size_t k=0; k<static_cast<size_t>(static_cast<size_t>(inArray.dimension(2))); k++)
1801 for (size_t l=0; l<static_cast<size_t>(static_cast<size_t>(inArray.dimension(3))); l++)
1802 for (size_t m=0; m<static_cast<size_t>(static_cast<size_t>(inArray.dimension(4))); m++){
1803 scaledArrayWrap(i,j,k,l,m) = scalar*inArrayWrap(i,j,k,l,m);
1804 }
1805 }else if(inArrayRank==4){
1806 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++)
1807 for (size_t j=0; j<static_cast<size_t>(static_cast<size_t>(inArray.dimension(1))); j++)
1808 for (size_t k=0; k<static_cast<size_t>(static_cast<size_t>(inArray.dimension(2))); k++)
1809 for (size_t l=0; l<static_cast<size_t>(static_cast<size_t>(inArray.dimension(3))); l++){
1810 scaledArrayWrap(i,j,k,l) = scalar*inArrayWrap(i,j,k,l);
1811 }
1812 }else if(inArrayRank==3){
1813 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++)
1814 for (size_t j=0; j<static_cast<size_t>(static_cast<size_t>(inArray.dimension(1))); j++)
1815 for (size_t k=0; k<static_cast<size_t>(static_cast<size_t>(inArray.dimension(2))); k++){
1816 scaledArrayWrap(i,j,k) = scalar*inArrayWrap(i,j,k);
1817 }
1818 }else if(inArrayRank==2){
1819 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++)
1820 for (size_t j=0; j<static_cast<size_t>(static_cast<size_t>(inArray.dimension(1))); j++){
1821 scaledArrayWrap(i,j) = scalar*inArrayWrap(i,j);
1822 }
1823 }else if(inArrayRank==1){
1824 for (size_t i=0; i<static_cast<size_t>(static_cast<size_t>(inArray.dimension(0))); i++){
1825 scaledArrayWrap(i) = scalar*inArrayWrap(i);
1826
1827 }
1828 }
1829}
1830
1831
1832
1833template<class Scalar>
1834template<class ArrayScaled>
1835void RealSpaceTools<Scalar>::scale(ArrayScaled & inoutScaledArray, const Scalar scalar) {
1836 // Intrepid::FieldContainer has size type int
1837 const int theSize = (int) inoutScaledArray.size();
1838 for (int i=0; i<theSize; i++) {
1839 inoutScaledArray[i] *= scalar;
1840 }
1841}
1842
1843
1844
1845
1846template<class Scalar>
1847Scalar RealSpaceTools<Scalar>::dot(const Scalar* inArray1, const Scalar* inArray2, const int size) {
1848 Scalar dot(0);
1849 for (int i=0; i<size; i++) {
1850 dot += inArray1[i]*inArray2[i];
1851 }
1852 return dot;
1853}
1854
1855
1856
1857template<class Scalar>
1858template<class ArrayVec1, class ArrayVec2>
1859Scalar RealSpaceTools<Scalar>::dot(const ArrayVec1 & inVec1, const ArrayVec2 & inVec2) {
1860#ifdef HAVE_INTREPID_DEBUG
1861 TEUCHOS_TEST_FOR_EXCEPTION( ( (getrank(inVec1) != 1) || (getrank(inVec2) != 1) ),
1862 std::invalid_argument,
1863 ">>> ERROR (RealSpaceTools::dot): Vector arguments must have rank 1!");
1864 TEUCHOS_TEST_FOR_EXCEPTION( ( static_cast<size_t>(inVec1.dimension(0)) != static_cast<size_t>(inVec2.dimension(0)) ),
1865 std::invalid_argument,
1866 ">>> ERROR (RealSpaceTools::dot): Dimensions of vector arguments must agree!");
1867#endif
1868 ArrayWrapper<Scalar,ArrayVec1, Rank<ArrayVec1 >::value, true>inVec1Wrap(inVec1);
1869 ArrayWrapper<Scalar,ArrayVec2, Rank<ArrayVec2 >::value, true>inVec2Wrap(inVec2);
1870 Scalar dot(0);
1871 for (size_t i=0; i<static_cast<size_t>(inVec1.dimension(0)); i++) {
1872 dot += inVec1Wrap(i)*inVec2Wrap(i);
1873 }
1874 return dot;
1875
1876}
1877
1878
1879
1880template<class Scalar>
1881template<class ArrayDot, class ArrayVec1, class ArrayVec2>
1882void RealSpaceTools<Scalar>::dot(ArrayDot & dotArray, const ArrayVec1 & inVecs1, const ArrayVec2 & inVecs2) {
1883
1884 size_t arrayRank = getrank(inVecs1);
1885
1886#ifdef HAVE_INTREPID_DEBUG
1887 TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != getrank(dotArray)+1 ),
1888 std::invalid_argument,
1889 ">>> ERROR (RealSpaceTools::dot): Ranks of norm and vector array arguments are incompatible!");
1890 TEUCHOS_TEST_FOR_EXCEPTION( ( arrayRank != getrank(inVecs2) ),
1891 std::invalid_argument,
1892 ">>> ERROR (RealSpaceTools::dot): Ranks of input vector arguments must be identical!");
1893 TEUCHOS_TEST_FOR_EXCEPTION( ( (arrayRank < 2) || (arrayRank > 3) ),
1894 std::invalid_argument,
1895 ">>> ERROR (RealSpaceTools::dot): Rank of input vector arguments must be 2 or 3!");
1896 for (size_t i=0; i<arrayRank; i++) {
1897 TEUCHOS_TEST_FOR_EXCEPTION( ( static_cast<size_t>(inVecs1.dimension(i)) != static_cast<size_t>(inVecs2.dimension(i)) ),
1898 std::invalid_argument,
1899 ">>> ERROR (RealSpaceTools::dot): Dimensions of input vector arguments do not agree!");
1900 }
1901 for (size_t i=0; i<arrayRank-1; i++) {
1902 TEUCHOS_TEST_FOR_EXCEPTION( ( static_cast<size_t>(inVecs1.dimension(i)) != static_cast<size_t>(dotArray.dimension(i)) ),
1903 std::invalid_argument,
1904 ">>> ERROR (RealSpaceTools::dot): Dimensions of dot-product and vector arrays do not agree!");
1905 }
1906#endif
1907 ArrayWrapper<Scalar,ArrayDot, Rank<ArrayDot >::value, false>dotArrayWrap(dotArray);
1908 ArrayWrapper<Scalar,ArrayVec1, Rank<ArrayVec1 >::value, true>inVecs1Wrap(inVecs1);
1909 ArrayWrapper<Scalar,ArrayVec2, Rank<ArrayVec2 >::value, true>inVecs2Wrap(inVecs2);
1910 size_t dim_i0 = 1; // first index dimension (e.g. cell index)
1911 size_t dim_i1 = 1; // second index dimension (e.g. point index)
1912 size_t dim = static_cast<size_t>(inVecs1.dimension(arrayRank-1)); // spatial dimension
1913
1914 // determine i0 and i1 dimensions
1915 switch(arrayRank) {
1916 case 3:
1917 dim_i0 = static_cast<size_t>(inVecs1.dimension(0));
1918 dim_i1 = static_cast<size_t>(inVecs1.dimension(1));
1919 for (size_t i0=0; i0<dim_i0; i0++) {
1920 for (size_t i1=0; i1<dim_i1; i1++) {
1921 Scalar dot(0);
1922 for (size_t i=0; i<dim; i++) {
1923 dot += inVecs1Wrap(i0,i1,i)*inVecs2Wrap(i0,i1,i);
1924 }
1925 dotArrayWrap(i0,i1) = dot;
1926 }
1927 }
1928 break;
1929 case 2:
1930 dim_i1 = static_cast<size_t>(inVecs1.dimension(0));
1931 for (size_t i1=0; i1<dim_i1; i1++) {
1932 Scalar dot(0);
1933 for (size_t i=0; i<dim; i++) {
1934 dot += inVecs1Wrap(i1,i)*inVecs2Wrap(i1,i);
1935 }
1936 dotArrayWrap(i1) = dot;
1937 }
1938 break;
1939 case 1:
1940 Scalar dot(0);
1941 for (size_t i=0; i<dim; i++) {
1942 dot += inVecs1Wrap(i)*inVecs2Wrap(i);
1943 }
1944 dotArrayWrap(0) = dot;
1945
1946 break;
1947 }
1948
1949}
1950
1951
1952
1953template<class Scalar>
1954void RealSpaceTools<Scalar>::matvec(Scalar* matVec, const Scalar* inMat, const Scalar* inVec, const size_t dim) {
1955 for (size_t i=0; i<dim; i++) {
1956 Scalar sumdot(0);
1957 for (size_t j=0; j<dim; j++) {
1958 sumdot += inMat[i*dim+j]*inVec[j];
1959 }
1960 matVec[i] = sumdot;
1961 }
1962}
1963
1964
1965
1966template<class Scalar>
1967template<class ArrayMatVec, class ArrayMat, class ArrayVec>
1968void RealSpaceTools<Scalar>::matvec(ArrayMatVec & matVecs, const ArrayMat & inMats, const ArrayVec & inVecs) {
1969 size_t matArrayRank = getrank(inMats);
1970
1971#ifdef HAVE_INTREPID_DEBUG
1972 TEUCHOS_TEST_FOR_EXCEPTION( ( matArrayRank != getrank(inVecs)+1 ),
1973 std::invalid_argument,
1974 ">>> ERROR (RealSpaceTools::matvec): Vector and matrix array arguments do not have compatible ranks!");
1975 TEUCHOS_TEST_FOR_EXCEPTION( ( (matArrayRank < 3) || (matArrayRank > 4) ),
1976 std::invalid_argument,
1977 ">>> ERROR (RealSpaceTools::matvec): Rank of matrix array must be 3 or 4!");
1978 TEUCHOS_TEST_FOR_EXCEPTION( ( getrank(matVecs) != getrank(inVecs) ),
1979 std::invalid_argument,
1980 ">>> ERROR (RealSpaceTools::matvec): Vector arrays must be have the same rank!");
1981 for (size_t i=0; i<matArrayRank-1; i++) {
1982 TEUCHOS_TEST_FOR_EXCEPTION( ( static_cast<size_t>(inMats.dimension(i)) != static_cast<size_t>(inVecs.dimension(i)) ),
1983 std::invalid_argument,
1984 ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector and matrix array arguments do not agree!");
1985 }
1986 for (size_t i=0; i<getrank(inVecs); i++) {
1987 TEUCHOS_TEST_FOR_EXCEPTION( ( static_cast<size_t>(matVecs.dimension(i)) != static_cast<size_t>(inVecs.dimension(i)) ),
1988 std::invalid_argument,
1989 ">>> ERROR (RealSpaceTools::matvec): Dimensions of vector array arguments do not agree!");
1990 }
1991 TEUCHOS_TEST_FOR_EXCEPTION( ( static_cast<size_t>(inMats.dimension(matArrayRank-2)) != static_cast<size_t>(inMats.dimension(matArrayRank-1)) ),
1992 std::invalid_argument,
1993 ">>> ERROR (RealSpaceTools::matvec): Matrices are not square!");
1994#endif
1995 ArrayWrapper<Scalar,ArrayMatVec, Rank<ArrayMatVec >::value, false>matVecsWrap(matVecs);
1996 ArrayWrapper<Scalar,ArrayMat, Rank<ArrayMat >::value, true>inMatsWrap(inMats);
1997 ArrayWrapper<Scalar,ArrayVec, Rank<ArrayVec >::value, true>inVecsWrap(inVecs);
1998 size_t dim_i0 = 1; // first index dimension (e.g. cell index)
1999 size_t dim_i1 = 1; // second index dimension (e.g. point index)
2000 size_t dim = static_cast<size_t>(inMats.dimension(matArrayRank-2)); // spatial dimension
2001
2002 // determine i0 and i1 dimensions
2003 switch(matArrayRank) {
2004 case 4:
2005 dim_i0 = static_cast<size_t>(inMats.dimension(0));
2006 dim_i1 = static_cast<size_t>(inMats.dimension(1));
2007 for (size_t i0=0; i0<dim_i0; i0++) {
2008 for (size_t i1=0; i1<dim_i1; i1++) {
2009 for (size_t i=0; i<dim; i++) {
2010 Scalar sumdot(0);
2011 for (size_t j=0; j<dim; j++) {
2012 sumdot += inMatsWrap(i0,i1,i,j)*inVecsWrap(i0,i1,j);
2013 }
2014 matVecsWrap(i0,i1,i) = sumdot;
2015 }
2016 }
2017 }
2018 break;
2019 case 3:
2020 dim_i1 = static_cast<size_t>(inMats.dimension(0));
2021
2022 for (size_t i1=0; i1<dim_i1; i1++) {
2023 for (size_t i=0; i<dim; i++) {
2024 Scalar sumdot(0);
2025 for (size_t j=0; j<dim; j++) {
2026 sumdot += inMatsWrap(i1,i,j)*inVecsWrap(i1,j);
2027 }
2028 matVecsWrap(i1,i) = sumdot;
2029 }
2030 }
2031
2032 break;
2033
2034 }
2035
2036}
2037template<class Scalar>
2038template<class ArrayVecProd, class ArrayIn1, class ArrayIn2>
2039void RealSpaceTools<Scalar>::vecprod(ArrayVecProd & vecProd, const ArrayIn1 & inLeft, const ArrayIn2 & inRight) {
2040 ArrayWrapper<Scalar,ArrayVecProd, Rank<ArrayVecProd >::value, false>vecProdWrap(vecProd);
2041 ArrayWrapper<Scalar,ArrayIn1, Rank<ArrayIn1 >::value, true>inLeftWrap(inLeft);
2042 ArrayWrapper<Scalar,ArrayIn2, Rank<ArrayIn2 >::value, true>inRightWrap(inRight);
2043#ifdef HAVE_INTREPID_DEBUG
2044 /*
2045 * Check array rank and spatial dimension range (if applicable)
2046 * (1) all array arguments are required to have matching dimensions and rank: (D), (I0,D) or (I0,I1,D)
2047 * (2) spatial dimension should be 2 or 3
2048 */
2049 std::string errmsg = ">>> ERROR (RealSpaceTools::vecprod):";
2050
2051 // (1) check rank range on inLeft and then compare the other arrays with inLeft
2052 TEUCHOS_TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inLeft, 1,3), std::invalid_argument, errmsg);
2053 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inLeft, inRight), std::invalid_argument, errmsg);
2054 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inLeft, vecProd), std::invalid_argument, errmsg);
2055
2056 // (2) spatial dimension ordinal = array rank - 1. Suffices to check only one array because we just
2057 // checked whether or not the arrays have matching dimensions.
2058 TEUCHOS_TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inLeft, getrank(inLeft) - 1, 2,3), std::invalid_argument, errmsg);
2059
2060#endif
2061
2062 int spaceDim = static_cast<size_t>(inLeft.dimension(getrank(inLeft) - 1));
2063
2064 switch(getrank(inLeft) ){
2065
2066 case 1:
2067 {
2068 vecProdWrap(0) = inLeftWrap(1)*inRightWrap(2) - inLeftWrap(2)*inRightWrap(1);
2069 vecProdWrap(1) = inLeftWrap(2)*inRightWrap(0) - inLeftWrap(0)*inRightWrap(2);
2070 vecProdWrap(2) = inLeftWrap(0)*inRightWrap(1) - inLeftWrap(1)*inRightWrap(0);
2071 }
2072 break;
2073
2074 case 2:
2075 {
2076 size_t dim0 = static_cast<size_t>(inLeft.dimension(0));
2077 if(spaceDim == 3) {
2078 for(size_t i0 = 0; i0 < dim0; i0++){
2079 vecProdWrap(i0, 0) = inLeftWrap(i0, 1)*inRightWrap(i0, 2) - inLeftWrap(i0, 2)*inRightWrap(i0, 1);
2080 vecProdWrap(i0, 1) = inLeftWrap(i0, 2)*inRightWrap(i0, 0) - inLeftWrap(i0, 0)*inRightWrap(i0, 2);
2081 vecProdWrap(i0, 2) = inLeftWrap(i0, 0)*inRightWrap(i0, 1) - inLeftWrap(i0, 1)*inRightWrap(i0, 0);
2082 }// i0
2083 } //spaceDim == 3
2084 else if(spaceDim == 2){
2085 for(size_t i0 = 0; i0 < dim0; i0++){
2086 // vecprod is scalar - do we still want result to be (i0,i1,D)?
2087 vecProdWrap(i0, 0) = inLeftWrap(i0, 0)*inRightWrap(i0, 1) - inLeftWrap(i0, 1)*inRightWrap(i0, 0);
2088 }// i0
2089 }// spaceDim == 2
2090 }// case 2
2091 break;
2092
2093 case 3:
2094 {
2095 size_t dim0 = static_cast<size_t>(inLeft.dimension(0));
2096 size_t dim1 = static_cast<size_t>(inLeft.dimension(1));
2097 if(spaceDim == 3) {
2098 for(size_t i0 = 0; i0 < dim0; i0++){
2099 for(size_t i1 = 0; i1 < dim1; i1++){
2100 vecProdWrap(i0, i1, 0) = inLeftWrap(i0, i1, 1)*inRightWrap(i0, i1, 2) - inLeftWrap(i0, i1, 2)*inRightWrap(i0, i1, 1);
2101 vecProdWrap(i0, i1, 1) = inLeftWrap(i0, i1, 2)*inRightWrap(i0, i1, 0) - inLeftWrap(i0, i1, 0)*inRightWrap(i0, i1, 2);
2102 vecProdWrap(i0, i1, 2) = inLeftWrap(i0, i1, 0)*inRightWrap(i0, i1, 1) - inLeftWrap(i0, i1, 1)*inRightWrap(i0, i1, 0);
2103 }// i1
2104 }// i0
2105 } //spaceDim == 3
2106 else if(spaceDim == 2){
2107 for(size_t i0 = 0; i0 < dim0; i0++){
2108 for(size_t i1 = 0; i1 < dim1; i1++){
2109 // vecprod is scalar - do we still want result to be (i0,i1,D)?
2110 vecProdWrap(i0, i1, 0) = inLeftWrap(i0, i1, 0)*inRightWrap(i0, i1, 1) - inLeftWrap(i0, i1, 1)*inRightWrap(i0, i1, 0);
2111 }// i1
2112 }// i0
2113 }// spaceDim == 2
2114 } // case 3
2115 break;
2116
2117 default:
2118 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
2119 ">>> ERROR (RealSpaceTools::vecprod): rank-1,2,3 arrays required");
2120 }
2121
2122}
2123
2124
2125} // namespace Intrepid
bool requireRankRange(std::string &errmsg, const Array &array, const int lowerBound, const int upperBound)
Checks if the rank of the array argument is in the required range.
bool requireDimensionRange(std::string &errmsg, const Array &array, const int dim, const int lowerBound, const int upperBound)
Checks if the specified array dimension is in the required range.
bool requireDimensionMatch(std::string &errmsg, const Array1 &array1, const int a1_dim0, const Array2 &array2, const int a2_dim0)
Checks arrays for a single matching dimension.
Implementation of basic linear algebra functionality in Euclidean space.
static void subtract(Scalar *diffArray, const Scalar *inArray1, const Scalar *inArray2, const int size)
Subtracts contiguous data inArray2 from inArray1 of size size: diffArray = inArray1 - inArray2.
static void inverse(Scalar *inverseMat, const Scalar *inMat, const size_t dim)
Computes inverse of the square matrix inMat of size dim by dim.
static Scalar dot(const Scalar *inArray1, const Scalar *inArray2, const int size)
Computes dot product of contiguous data inArray1 and inArray2 of size size.
static void matvec(Scalar *matVec, const Scalar *inMat, const Scalar *inVec, const size_t dim)
Matrix-vector left multiply using contiguous data: matVec = inMat * inVec.
static void vecprod(ArrayVecProd &vecProd, const ArrayIn1 &inLeft, const ArrayIn2 &inRight)
Vector product using multidimensional arrays: vecProd = inVecLeft x inVecRight
static void scale(Scalar *scaledArray, const Scalar *inArray, const int size, const Scalar scalar)
Multiplies contiguous data inArray of size size by a scalar (componentwise): scaledArray = scalar * ...
static Scalar vectorNorm(const Scalar *inVec, const size_t dim, const ENorm normType)
Computes norm (1, 2, infinity) of the vector inVec of size dim.
static void add(Scalar *sumArray, const Scalar *inArray1, const Scalar *inArray2, const int size)
Adds contiguous data inArray1 and inArray2 of size size: sumArray = inArray1 + inArray2.
static void transpose(Scalar *transposeMat, const Scalar *inMat, const size_t dim)
Computes transpose of the square matrix inMat of size dim by dim.
static void absval(Scalar *absArray, const Scalar *inArray, const int size)
Computes absolute value of contiguous input data inArray of size size.