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 
50 namespace Intrepid {
51 
52 
53 
54 template<class Scalar>
55 void 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 
63 template<class Scalar>
64 void 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 
72 template<class Scalar>
73 template<class ArrayAbs, class ArrayIn>
74 void 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 
128 template<class Scalar>
129 template<class ArrayInOut>
130 void 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 
138 template<class Scalar>
139 Scalar 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 
168 template<class Scalar>
169 template<class ArrayIn>
170 Scalar 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 /*
297 template<class Scalar>
298 template<class ArrayIn>
299 Scalar 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 }*/
336 template<class Scalar>
337 template<class ArrayNorm, class ArrayIn>
338 void 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 
473 template<class Scalar>
474 template<class ArrayNorm, class ArrayIn>
475 void 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 
575 template<class Scalar>
576 void 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 
588 template<class Scalar>
589 template<class ArrayTranspose, class ArrayIn>
590 void 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 
652 template<class Scalar>
653 void 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 
755 template<class Scalar>
756 template<class ArrayInverse, class ArrayIn>
757 void 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 
1171 template<class Scalar>
1172 Scalar 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 
1238 template<class Scalar>
1239 template<class ArrayIn>
1240 Scalar 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 
1319 template<class Scalar>
1320 template<class ArrayDet, class ArrayIn>
1321 void 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 }
1497 break;
1498 
1499  } // switch (dim)
1500 
1501  }
1502 
1503 
1504 
1505 template<class Scalar>
1506 void 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 
1514 template<class Scalar>
1515 void 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 
1523 template<class Scalar>
1524 template<class ArraySum, class ArrayIn1, class ArrayIn2>
1525 void 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 
1582 template<class Scalar>
1583 template<class ArraySum, class ArrayIn>
1584 void 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 
1637 template<class Scalar>
1638 void 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 
1646 template<class Scalar>
1647 void 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 
1655 template<class Scalar>
1656 template<class ArrayDiff, class ArrayIn1, class ArrayIn2>
1657 void 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 
1708 template<class Scalar>
1709 template<class ArrayDiff, class ArrayIn>
1710 void 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 
1761 template<class Scalar>
1762 void 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 
1770 template<class Scalar>
1771 void 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 
1779 template<class Scalar>
1780 template<class ArrayScaled, class ArrayIn>
1781 void 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 
1833 template<class Scalar>
1834 template<class ArrayScaled>
1835 void 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 
1846 template<class Scalar>
1847 Scalar 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 
1857 template<class Scalar>
1858 template<class ArrayVec1, class ArrayVec2>
1859 Scalar 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 
1880 template<class Scalar>
1881 template<class ArrayDot, class ArrayVec1, class ArrayVec2>
1882 void 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 
1953 template<class Scalar>
1954 void 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 
1966 template<class Scalar>
1967 template<class ArrayMatVec, class ArrayMat, class ArrayVec>
1968 void 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 }
2037 template<class Scalar>
2038 template<class ArrayVecProd, class ArrayIn1, class ArrayIn2>
2039 void 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
Intrepid::RealSpaceTools
Implementation of basic linear algebra functionality in Euclidean space.
Definition: Intrepid_RealSpaceTools.hpp:65
Intrepid::RealSpaceTools::transpose
static void transpose(Scalar *transposeMat, const Scalar *inMat, const size_t dim)
Computes transpose of the square matrix inMat of size dim by dim.
Definition: Intrepid_RealSpaceToolsDef.hpp:576
Intrepid::RealSpaceTools::add
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.
Definition: Intrepid_RealSpaceToolsDef.hpp:1506
Intrepid::RealSpaceTools::dot
static Scalar dot(const Scalar *inArray1, const Scalar *inArray2, const int size)
Computes dot product of contiguous data inArray1 and inArray2 of size size.
Definition: Intrepid_RealSpaceToolsDef.hpp:1847
Intrepid::RealSpaceTools::vectorNorm
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.
Definition: Intrepid_RealSpaceToolsDef.hpp:139
Intrepid::RealSpaceTools::scale
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 * ...
Definition: Intrepid_RealSpaceToolsDef.hpp:1762
Intrepid::RealSpaceTools::inverse
static void inverse(Scalar *inverseMat, const Scalar *inMat, const size_t dim)
Computes inverse of the square matrix inMat of size dim by dim.
Definition: Intrepid_RealSpaceToolsDef.hpp:653
Intrepid::RealSpaceTools::absval
static void absval(Scalar *absArray, const Scalar *inArray, const int size)
Computes absolute value of contiguous input data inArray of size size.
Definition: Intrepid_RealSpaceToolsDef.hpp:55
Intrepid::RealSpaceTools::matvec
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.
Definition: Intrepid_RealSpaceToolsDef.hpp:1954
Intrepid::RealSpaceTools::vecprod
static void vecprod(ArrayVecProd &vecProd, const ArrayIn1 &inLeft, const ArrayIn2 &inRight)
Vector product using multidimensional arrays: vecProd = inVecLeft x inVecRight
Definition: Intrepid_RealSpaceToolsDef.hpp:2039
Intrepid::RealSpaceTools::subtract
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...
Definition: Intrepid_RealSpaceToolsDef.hpp:1638
Generated by   doxygen 1.8.5