Intrepid
test_04_kokkos.cpp
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 
48 #include "Intrepid_Utils.hpp"
49 #include "Intrepid_ArrayTools.hpp"
52 #include "Teuchos_oblackholestream.hpp"
53 #include "Teuchos_RCP.hpp"
54 #include "Teuchos_ScalarTraits.hpp"
55 #include <Kokkos_Core.hpp>
56 
57 using namespace std;
58 using namespace Intrepid;
59 
60 #define INTREPID_TEST_COMMAND( S ) \
61 { \
62  try { \
63  S ; \
64  } \
65  catch (std::logic_error err) { \
66  *outStream << "Expected Error ----------------------------------------------------------------\n"; \
67  *outStream << err.what() << '\n'; \
68  *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
69  }; \
70 }
71 
72 
73 int main(int argc, char *argv[]) {
74  Kokkos::initialize();
75  // This little trick lets us print to std::cout only if
76  // a (dummy) command-line argument is provided.
77  int iprint = argc - 1;
78  Teuchos::RCP<std::ostream> outStream;
79  Teuchos::oblackholestream bhs; // outputs nothing
80  if (iprint > 0)
81  outStream = Teuchos::rcp(&std::cout, false);
82  else
83  outStream = Teuchos::rcp(&bhs, false);
84 
85  // Save the format state of the original std::cout.
86  Teuchos::oblackholestream oldFormatState;
87  oldFormatState.copyfmt(std::cout);
88 
89  *outStream \
90  << "===============================================================================\n" \
91  << "| |\n" \
92  << "| Unit Test (ArrayTools) |\n" \
93  << "| |\n" \
94  << "| 1) Array operations: multiplication, contractions |\n" \
95  << "| |\n" \
96  << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
97  << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
98  << "| |\n" \
99  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
100  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
101  << "| |\n" \
102  << "===============================================================================\n";
103 
104 
105  int errorFlag = 0;
106 #ifdef HAVE_INTREPID_DEBUG
107  int beginThrowNumber = Teuchos::TestForException_getThrowNumber();
108  int endThrowNumber = beginThrowNumber + 39 + 18 + 28 + 26 + 34 + 37 + 46 + 45;
109 #endif
110  typedef ArrayTools art;
111  typedef RealSpaceTools<double> rst;
112  #ifdef HAVE_INTREPID_DEBUG
113  ArrayTools atools;
114  /************************************************************************************************
115  * *
116  * Exception tests: should only run when compiled in DEBUG mode *
117  * *
118  ***********************************************************************************************/
119  int C = 12, C1 = 10;
120  int P = 8, P1 = 16;
121  int F = 4, F1 = 8;
122  int D1 = 1, D2 = 2, D3 = 3, D4 = 4;
123 
124  Kokkos::View<double*> fc_P ("fc_P",P);
125  Kokkos::View<double**> fc_P_D1 ("fc_P_D1",P, D1);
126  Kokkos::View<double**> fc_P_D2 ("fc_P_D2",P, D2);
127  Kokkos::View<double**> fc_P_D3 ("fc_P_D3",P, D3);
128  Kokkos::View<double***> fc_P_D2_D2 ("fc_P_D2_D2",P, D2, D2);
129  Kokkos::View<double**> fc_P1_D3 ("fc_P1_D3",P1, D3);
130  Kokkos::View<double***> fc_P1_D2_D2 ("fc_P1_D2_D2",P1, D2, D2);
131  Kokkos::View<double***> fc_P1_D3_D3 ("fc_P1_D3_D3",P1, D3, D3);
132 
133  Kokkos::View<double*> fc_C ("fc_C",C);
134  Kokkos::View<double**> fc_C1_P ("fc_C1_P",C1,P);
135  Kokkos::View<double**> fc_C_P1 ("fc_C_P1",C, P1);
136  Kokkos::View<double**> fc_C_P ("fc_C_P",C, P);
137  Kokkos::View<double***> fc_C_P_D1 ("fc_C_P_D1",C, P, D1);
138  Kokkos::View<double***> fc_C_P_D2 ("fc_C_P_D2",C, P, D2);
139  Kokkos::View<double***> fc_C_P_D3 ("fc_C_P_D3",C, P, D3);
140  Kokkos::View<double***> fc_C_P_D4 ("fc_C_P_D4",C, P, D4);
141  Kokkos::View<double***> fc_C1_P_D2 ("fc_C1_P_D2",C1,P, D2);
142  Kokkos::View<double***> fc_C1_P_D3 ("fc_C1_P_D3",C1,P, D3);
143  Kokkos::View<double***> fc_C_P1_D1 ("fc_C_P1_D1",C, P1,D1);
144  Kokkos::View<double***> fc_C_P1_D2 ("fc_C_P1_D2",C, P1,D2);
145  Kokkos::View<double***> fc_C_P1_D3 ("fc_C_P1_D3",C, P1,D3);
146  Kokkos::View<double****> fc_C_P_D1_D1 ("fc_C_P_D1_D1",C, P, D1, D1);
147  Kokkos::View<double****> fc_C_P_D1_D2 ("fc_C_P_D1_D2",C, P, D1, D2);
148  Kokkos::View<double****> fc_C_P_D1_D3 ("fc_C_P_D1_D3",C, P, D1, D3);
149  Kokkos::View<double****> fc_C_P_D2_D2 ("fc_C_P_D2_D2",C, P, D2, D2);
150  Kokkos::View<double****> fc_C_P_D3_D1 ("fc_C_P_D3_D1",C, P, D3, D1);
151  Kokkos::View<double****> fc_C_P_D3_D2 ("fc_C_P_D3_D2",C, P, D3, D2);
152  Kokkos::View<double****> fc_C_P_D2_D3 ("fc_C_P_D2_D3",C, P, D2, D3);
153  Kokkos::View<double****> fc_C_P_D3_D3 ("fc_C_P_D3_D3",C, P, D3, D3);
154  Kokkos::View<double****> fc_C1_P_D3_D3 ("fc_C1_P_D3_D3",C1,P, D3, D3);
155  Kokkos::View<double****> fc_C1_P_D2_D2 ("fc_C1_P_D2_D2",C1,P, D2, D2);
156  Kokkos::View<double****> fc_C_P1_D2_D2 ("fc_C_P1_D2_D2",C, P1,D2, D2);
157  Kokkos::View<double****> fc_C_P1_D3_D3 ("fc_C_P1_D3_D3",C, P1,D3, D3);
158  Kokkos::View<double*****> fc_C_P_D3_D3_D3 ("fc_C_P_D3_D3_D3",C, P, D3, D3, D3);
159 
160  Kokkos::View<double**> fc_F_P ("fc_F_P",F, P);
161  Kokkos::View<double***> fc_F_P_D1 ("fc_F_P_D1",F, P, D1);
162  Kokkos::View<double***> fc_F_P_D2 ("fc_F_P_D2",F, P, D2);
163  Kokkos::View<double***> fc_F_P1_D2 ("fc_F_P1_D2",F, P1, D2);
164  Kokkos::View<double***> fc_F_P_D3 ("fc_F_P_D3",F, P, D3);
165  Kokkos::View<double****> fc_F_P_D3_D3 ("fc_F_P_D3_D3",F, P, D3, D3);
166  Kokkos::View<double***> fc_F1_P_D2 ("fc_F1_P_D2",F1,P, D2);
167  Kokkos::View<double***> fc_F1_P_D3 ("fc_F1_P_D3",F1,P, D3);
168  Kokkos::View<double****> fc_F1_P_D3_D3 ("fc_F1_P_D3_D3",F1,P, D3, D3);
169  Kokkos::View<double***> fc_F_P1_D3 ("fc_F_P1_D3",F, P1,D3);
170  Kokkos::View<double****> fc_F_P1_D3_D3 ("fc_F_P1_D3_D3",F, P1,D3, D3);
171  Kokkos::View<double****> fc_F_P_D2_D2 ("fc_F_P_D2_D2",F, P, D2, D2);
172  Kokkos::View<double****> fc_F_P1_D2_D2 ("fc_F_P1_D2_D2",F, P1,D2, D2);
173  Kokkos::View<double***> fc_C_F_P ("fc_C_F_P",C, F, P);
174  Kokkos::View<double***> fc_C1_F_P ("fc_C1_F_P",C1, F, P);
175  Kokkos::View<double***> fc_C_F1_P ("fc_C_F1_P",C, F1,P);
176  Kokkos::View<double***> fc_C_F_P1 ("fc_C_F_P1",C, F, P1);
177  Kokkos::View<double****> fc_C_F_P_D1 ("fc_C_F_P_D1",C, F, P, D1);
178  Kokkos::View<double****> fc_C_F_P_D2 ("fc_C_F_P_D2",C, F, P, D2);
179  Kokkos::View<double****> fc_C_F_P_D3 ("fc_C_F_P_D3",C, F, P, D3);
180  Kokkos::View<double****> fc_C1_F_P_D2 ("fc_C1_F_P_D2",C1, F, P,D2);
181  Kokkos::View<double****> fc_C1_F_P_D3 ("fc_C1_F_P_D3",C1, F, P,D3);
182  Kokkos::View<double****> fc_C_F1_P_D2 ("fc_C_F1_P_D2",C, F1,P, D2);
183  Kokkos::View<double****> fc_C_F1_P_D3 ("fc_C_F1_P_D3",C, F1,P, D3);
184  Kokkos::View<double*****> fc_C_F1_P_D3_D3 ("fc_C_F1_P_D3_D3",C, F1,P, D3, D3);
185  Kokkos::View<double****> fc_C_F_P1_D2 ("fc_C_F_P1_D2",C, F, P1,D2);
186  Kokkos::View<double****> fc_C_F_P1_D3 ("fc_C_F_P1_D3",C, F, P1,D3);
187  Kokkos::View<double*****> fc_C_F_P_D1_D1 ("fc_C_F_P_D1_D1",C, F, P, D1, D1);
188  Kokkos::View<double*****> fc_C_F_P_D2_D2 ("fc_C_F_P_D2_D2",C, F, P, D2, D2);
189  Kokkos::View<double*****> fc_C_F_P_D1_D3 ("fc_C_F_P_D1_D3",C, F, P, D1, D3);
190  Kokkos::View<double*****> fc_C_F_P_D2_D3 ("fc_C_F_P_D2_D3",C, F, P, D2, D3);
191  Kokkos::View<double*****> fc_C_F_P_D3_D1 ("fc_C_F_P_D3_D1",C, F, P, D3, D1);
192  Kokkos::View<double*****> fc_C_F_P_D3_D2 ("fc_C_F_P_D3_D2",C, F, P, D3, D2);
193  Kokkos::View<double*****> fc_C_F_P_D3_D3 ("fc_C_F_P_D3_D3",C, F, P, D3, D3);
194  Kokkos::View<double*****> fc_C_F_P1_D2_D2 ("fc_C_F_P1_D2_D2",C, F, P1,D2, D2);
195  Kokkos::View<double*****> fc_C_F_P1_D3_D3 ("fc_C_F_P1_D3_D3",C, F, P1,D3, D3);
196  Kokkos::View<double*****> fc_C1_F_P_D2_D2 ("fc_C1_F_P_D2_D2",C1,F, P, D2, D2);
197  Kokkos::View<double*****> fc_C1_F_P1_D2_D2 ("fc_C1_F_P1_D2_D2",C1,F, P1,D2, D2);
198  Kokkos::View<double*****> fc_C1_F_P_D3_D3 ("fc_C1_F_P_D3_D3",C1,F, P, D3, D3);
199 
200  Teuchos::Array<int> dimensions(6);
201  dimensions[0] = C;
202  dimensions[1] = F;
203  dimensions[2] = P;
204  dimensions[3] = D3;
205  dimensions[4] = D3;
206  dimensions[5] = D3;
207  FieldContainer<double> fc_C_F_P_D3_D3_D3 (dimensions);
208 
209 
210  *outStream \
211  << "\n"
212  << "===============================================================================\n"\
213  << "| TEST 1: crossProductDataField exceptions |\n"\
214  << "===============================================================================\n";
215 
216  try{
217  // 39 exceptions
218  // Test rank and D dimension of inputData
219  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P, fc_C_F_P_D3) );
220  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D2_D2, fc_C_F_P_D3) );
221  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D1, fc_C_F_P_D3) );
222 
223  // Test rank and D dimension of inputFields
224  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_F_P) );
225  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_C_F_P_D3_D3) );
226  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_C_F_P_D1) );
227  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_F_P_D1) );
228 
229  // Test rank of outputFields
230  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D2, fc_C_F_P_D2) );
231  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P_D3, fc_C_F_P_D3) );
232 
233  // Dimension cross-check: (1) inputData vs. inputFields
234  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_C1_F_P_D3) );
235  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_C_F_P1_D3) );
236  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_C_F_P_D2) );
237  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P_D2, fc_C1_F_P_D2) );
238  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P_D2, fc_C_F_P1_D2) );
239  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P_D2, fc_C_F_P_D3) );
240  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_F_P1_D3) );
241  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_F_P_D2) );
242  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P_D2, fc_F_P1_D2) );
243  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P_D2, fc_F_P_D3) );
244 
245  // Dimension cross-check: (2) outputFields vs. inputData
246  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C1_F_P, fc_C_P_D2, fc_C_F_P_D2) );
247  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P1, fc_C_P_D2, fc_C_F_P_D2) );
248  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C1_F_P, fc_C_P_D2, fc_F_P_D2) );
249  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P1, fc_C_P_D2, fc_F_P_D2) );
250  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C1_F_P_D3, fc_C_P_D3, fc_C_F_P_D3) );
251  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P1_D3, fc_C_P_D3, fc_C_F_P_D3) );
252  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D2, fc_C_P_D3, fc_C_F_P_D3) );
253  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C1_F_P_D3, fc_C_P_D3, fc_F_P_D3) );
254  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P1_D3, fc_C_P_D3, fc_F_P_D3) );
255  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D2, fc_C_P_D3, fc_F_P_D3) );
256 
257  // Dimension cross-check: (3) outputFields vs. inputFields
258  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C1_P_D2, fc_C1_F_P_D2) );
259  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P1_D2, fc_C_F_P1_D2) );
260  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P_D2, fc_C_F1_P_D2) );
261  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P1_D2, fc_F_P1_D2) );
262  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P, fc_C_P_D2, fc_F1_P_D2) );
263  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C1_P_D3, fc_C1_F_P_D3) );
264  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P1_D3, fc_C_F_P1_D3) );
265  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_C_F1_P_D3) );
266  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P1_D3, fc_F_P1_D3) );
267  INTREPID_TEST_COMMAND(atools.crossProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_F1_P_D3) );
268 
269  *outStream \
270  << "\n"
271  << "===============================================================================\n"\
272  << "| TEST 2: crossProductDataData exceptions |\n"\
273  << "===============================================================================\n";
274 
275  // 18 exceptions
276  // inputDataL is (C, P, D) and 2 <= D <= 3 is required
277  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P, fc_C_P_D3) );
278  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D2_D2, fc_C_P_D3) );
279  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D1, fc_C_P_D3) );
280 
281  // inputDataRight is (C, P, D) or (P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required.
282  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_C) );
283  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_C_P_D2_D2) );
284  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_C_P_D1) );
285 
286  // outputData is (C,P,D) in 3D and (C,P) in 2D => rank = inputDataLeft.dimension(2)
287  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D2, fc_C_P_D2) );
288  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P, fc_C_P_D3, fc_C_P_D3) );
289 
290  // Dimension cross-check (1):
291  // inputDataLeft(C,P,D) vs. inputDataRight(C,P,D): C, P, D must match
292  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C1_P_D3, fc_C_P_D3) );
293  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P1_D3, fc_C_P_D3) );
294  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_C_P_D2) );
295  // inputDataLeft(C, P,D) vs. inputDataRight(P,D): P, D must match
296  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P1_D3, fc_P_D3) );
297  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_P_D2) );
298 
299  // Dimension cross-check (2):
300  // in 2D: outputData(C,P) vs. inputDataLeft(C,P,D): dimensions C, P must match
301  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C1_P, fc_C_P_D2, fc_P_D2) );
302  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P1, fc_C_P_D2, fc_P_D2) );
303  // in 3D: outputData(C,P,D) vs. inputDataLeft(C,P,D): all dimensions C, P, D must match
304  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C1_P_D3, fc_C_P_D3, fc_P_D3) );
305  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P1_D3, fc_C_P_D3, fc_P_D3) );
306  INTREPID_TEST_COMMAND(atools.crossProductDataData<double>(fc_C_P_D2, fc_C_P_D3, fc_P_D3) );
307 
308  *outStream \
309  << "\n"
310  << "===============================================================================\n"\
311  << "| TEST 3: outerProductDataField exceptions |\n"\
312  << "===============================================================================\n";
313  // 28 exceptions
314  // Test rank and D dimension: inputData(C, P, D) and 2 <= D <= 3 is required
315  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P, fc_C_F_P_D3) );
316  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D2_D2, fc_C_F_P_D3) );
317  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D1, fc_C_F_P_D3) );
318 
319  // Test rank and D dimension: inputFields(C,F,P,D)/(F,P,D) and 2 <= (D=dimension(rank - 1)) <= 3 is required.
320  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_F_P) );
321  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_C_F_P_D3_D3) );
322  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_C_F_P_D1) );
323 
324  // Test rank and D dimension: outputFields(C,F,P,D,D)
325  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P, fc_C_P_D3, fc_C_F_P_D3) );
326  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_C_F_P_D3) );
327  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3_D3,fc_C_P_D3, fc_C_F_P_D3) );
328  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D1_D3, fc_C_P_D3, fc_C_F_P_D3) );
329  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D1, fc_C_P_D3, fc_C_F_P_D3) );
330  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D1_D1, fc_C_P_D3, fc_C_F_P_D3) );
331 
332  // Cross-check (2): outputFields(C,F,P,D,D) vs. inputData(C,P,D): dimensions C, P, D must match
333  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C1_P_D3, fc_C_F_P_D3) );
334  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1_D3, fc_C_F_P_D3) );
335  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D2, fc_C_F_P_D2) );
336 
337  // Cross-check (1): inputData(C,P,D) vs. inputFields(C,F,P,D): dimensions C, P, D must match
338  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_C1_F_P_D3) );
339  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_C_F_P1_D3) );
340  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_C_F_P_D2) );
341  // Cross-check (1): inputData(C,P,D) vs. inputFields(F,P,D): dimensions P, D must match
342  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_F_P1_D3) );
343  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_F_P_D2) );
344 
345  // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(C,F,P,D): dimensions C, F, P, D must match
346  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C1_P_D3, fc_C1_F_P_D3) );
347  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_C_F1_P_D3) );
348  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1_D3, fc_C_F_P1_D3) );
349  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D2, fc_C_F_P_D2) );
350  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D2_D3, fc_C_P_D2, fc_C_F_P_D2) );
351  // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(F,P,D): dimensions F, P, D must match
352  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_F1_P_D3) );
353  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1_D3, fc_F_P1_D3) );
354  INTREPID_TEST_COMMAND(atools.outerProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D2, fc_F_P_D2) );
355  *outStream \
356  << "\n"
357  << "===============================================================================\n"\
358  << "| TEST 4: outerProductDataData exceptions |\n"\
359  << "===============================================================================\n";
360  // 26 exceptions
361  // (1) inputDataLeft is (C, P, D) and 2 <= D <= 3 is required
362  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P, fc_C_P_D3) );
363  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P_D3) );
364  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D1, fc_C_P_D3) );
365 
366  // (2) inputDataRight is (C, P, D) or (P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required.
367  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C) );
368  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C_P_D3_D3) );
369  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C_P_D1) );
370  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_P_D1) );
371 
372  // (3) outputData is (C,P,D,D)
373  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_C_P_D3) );
374  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3_D3,fc_C_P_D3, fc_C_P_D3) );
375  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D1, fc_C_P_D3, fc_C_P_D3) );
376  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D1_D2, fc_C_P_D3, fc_C_P_D3) );
377 
378  // Cross-check (2): outputData(C,P,D,D) vs. inputDataLeft(C,P,D): dimensions C, P, D must match
379  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C1_P_D3, fc_P_D3) );
380  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P1_D3, fc_P_D3) );
381  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D2, fc_P_D3) );
382  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D2_D3, fc_C_P_D2, fc_P_D3) );
383  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D2, fc_C_P_D2, fc_P_D3) );
384 
385  // Cross-check (1): inputDataLeft(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match
386  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C1_P_D3) );
387  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C_P1_D3) );
388  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C_P_D2) );
389  // Cross-check (1): inputDataLeft(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match
390  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_P1_D3) );
391  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_P_D2) );
392 
393  // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(C,P,D): dimensions C, P, D must match
394  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C1_P_D3, fc_C1_P_D3) );
395  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P1_D3, fc_C_P1_D3) );
396  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D2, fc_C_P_D2) );
397  // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(P,D): dimensions P, D must match
398  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P1_D3, fc_P1_D3) );
399  INTREPID_TEST_COMMAND(atools.outerProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D2, fc_P_D2) );
400 
401  *outStream \
402  << "\n"
403  << "===============================================================================\n"\
404  << "| TEST 5: matvecProductDataField exceptions |\n"\
405  << "===============================================================================\n";
406  // 34 exceptions
407  // (1) inputData is (C,P), (C,P,D) or (C, P, D, D) and 1 <= D <= 3 is required
408  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C, fc_C_F_P_D3) );
409  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D4, fc_C_F_P_D3) );
410  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3_D3, fc_C_F_P_D3) );
411  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D1, fc_C_F_P_D3) );
412  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D1_D3, fc_C_F_P_D3) );
413 
414  // (2) inputFields is (C, F, P, D) or (F, P, D) and 1 <= (D=dimension(rank - 1)) <= 3 is required.
415  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_P_D3) );
416  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C_F_P_D3_D3) );
417  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C_F_P_D1) );
418  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_F_P_D1) );
419  // (3) outputFields is (C,F,P,D)
420  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P_D3) );
421  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P, fc_C_P_D3_D3, fc_C_F_P_D3) );
422  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D1, fc_C_P_D3_D3, fc_C_F_P_D3) );
423 
424  // Cross-check (2): outputFields(C,F,P,D) vs. inputData(C,P,D) and (C,P,D,D): dimensions C, P, D must match
425  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C1_P_D3_D3, fc_C_F_P_D3) );
426  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P1_D3_D3, fc_C_F_P_D3) );
427  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D2_D2, fc_C_F_P_D3) );
428  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C1_P, fc_C_F_P_D3) );
429  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P1, fc_C_F_P_D3) );
430  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C1_P_D3, fc_C_F_P_D3) );
431  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P1_D3, fc_C_F_P_D3) );
432  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D2, fc_C_F_P_D3) );
433 
434  // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(C,F,P,D): dimensions C, P, D must match
435  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C1_F_P_D3) );
436  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C_F_P1_D3) );
437  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C_F_P_D2) );
438  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_C_F_P_D2) );
439 
440  // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(F,P,D): dimensions P, D must match
441  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_F_P1_D3) );
442  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_F_P_D2) );
443  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3, fc_F_P_D2) );
444 
445  // Cross-check (3): outputFields(C,F,P,D) vs. inputFields(C,F,P,D): all dimensions C, F, P, D must match
446  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C1_F_P_D3) );
447  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C_F1_P_D3) );
448  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P1_D3_D3, fc_C_F_P1_D3) );
449  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D2, fc_C_P_D2_D2, fc_C_F_P_D3) );
450 
451  // Cross-check (3): outputFields(C,F,P,D) vs. inputFields(F,P,D): dimensions F, P, D must match
452  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C_F1_P_D3) );
453  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P1_D3_D3, fc_C_F_P1_D3) );
454  INTREPID_TEST_COMMAND(atools.matvecProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C_F_P_D2) );
455 
456  *outStream \
457  << "\n"
458  << "===============================================================================\n"\
459  << "| TEST 6: matvecProductDataData exceptions |\n"\
460  << "===============================================================================\n";
461  // 37 exceptions
462  // (1) inputDataLeft is (C,P), (C,P,D) or (C,P,D,D) and 1 <= D <= 3 is required
463  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C, fc_C_P_D3) );
464  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3_D3, fc_C_P_D3) );
465  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D1, fc_C_P_D3) );
466  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D1_D3, fc_C_P_D3) );
467 
468  // (2) inputDataRight is (C, P, D) or (P, D) and 1 <= (D=dimension(rank - 1)) <= 3 is required.
469  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_C_P_D3_D3) );
470  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_P) );
471  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_P_D1) );
472  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_C_P_D1) );
473 
474  // (3) outputData is (C,P,D)
475  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P, fc_C_P_D3_D3, fc_C_P_D3) );
476  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P_D3) );
477  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D1, fc_C_P_D3_D3, fc_C_P_D3) );
478 
479  // Cross-check (2): outputData(C,P,D) vs. inputDataLeft(C,P), (C,P,D), (C,P,D,D): dimensions C, P, D must match
480  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C1_P_D3, fc_C_P_D3_D3, fc_C_P_D3) );
481  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P1_D3, fc_C_P_D3_D3, fc_C_P_D3) );
482  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D2, fc_C_P_D3_D3, fc_C_P_D3) );
483  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C1_P_D3, fc_C_P, fc_C_P_D3) );
484  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P1_D3, fc_C_P, fc_C_P_D3) );
485  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C1_P_D3, fc_C_P_D3, fc_C_P_D3) );
486  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P1_D3, fc_C_P_D3, fc_C_P_D3) );
487  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D2, fc_C_P_D3, fc_C_P_D3) );
488 
489  // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(C,P,D): dimensions C, P, D must match
490  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_C1_P_D3) );
491  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_C_P1_D3) );
492  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_C_P_D2) );
493  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P, fc_C1_P_D3) );
494  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P, fc_C_P1_D3) );
495  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_C1_P_D3) );
496  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_C_P1_D3) );
497  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_C_P_D2) );
498 
499  // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(P,D): dimensions P, D must match
500  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_P1_D3) );
501  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_P_D2) );
502  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P, fc_P1_D3) );
503  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_P1_D3) );
504  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3, fc_P_D2) );
505 
506  // Cross-check (3): outputData(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match
507  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C1_P_D3, fc_C1_P_D3_D3, fc_C_P_D3) );
508  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P1_D3, fc_C_P1_D3_D3, fc_C_P_D3) );
509  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D2, fc_C_P_D3_D3, fc_C_P_D3) );
510 
511  // Cross-check (3): outputData(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match
512  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P1_D3, fc_C_P1_D3_D3, fc_P_D3) );
513  INTREPID_TEST_COMMAND(atools.matvecProductDataData<double>(fc_C_P_D3, fc_C_P_D3_D3, fc_P_D2) );
514 
515  *outStream \
516  << "\n"
517  << "===============================================================================\n"\
518  << "| TEST 7: matmatProductDataField exceptions |\n"\
519  << "===============================================================================\n";
520  // 46 exceptions
521  // (1) inputData is (C,P), (C,P,D), or (C,P,D,D) and 1 <= D <= 3 is required
522  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C, fc_C_F_P_D3_D3) );
523  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3_D3, fc_C_F_P_D3_D3) );
524  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D1_D3, fc_C_F_P_D3_D3) );
525  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D1, fc_C_F_P_D3_D3) );
526 
527  // (2) inputFields is (C,F,P,D,D) or (F,P,D,D) and 1 <= (dimension(rank-1), (rank-2)) <= 3 is required.
528  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P) );
529  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P_D3_D3_D3) );
530  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P_D1_D3) );
531  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P_D3_D1) );
532 
533  // (3) outputFields is (C,F,P,D,D)
534  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3, fc_C_P_D3_D3, fc_C_F_P_D3_D3) );
535  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3_D3, fc_C_P_D3_D3, fc_C_F_P_D3_D3) );
536  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D1_D3, fc_C_P_D3_D3, fc_C_F_P_D3_D3) );
537  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D1, fc_C_P_D3_D3, fc_C_F_P_D3_D3) );
538 
539  // Cross-check (2): outputFields(C,F,P,D,D) vs. inputData(C,P), (C,P,D), or (C,P,D,D): dimensions C, P, D must match
540  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C1_P_D3_D3, fc_C_F_P_D3_D3) );
541  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1_D3_D3, fc_C_F_P_D3_D3) );
542  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D2_D2, fc_C_F_P_D3_D3) );
543  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1_D2_D2, fc_C_F_P_D3_D3) );
544  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D2_D2, fc_C_P_D3_D3, fc_C_F_P_D3_D3) );
545  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C1_P, fc_C_F_P_D3_D3) );
546  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1, fc_C_F_P_D3_D3) );
547  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C1_P_D3, fc_C_F_P_D3_D3) );
548  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P1_D3, fc_C_F_P_D3_D3) );
549  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D1, fc_C_F_P_D3_D3) );
550 
551  // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(C,F,P,D,D): dimensions C, P, D must match
552  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C1_F_P_D3_D3) );
553  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P1_D3_D3) );
554  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P_D2_D2) );
555  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P1_D2_D2) );
556  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C1_F_P_D2_D2) );
557  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C1_F_P1_D2_D2) );
558  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P, fc_C1_F_P_D3_D3) );
559  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P, fc_C_F_P1_D3_D3) );
560  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_C1_F_P_D3_D3) );
561  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_C_F_P1_D3_D3) );
562  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_C_F_P_D2_D2) );
563 
564  // Cross-check (1): inputData(C,P), (C,P,D), or (C,P,D,D) vs. inputFields(F,P,D,D): dimensions P, D must match
565  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_F_P1_D3_D3) );
566  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_F_P_D2_D2) );
567  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_F_P1_D2_D2) );
568  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P, fc_F_P1_D3_D3) );
569  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_F_P1_D3_D3) );
570  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3, fc_F_P_D2_D2) );
571 
572  // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(C,F,P,D,D): all dimensions C, F, P, D must match
573  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C1_F_P_D3_D3) );
574  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F1_P_D3_D3) );
575  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P1_D3_D3) );
576  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_C_F_P_D2_D2) );
577 
578  // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(F,P,D,D): dimensions F, P, D must match
579  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_F1_P_D3_D3) );
580  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_F_P1_D3_D3) );
581  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_F_P_D3_D3, fc_C_P_D3_D3, fc_F_P_D2_D2) );
582  *outStream \
583  << "\n"
584  << "===============================================================================\n"\
585  << "| TEST 8: matmatProductDataData exceptions |\n"\
586  << "===============================================================================\n";
587  // 45 exceptions
588  // (1) inputDataLeft is (C,P), (C,P,D), or (C,P,D,D) and 2 <= D <= 3 is required
589  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C, fc_C_P_D3_D3) );
590  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3_D3, fc_C_P_D3_D3) );
591  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D1_D3, fc_C_P_D3_D3) );
592  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D1, fc_C_P_D3_D3) );
593 
594  // (2) inputDataRight is (C,P,D,D) or (P,D,D) and 1 <= (dimension(rank-1), (rank-2)) <= 3 is required.
595  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P, fc_C_P) );
596  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C_P_D3_D3_D3) );
597  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P_D1_D3) );
598  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P_D3_D1) );
599 
600  // (3) outputData is (C,P,D,D)
601  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3, fc_C_P, fc_C_P_D3_D3) );
602  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3_D3, fc_C_P_D3, fc_C_P_D3_D3) );
603  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D1_D3, fc_C_P_D3_D3, fc_C_P_D3_D3) );
604  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D1, fc_C_P_D3_D3, fc_C_P_D3_D3) );
605 
606  // Cross-check (2): outputData(C,P,D,D) vs. inputDataLeft(C,P), (C,P,D), or (C,P,D,D): dimensions C, P, D must match
607  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C1_P_D3_D3, fc_C_P_D3_D3) );
608  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P1_D3_D3, fc_C_P_D3_D3) );
609  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D2_D2, fc_C_P_D3_D3) );
610  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P1_D2_D2, fc_C_P_D3_D3) );
611  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D2_D2, fc_C_P_D3_D3, fc_C_P_D3_D3) );
612  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C1_P, fc_C_P_D3_D3) );
613  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P1, fc_C_P_D3_D3) );
614  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C1_P_D3, fc_C_P_D3_D3) );
615  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P1_D3, fc_C_P_D3_D3) );
616  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D1, fc_C_P_D3_D3) );
617 
618  // Cross-check (1): inputDataLeft(C,P), (C,P,D) or (C,P,D,D) vs. inputDataRight(C,P,D,D): dimensions C, P, D must match
619  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C1_P_D3_D3) );
620  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P1_D3_D3) );
621  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P_D2_D2) );
622  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P1_D2_D2) );
623  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C1_P_D2_D2) );
624  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P1_D2_D2) );
625  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P, fc_C1_P_D3_D3) );
626  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P, fc_C_P1_D3_D3) );
627  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C1_P_D3_D3) );
628  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C_P1_D3_D3) );
629  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_C_P_D2_D2) );
630 
631  // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(P,D,D): dimensions P, D must match
632  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_P1_D3_D3) );
633  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_P_D2_D2) );
634  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_P1_D2_D2) );
635  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_P_D3_D3, fc_C_P, fc_P1_D3_D3) );
636  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_P1_D3_D3) );
637  INTREPID_TEST_COMMAND(atools.matmatProductDataField<double>(fc_C_P_D3_D3, fc_C_P_D3, fc_P_D2_D2) );
638 
639  // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(C,P,D,D): all dimensions C, P, D must match
640  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C1_P_D3_D3) );
641  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C1_P_D3_D3) );
642  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P1_D3_D3) );
643  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_C_P_D2_D2) );
644 
645  // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(P,D,D): dimensions P, D must match
646  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_P_D2_D2) );
647  INTREPID_TEST_COMMAND(atools.matmatProductDataData<double>(fc_C_P_D3_D3, fc_C_P_D3_D3, fc_P1_D3_D3) );
648  }
649 
650  catch (std::logic_error err) {
651  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
652  *outStream << err.what() << '\n';
653  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
654  errorFlag = -1000;
655  };
656 
657  if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
658  errorFlag++;
659 #endif
660 
661 
662  /************************************************************************************************
663  * *
664  * Operation corectness tests *
665  * *
666  ***********************************************************************************************/
667 
668  try{
669  *outStream \
670  << "\n"
671  << "===============================================================================\n"\
672  << "| TEST 1.a: 3D crossProductDataField operations: (C,P,D) and (C,F,P,D) |\n"\
673  << "===============================================================================\n";
674  /*
675  * ijkData_1a ijkFields_1a Expected result in (C,F,P,D) array
676  * i,i i,i j,j, k,k 0, 0 k, k -j,-j
677  * j,j i,i j,j, k,k -k,-k 0, 0 i, i
678  * k,k i,i j,j, k,k j, j -i,-i 0, 0
679  */
680 
681  // input data is (C,P,D)
682  Kokkos::View<double***> ijkData_1a("ijkData_1a", 3, 2, 3);
683  // C=0 contains i
684  ijkData_1a(0, 0, 0) = 1.0; ijkData_1a(0, 0, 1) = 0.0; ijkData_1a(0, 0, 2) = 0.0;
685  ijkData_1a(0, 1, 0) = 1.0; ijkData_1a(0, 1, 1) = 0.0; ijkData_1a(0, 1, 2) = 0.0;
686  // C=1 contains j
687  ijkData_1a(1, 0, 0) = 0.0; ijkData_1a(1, 0, 1) = 1.0; ijkData_1a(1, 0, 2) = 0.0;
688  ijkData_1a(1, 1, 0) = 0.0; ijkData_1a(1, 1, 1) = 1.0; ijkData_1a(1, 1, 2) = 0.0;
689  // C=2 contains k
690  ijkData_1a(2, 0, 0) = 0.0; ijkData_1a(2, 0, 1) = 0.0; ijkData_1a(2, 0, 2) = 1.0;
691  ijkData_1a(2, 1, 0) = 0.0; ijkData_1a(2, 1, 1) = 0.0; ijkData_1a(2, 1, 2) = 1.0;
692 
693 
694  Kokkos::View<double****> ijkFields_1a("ijkFields_1a", 3, 3, 2, 3);
695  // C=0, F=0 is i
696  ijkFields_1a(0, 0, 0, 0) = 1.0; ijkFields_1a(0, 0, 0, 1) = 0.0; ijkFields_1a(0, 0, 0, 2) = 0.0;
697  ijkFields_1a(0, 0, 1, 0) = 1.0; ijkFields_1a(0, 0, 1, 1) = 0.0; ijkFields_1a(0, 0, 1, 2) = 0.0;
698  // C=0, F=1 is j
699  ijkFields_1a(0, 1, 0, 0) = 0.0; ijkFields_1a(0, 1, 0, 1) = 1.0; ijkFields_1a(0, 1, 0, 2) = 0.0;
700  ijkFields_1a(0, 1, 1, 0) = 0.0; ijkFields_1a(0, 1, 1, 1) = 1.0; ijkFields_1a(0, 1, 1, 2) = 0.0;
701  // C=0, F=2 is k
702  ijkFields_1a(0, 2, 0, 0) = 0.0; ijkFields_1a(0, 2, 0, 1) = 0.0; ijkFields_1a(0, 2, 0, 2) = 1.0;
703  ijkFields_1a(0, 2, 1, 0) = 0.0; ijkFields_1a(0, 2, 1, 1) = 0.0; ijkFields_1a(0, 2, 1, 2) = 1.0;
704 
705  // C=1, F=0 is i
706  ijkFields_1a(1, 0, 0, 0) = 1.0; ijkFields_1a(1, 0, 0, 1) = 0.0; ijkFields_1a(1, 0, 0, 2) = 0.0;
707  ijkFields_1a(1, 0, 1, 0) = 1.0; ijkFields_1a(1, 0, 1, 1) = 0.0; ijkFields_1a(1, 0, 1, 2) = 0.0;
708  // C=1, F=1 is j
709  ijkFields_1a(1, 1, 0, 0) = 0.0; ijkFields_1a(1, 1, 0, 1) = 1.0; ijkFields_1a(1, 1, 0, 2) = 0.0;
710  ijkFields_1a(1, 1, 1, 0) = 0.0; ijkFields_1a(1, 1, 1, 1) = 1.0; ijkFields_1a(1, 1, 1, 2) = 0.0;
711  // C=1, F=2 is k
712  ijkFields_1a(1, 2, 0, 0) = 0.0; ijkFields_1a(1, 2, 0, 1) = 0.0; ijkFields_1a(1, 2, 0, 2) = 1.0;
713  ijkFields_1a(1, 2, 1, 0) = 0.0; ijkFields_1a(1, 2, 1, 1) = 0.0; ijkFields_1a(1, 2, 1, 2) = 1.0;
714 
715  // C=2, F=0 is i
716  ijkFields_1a(2, 0, 0, 0) = 1.0; ijkFields_1a(2, 0, 0, 1) = 0.0; ijkFields_1a(2, 0, 0, 2) = 0.0;
717  ijkFields_1a(2, 0, 1, 0) = 1.0; ijkFields_1a(2, 0, 1, 1) = 0.0; ijkFields_1a(2, 0, 1, 2) = 0.0;
718  // C=2, F=1 is j
719  ijkFields_1a(2, 1, 0, 0) = 0.0; ijkFields_1a(2, 1, 0, 1) = 1.0; ijkFields_1a(2, 1, 0, 2) = 0.0;
720  ijkFields_1a(2, 1, 1, 0) = 0.0; ijkFields_1a(2, 1, 1, 1) = 1.0; ijkFields_1a(2, 1, 1, 2) = 0.0;
721  // C=2, F=2 is k
722  ijkFields_1a(2, 2, 0, 0) = 0.0; ijkFields_1a(2, 2, 0, 1) = 0.0; ijkFields_1a(2, 2, 0, 2) = 1.0;
723  ijkFields_1a(2, 2, 1, 0) = 0.0; ijkFields_1a(2, 2, 1, 1) = 0.0; ijkFields_1a(2, 2, 1, 2) = 1.0;
724 
725 
726  Kokkos::View<double****> outFields4("outFields4", 3, 3, 2, 3);
727  art::crossProductDataField<double>(outFields4, ijkData_1a, ijkFields_1a);
728 
729  // checks for C = 0
730  if( !(outFields4(0,0,0,0)==0.0 && outFields4(0,0,0,1)==0.0 && outFields4(0,0,0,2)==0.0 &&
731  outFields4(0,0,1,0)==0.0 && outFields4(0,0,1,1)==0.0 && outFields4(0,0,1,2)==0.0 ) ) {
732  *outStream << "\n\nINCORRECT crossProductDataField (1): i x i != 0; ";
733  errorFlag++;
734  }
735  if( !(outFields4(0,1,0,0)==0.0 && outFields4(0,1,0,1)==0.0 && outFields4(0,1,0,2)==1.0 &&
736  outFields4(0,1,1,0)==0.0 && outFields4(0,1,1,1)==0.0 && outFields4(0,1,1,2)==1.0 ) ) {
737  *outStream << "\n\nINCORRECT crossProductDataField (2): i x j != k; ";
738  errorFlag++;
739  }
740  if( !(outFields4(0,2,0,0)==0.0 && outFields4(0,2,0,1)==-1.0 && outFields4(0,2,0,2)==0.0 &&
741  outFields4(0,2,1,0)==0.0 && outFields4(0,2,1,1)==-1.0 && outFields4(0,2,1,2)==0.0 ) ) {
742  *outStream << "\n\nINCORRECT crossProductDataField (3): i x k != -j; ";
743  errorFlag++;
744  }
745 
746  // checks for C = 1
747  if( !(outFields4(1,0,0,0)==0.0 && outFields4(1,0,0,1)==0.0 && outFields4(1,0,0,2)==-1.0 &&
748  outFields4(1,0,1,0)==0.0 && outFields4(1,0,1,1)==0.0 && outFields4(1,0,1,2)==-1.0 ) ) {
749  *outStream << "\n\nINCORRECT crossProductDataField (4): j x i != -k; ";
750  errorFlag++;
751  }
752  if( !(outFields4(1,1,0,0)==0.0 && outFields4(1,1,0,1)==0.0 && outFields4(1,1,0,2)==0.0 &&
753  outFields4(1,1,1,0)==0.0 && outFields4(1,1,1,1)==0.0 && outFields4(1,1,1,2)==0.0 ) ) {
754  *outStream << "\n\nINCORRECT crossProductDataField (5): j x j != 0; ";
755  errorFlag++;
756  }
757  if( !(outFields4(1,2,0,0)==1.0 && outFields4(1,2,0,1)==0.0 && outFields4(1,2,0,2)==0.0 &&
758  outFields4(1,2,1,0)==1.0 && outFields4(1,2,1,1)==0.0 && outFields4(1,2,1,2)==0.0 ) ) {
759  *outStream << "\n\nINCORRECT crossProductDataField (6): j x k != i; ";
760  errorFlag++;
761  }
762 
763  // checks for C = 2
764  if( !(outFields4(2,0,0,0)==0.0 && outFields4(2,0,0,1)==1.0 && outFields4(2,0,0,2)==0.0 &&
765  outFields4(2,0,1,0)==0.0 && outFields4(2,0,1,1)==1.0 && outFields4(2,0,1,2)==0.0 ) ) {
766  *outStream << "\n\nINCORRECT crossProductDataField (7): k x i != j; ";
767  errorFlag++;
768  }
769  if( !(outFields4(2,1,0,0)==-1.0 && outFields4(2,1,0,1)==0.0 && outFields4(2,1,0,2)==0.0 &&
770  outFields4(2,1,1,0)==-1.0 && outFields4(2,1,1,1)==0.0 && outFields4(2,1,1,2)==0.0 ) ) {
771  *outStream << "\n\nINCORRECT crossProductDataField (8): k x j != -i; ";
772  errorFlag++;
773  }
774  if( !(outFields4(2,2,0,0)==0.0 && outFields4(2,2,0,1)==0.0 && outFields4(2,2,0,2)==0.0 &&
775  outFields4(2,2,1,0)==0.0 && outFields4(2,2,1,1)==0.0 && outFields4(2,2,1,2)==0.0 ) ) {
776  *outStream << "\n\nINCORRECT crossProductDataField (9): k x k != 0; ";
777  errorFlag++;
778  }
779 
780  *outStream \
781  << "\n"
782  << "===============================================================================\n"\
783  << "| TEST 1.b: 3D crossProductDataField operations: (C,P,D) and (F,P,D) |\n"\
784  << "===============================================================================\n";
785  /*
786  * ijkData_1b ijkFields_1b expected result in (C,F,P,D) array
787  * i,i,i i,j,k 0, k,-j
788  * j,j,j -k, 0, i
789  * k,k,k j,-i, 0
790  */
791 
792  // input data is (C,P,D)
793  Kokkos::View<double***> ijkData_1b("ijkData_1b", 3, 3, 3);
794  // C=0 contains i
795  ijkData_1b(0, 0, 0) = 1.0; ijkData_1b(0, 0, 1) = 0.0; ijkData_1b(0, 0, 2) = 0.0;
796  ijkData_1b(0, 1, 0) = 1.0; ijkData_1b(0, 1, 1) = 0.0; ijkData_1b(0, 1, 2) = 0.0;
797  ijkData_1b(0, 2, 0) = 1.0; ijkData_1b(0, 2, 1) = 0.0; ijkData_1b(0, 2, 2) = 0.0;
798  // C=1 contains j
799  ijkData_1b(1, 0, 0) = 0.0; ijkData_1b(1, 0, 1) = 1.0; ijkData_1b(1, 0, 2) = 0.0;
800  ijkData_1b(1, 1, 0) = 0.0; ijkData_1b(1, 1, 1) = 1.0; ijkData_1b(1, 1, 2) = 0.0;
801  ijkData_1b(1, 2, 0) = 0.0; ijkData_1b(1, 2, 1) = 1.0; ijkData_1b(1, 2, 2) = 0.0;
802  // C=2 contains k
803  ijkData_1b(2, 0, 0) = 0.0; ijkData_1b(2, 0, 1) = 0.0; ijkData_1b(2, 0, 2) = 1.0;
804  ijkData_1b(2, 1, 0) = 0.0; ijkData_1b(2, 1, 1) = 0.0; ijkData_1b(2, 1, 2) = 1.0;
805  ijkData_1b(2, 2, 0) = 0.0; ijkData_1b(2, 2, 1) = 0.0; ijkData_1b(2, 2, 2) = 1.0;
806 
807  // input fields are (F,P,D)
808  Kokkos::View<double***> ijkFields_1b("ijkFields_1b", 1, 3, 3);
809  // F=0 at 3 points is (i,j,k)
810  ijkFields_1b(0, 0, 0) = 1.0; ijkFields_1b(0, 0, 1) = 0.0; ijkFields_1b(0, 0, 2) = 0.0;
811  ijkFields_1b(0, 1, 0) = 0.0; ijkFields_1b(0, 1, 1) = 1.0; ijkFields_1b(0, 1, 2) = 0.0;
812  ijkFields_1b(0, 2, 0) = 0.0; ijkFields_1b(0, 2, 1) = 0.0; ijkFields_1b(0, 2, 2) = 1.0;
813 
814  // Output array is (C,F,P,D)
815  Kokkos::realloc(outFields4, 3, 1, 3, 3);
816  art::crossProductDataField<double>(outFields4, ijkData_1b, ijkFields_1b);
817 
818  // checks for C = 0
819  if( !(outFields4(0,0,0,0)==0.0 && outFields4(0,0,0,1)==0.0 && outFields4(0,0,0,2)==0.0) ) {
820  *outStream << "\n\nINCORRECT crossProductDataField (10): i x i != 0; ";
821  errorFlag++;
822  }
823  if( !(outFields4(0,0,1,0)==0.0 && outFields4(0,0,1,1)==0.0 && outFields4(0,0,1,2)==1.0) ) {
824  *outStream << "\n\nINCORRECT crossProductDataField (11): i x j != k; ";
825  errorFlag++;
826  }
827  if( !(outFields4(0,0,2,0)==0.0 && outFields4(0,0,2,1)==-1.0 && outFields4(0,0,2,2)==0.0) ) {
828  *outStream << "\n\nINCORRECT crossProductDataField (12): i x k != -j; ";
829  errorFlag++;
830  }
831 
832  // checks for C = 1
833  if( !(outFields4(1,0,0,0)==0.0 && outFields4(1,0,0,1)==0.0 && outFields4(1,0,0,2)==-1.0) ) {
834  *outStream << "\n\nINCORRECT crossProductDataField (13): j x i != -k; ";
835  errorFlag++;
836  }
837  if( !(outFields4(1,0,1,0)==0.0 && outFields4(1,0,1,1)==0.0 && outFields4(1,0,1,2)==0.0) ) {
838  *outStream << "\n\nINCORRECT crossProductDataField (14): j x j != 0; ";
839  errorFlag++;
840  }
841  if( !(outFields4(1,0,2,0)==1.0 && outFields4(1,0,2,1)==0.0 && outFields4(1,0,2,2)==0.0) ) {
842  *outStream << "\n\nINCORRECT crossProductDataField (15): j x k != i; ";
843  errorFlag++;
844  }
845 
846  // checks for C = 2
847  if( !(outFields4(2,0,0,0)==0.0 && outFields4(2,0,0,1)==1.0 && outFields4(2,0,0,2)==0.0) ) {
848  *outStream << "\n\nINCORRECT crossProductDataField (16): k x i != j; ";
849  errorFlag++;
850  }
851  if( !(outFields4(2,0,1,0)==-1.0 && outFields4(2,0,1,1)==0.0 && outFields4(2,0,1,2)==0.0) ) {
852  *outStream << "\n\nINCORRECT crossProductDataField (17): k x j != -i; ";
853  errorFlag++;
854  }
855  if( !(outFields4(2,0,2,0)==0.0 && outFields4(2,0,2,1)==0.0 && outFields4(2,0,2,2)==0.0) ) {
856  *outStream << "\n\nINCORRECT crossProductDataField (18): k x k != 0; ";
857  errorFlag++;
858  }
859 
860  *outStream \
861  << "\n"
862  << "===============================================================================\n"\
863  << "| TEST 1.c: 2D crossProductDataField operations: (C,P,D) and (C,F,P,D) |\n"\
864  << "===============================================================================\n";
865 
866  /*
867  * ijData_1c ijFields_1c expected result in (C,F,P) array
868  * i,i i,i j,j 0, 0 1, 1
869  * j,j i,i j,j -1,-1 0, 0
870  */
871  // input data is (C,P,D)
872  Kokkos::View<double***> ijData_1c("ijData_1c", 2, 2, 2);
873  // C=0 contains i
874  ijData_1c(0, 0, 0) = 1.0; ijData_1c(0, 0, 1) = 0.0;
875  ijData_1c(0, 1, 0) = 1.0; ijData_1c(0, 1, 1) = 0.0;
876  // C=1 contains j
877  ijData_1c(1, 0, 0) = 0.0; ijData_1c(1, 0, 1) = 1.0;
878  ijData_1c(1, 1, 0) = 0.0; ijData_1c(1, 1, 1) = 1.0;
879 
880 
881  Kokkos::View<double****> ijFields_1c("ijFields_1c", 2, 2, 2, 2);
882  // C=0, F=0 is i
883  ijFields_1c(0, 0, 0, 0) = 1.0; ijFields_1c(0, 0, 0, 1) = 0.0;
884  ijFields_1c(0, 0, 1, 0) = 1.0; ijFields_1c(0, 0, 1, 1) = 0.0;
885  // C=0, F=1 is j
886  ijFields_1c(0, 1, 0, 0) = 0.0; ijFields_1c(0, 1, 0, 1) = 1.0;
887  ijFields_1c(0, 1, 1, 0) = 0.0; ijFields_1c(0, 1, 1, 1) = 1.0;
888 
889  // C=1, F=0 is i
890  ijFields_1c(1, 0, 0, 0) = 1.0; ijFields_1c(1, 0, 0, 1) = 0.0;
891  ijFields_1c(1, 0, 1, 0) = 1.0; ijFields_1c(1, 0, 1, 1) = 0.0;
892  // C=1, F=1 is j
893  ijFields_1c(1, 1, 0, 0) = 0.0; ijFields_1c(1, 1, 0, 1) = 1.0;
894  ijFields_1c(1, 1, 1, 0) = 0.0; ijFields_1c(1, 1, 1, 1) = 1.0;
895 
896  // Output array is (C,F,P)
897  Kokkos::View<double***>outFields3("outFields",2,2,2);
898  //Kokkos::realloc(outFields, 2, 2, 2);
899  art::crossProductDataField<double>(outFields3, ijData_1c, ijFields_1c);
900 
901  if( !(outFields3(0,0,0)==0.0 && outFields3(0,0,1)==0.0 ) ) {
902  *outStream << "\n\nINCORRECT crossProductDataField (19): i x i != 0; ";
903  errorFlag++;
904  }
905  if( !(outFields3(0,1,0)==1.0 && outFields3(0,1,1)==1.0 ) ) {
906  *outStream << "\n\nINCORRECT crossProductDataField (20): i x j != 1; ";
907  errorFlag++;
908  }
909 
910  if( !(outFields3(1,0,0)==-1.0 && outFields3(1,0,1)==-1.0 ) ) {
911  *outStream << "\n\nINCORRECT crossProductDataField (21): j x i != -1; ";
912  errorFlag++;
913  }
914  if( !(outFields3(1,1,0)==0.0 && outFields3(1,1,1)==0.0 ) ) {
915  *outStream << "\n\nINCORRECT crossProductDataField (22): j x j != 0; ";
916  errorFlag++;
917  }
918 
919  *outStream \
920  << "\n"
921  << "===============================================================================\n"\
922  << "| TEST 1.d: 2D crossProductDataField operations: (C,P,D) and (F,P,D) |\n"\
923  << "===============================================================================\n";
924  /*
925  * ijData_1c ijFields_1d expected result in (C,F,P) array
926  * i,i i,j 0, 1
927  * j,j -1, 0
928  */
929  // inputFields is (F,P,D)
930  Kokkos::View<double***> ijFields_1d("ijFields_1d", 1, 2, 2);
931  // F=0 at 2 points is i,j
932  ijFields_1d(0, 0, 0) = 1.0; ijFields_1d(0, 0, 1) = 0.0;
933  ijFields_1d(0, 1, 0) = 0.0; ijFields_1d(0, 1, 1) = 1.0;
934 
935  // Output array is (C,F,P)
936  Kokkos::realloc(outFields3, 2, 1, 2);
937  art::crossProductDataField<double>(outFields3, ijData_1c, ijFields_1d);
938 
939  if( !(outFields3(0,0,0)==0.0 ) ) {
940  *outStream << "\n\nINCORRECT crossProductDataField (23): i x i != 0; ";
941  errorFlag++;
942  }
943  if( !(outFields3(0,0,1)==1.0 ) ) {
944  *outStream << "\n\nINCORRECT crossProductDataField (24): i x j != 1; ";
945  errorFlag++;
946  }
947  if( !(outFields3(1,0,0)==-1.0 ) ) {
948  *outStream << "\n\nINCORRECT crossProductDataField (25): j x i != -1; ";
949  errorFlag++;
950  }
951  if( !(outFields3(1,0,1)==0.0 ) ) {
952  *outStream << "\n\nINCORRECT crossProductDataField (26): j x j != 0; ";
953  errorFlag++;
954  }
955 
956 
957  *outStream \
958  << "\n"
959  << "===============================================================================\n"\
960  << "| TEST 2.a: 3D crossProductDataData operations: (C,P,D) and (C,P,D) |\n"\
961  << "===============================================================================\n";
962  /*
963  * ijkData_1a jkiData_2a kijData_2a
964  * i,i j,j k,k
965  * j,j k,k i,i
966  * k,k i,i j,j
967  */
968  Kokkos::View<double***> jkiData_2a("jkiData_2a", 3, 2, 3);
969  // C=0 contains j
970  jkiData_2a(0, 0, 0) = 0.0; jkiData_2a(0, 0, 1) = 1.0; jkiData_2a(0, 0, 2) = 0.0;
971  jkiData_2a(0, 1, 0) = 0.0; jkiData_2a(0, 1, 1) = 1.0; jkiData_2a(0, 1, 2) = 0.0;
972  // C=1 contains k
973  jkiData_2a(1, 0, 0) = 0.0; jkiData_2a(1, 0, 1) = 0.0; jkiData_2a(1, 0, 2) = 1.0;
974  jkiData_2a(1, 1, 0) = 0.0; jkiData_2a(1, 1, 1) = 0.0; jkiData_2a(1, 1, 2) = 1.0;
975  // C=2 contains i
976  jkiData_2a(2, 0, 0) = 1.0; jkiData_2a(2, 0, 1) = 0.0; jkiData_2a(2, 0, 2) = 0.0;
977  jkiData_2a(2, 1, 0) = 1.0; jkiData_2a(2, 1, 1) = 0.0; jkiData_2a(2, 1, 2) = 0.0;
978 
979  Kokkos::View<double***> kijData_2a("kijData_2a", 3, 2, 3);
980  // C=0 contains k
981  kijData_2a(0, 0, 0) = 0.0; kijData_2a(0, 0, 1) = 0.0; kijData_2a(0, 0, 2) = 1.0;
982  kijData_2a(0, 1, 0) = 0.0; kijData_2a(0, 1, 1) = 0.0; kijData_2a(0, 1, 2) = 1.0;
983  // C=1 contains i
984  kijData_2a(1, 0, 0) = 1.0; kijData_2a(1, 0, 1) = 0.0; kijData_2a(1, 0, 2) = 0.0;
985  kijData_2a(1, 1, 0) = 1.0; kijData_2a(1, 1, 1) = 0.0; kijData_2a(1, 1, 2) = 0.0;
986  // C=2 contains j
987  kijData_2a(2, 0, 0) = 0.0; kijData_2a(2, 0, 1) = 1.0; kijData_2a(2, 0, 2) = 0.0;
988  kijData_2a(2, 1, 0) = 0.0; kijData_2a(2, 1, 1) = 1.0; kijData_2a(2, 1, 2) = 0.0;
989 
990 
991  // ijkData_1a x ijkData_1a: outData should contain ixi=0, jxj=0, kxk=0
992  Kokkos::View<double***> outData3("outData", 3,2,3);
993  art::crossProductDataData<double>(outData3, ijkData_1a, ijkData_1a);
994 
995  for(unsigned int i = 0; i < outData3.dimension(0); i++)
996  for(unsigned int j = 0; j < outData3.dimension(1); j++)
997  for(unsigned int k = 0; k < outData3.dimension(2); k++){
998  if(outData3(i,j,k) != 0) {
999  *outStream << "\n\nINCORRECT crossProductDataData (1): i x i, j x j, or k x k != 0; ";
1000  errorFlag++;
1001  }
1002  }
1003 
1004 
1005  // ijkData_1a x jkiData_2a
1006  art::crossProductDataData<double>(outData3, ijkData_1a, jkiData_2a);
1007 
1008  // cell 0 should contain i x j = k
1009  if( !( outData3(0,0,0)==0.0 && outData3(0,0,1)==0.0 && outData3(0,0,2)==1.0 &&
1010  outData3(0,1,0)==0.0 && outData3(0,1,1)==0.0 && outData3(0,1,2)==1.0) ) {
1011  *outStream << "\n\nINCORRECT crossProductDataData (2): i x j != k; ";
1012  errorFlag++;
1013  }
1014 
1015  // cell 1 should contain j x k = i
1016  if( !( outData3(1,0,0)==1.0 && outData3(1,0,1)==0.0 && outData3(1,0,2)==0.0 &&
1017  outData3(1,1,0)==1.0 && outData3(1,1,1)==0.0 && outData3(1,1,2)==0.0) ) {
1018  *outStream << "\n\nINCORRECT crossProductDataData (3): j x k != i; ";
1019  errorFlag++;
1020  }
1021 
1022  // cell 2 should contain k x i = j
1023  if( !( outData3(2,0,0)==0.0 && outData3(2,0,1)==1.0 && outData3(2,0,2)==0.0 &&
1024  outData3(2,1,0)==0.0 && outData3(2,1,1)==1.0 && outData3(2,1,2)==0.0) ) {
1025  *outStream << "\n\nINCORRECT crossProductDataData (4): k x i != j; ";
1026  errorFlag++;
1027  }
1028 
1029 
1030  // ijkData_1a x kijData_2a
1031  art::crossProductDataData<double>(outData3, ijkData_1a, kijData_2a);
1032 
1033  // cell 0 should contain i x k = -j
1034  if( !( outData3(0,0,0)==0.0 && outData3(0,0,1)==-1.0 && outData3(0,0,2)==0.0 &&
1035  outData3(0,1,0)==0.0 && outData3(0,1,1)==-1.0 && outData3(0,1,2)==0.0) ) {
1036  *outStream << "\n\nINCORRECT crossProductDataData (5): i x k != -j; ";
1037  errorFlag++;
1038  }
1039 
1040  // cell 1 should contain j x i = -k
1041  if( !( outData3(1,0,0)==0.0 && outData3(1,0,1)==0.0 && outData3(1,0,2)==-1.0 &&
1042  outData3(1,1,0)==0.0 && outData3(1,1,1)==0.0 && outData3(1,1,2)==-1.0) ) {
1043  *outStream << "\n\nINCORRECT crossProductDataData (6): j x i != -k; ";
1044  errorFlag++;
1045  }
1046 
1047  // cell 2 should contain k x j = -i
1048  if( !( outData3(2,0,0)==-1.0 && outData3(2,0,1)==0.0 && outData3(2,0,2)==0.0 &&
1049  outData3(2,1,0)==-1.0 && outData3(2,1,1)==0.0 && outData3(2,1,2)==0.0) ) {
1050  *outStream << "\n\nINCORRECT crossProductDataData (7): k x j != -i; ";
1051  errorFlag++;
1052  }
1053 
1054 
1055  *outStream \
1056  << "\n"
1057  << "===============================================================================\n"\
1058  << "| TEST 2.b: 3D crossProductDataData operations: (C,P,D) and (P,D) |\n"\
1059  << "===============================================================================\n";
1060  /*
1061  * ijkData_1b ijkData_2b expected result in (C,P,D) array
1062  * i,i,i i,j,k 0, k,-j
1063  * j,j,j -k, 0, i
1064  * k,k,k j,-i, 0
1065  */
1066  // input data is (P,D)
1067  Kokkos::View<double**> ijkData_2b("ijkData_2b", 3, 3);
1068  // F=0 at 3 points is (i,j,k)
1069  ijkData_2b(0, 0) = 1.0; ijkData_2b(0, 1) = 0.0; ijkData_2b(0, 2) = 0.0;
1070  ijkData_2b(1, 0) = 0.0; ijkData_2b(1, 1) = 1.0; ijkData_2b(1, 2) = 0.0;
1071  ijkData_2b(2, 0) = 0.0; ijkData_2b(2, 1) = 0.0; ijkData_2b(2, 2) = 1.0;
1072 
1073  // Output array is (C,P,D)
1074  Kokkos::realloc(outData3, 3, 3, 3);
1075  art::crossProductDataData<double>(outData3, ijkData_1b, ijkData_2b);
1076 
1077  // checks for C = 0
1078  if( !(outData3(0,0,0)==0.0 && outData3(0,0,1)==0.0 && outData3(0,0,2)==0.0) ) {
1079  *outStream << "\n\nINCORRECT crossProductDataData (5): i x i != 0; ";
1080  errorFlag++;
1081  }
1082  if( !(outData3(0,1,0)==0.0 && outData3(0,1,1)==0.0 && outData3(0,1,2)==1.0) ) {
1083  *outStream << "\n\nINCORRECT crossProductDataData (6): i x j != k; ";
1084  errorFlag++;
1085  }
1086  if( !(outData3(0,2,0)==0.0 && outData3(0,2,1)==-1.0 && outData3(0,2,2)==0.0) ) {
1087  *outStream << "\n\nINCORRECT crossProductDataData (7): i x k != -j; ";
1088  errorFlag++;
1089  }
1090 
1091  // checks for C = 1
1092  if( !(outData3(1,0,0)==0.0 && outData3(1,0,1)==0.0 && outData3(1,0,2)==-1.0) ) {
1093  *outStream << "\n\nINCORRECT crossProductDataData (8): j x i != -k; ";
1094  errorFlag++;
1095  }
1096  if( !(outData3(1,1,0)==0.0 && outData3(1,1,1)==0.0 && outData3(1,1,2)==0.0) ) {
1097  *outStream << "\n\nINCORRECT crossProductDataData (9): j x j != 0; ";
1098  errorFlag++;
1099  }
1100  if( !(outData3(1,2,0)==1.0 && outData3(1,2,1)==0.0 && outData3(1,2,2)==0.0) ) {
1101  *outStream << "\n\nINCORRECT crossProductDataData (10): j x k != i; ";
1102  errorFlag++;
1103  }
1104 
1105  // checks for C = 2
1106  if( !(outData3(2,0,0)==0.0 && outData3(2,0,1)==1.0 && outData3(2,0,2)==0.0) ) {
1107  *outStream << "\n\nINCORRECT crossProductDataData (11): k x i != j; ";
1108  errorFlag++;
1109  }
1110  if( !(outData3(2,1,0)==-1.0 && outData3(2,1,1)==0.0 && outData3(2,1,2)==0.0) ) {
1111  *outStream << "\n\nINCORRECT crossProductDataData (12): k x j != -i; ";
1112  errorFlag++;
1113  }
1114  if( !(outData3(2,2,0)==0.0 && outData3(2,2,1)==0.0 && outData3(2,2,2)==0.0) ) {
1115  *outStream << "\n\nINCORRECT crossProductDataData (13): k x k != 0; ";
1116  errorFlag++;
1117  }
1118 
1119 
1120  *outStream \
1121  << "\n"
1122  << "===============================================================================\n"\
1123  << "| TEST 2.c: 2D crossProductDataData operations: (C,P,D) and (C,P,D) |\n"\
1124  << "===============================================================================\n";
1125  /*
1126  * ijData_1c jiData_2c
1127  * i,i j,j
1128  * j,j i,i
1129  */
1130  Kokkos::View<double***> jiData_2c("jiData_2c", 2, 2, 2);
1131  // C=0 contains j
1132  jiData_2c(0, 0, 0) = 0.0; jiData_2c(0, 0, 1) = 1.0;
1133  jiData_2c(0, 1, 0) = 0.0; jiData_2c(0, 1, 1) = 1.0;
1134  // C=1 contains i
1135  jiData_2c(1, 0, 0) = 1.0; jiData_2c(1, 0, 1) = 0.0;
1136  jiData_2c(1, 1, 0) = 1.0; jiData_2c(1, 1, 1) = 0.0;
1137 
1138 
1139  // ijData_1c x ijData_1c: outData should contain ixi=0, jxj=0
1140  Kokkos::View<double**> outData2("outputData2",2,2);
1141  art::crossProductDataData<double>(outData2, ijData_1c, ijData_1c);
1142 
1143  for(unsigned int i = 0; i < outData2.dimension(0); i++){
1144  for(unsigned int j = 0; j < outData2.dimension(1); j++){
1145  if(outData2(i,j) != 0) {
1146  *outStream << "\n\nINCORRECT crossProductDataData (14): i x i or j x j != 0; ";
1147  errorFlag++;
1148  }
1149  }
1150  }
1151  // ijData_1c x jiData_1c: outData should contain ixi=0, jxj=0
1152  art::crossProductDataData<double>(outData2, ijData_1c, jiData_2c);
1153 
1154  if( !(outData2(0,0)==1.0 && outData2(0,1)==1.0 ) ) {
1155  *outStream << "\n\nINCORRECT crossProductDataData (15): i x j != 1; ";
1156  errorFlag++;
1157  }
1158  if( !(outData2(1,0)==-1.0 && outData2(1,1)==-1.0 ) ) {
1159  *outStream << "\n\nINCORRECT crossProductDataData (16): j x i != -1; ";
1160  errorFlag++;
1161  }
1162 
1163  *outStream \
1164  << "\n"
1165  << "===============================================================================\n"\
1166  << "| TEST 2.d: 2D crossProductDataData operations: (C,P,D) and (P,D) |\n"\
1167  << "===============================================================================\n";
1168  /*
1169  * ijData_1c ijData_2d expected result in (C,P) array
1170  * i,i i,j 0, 1
1171  * j,j -1, 0
1172  */
1173  Kokkos::View<double**> ijData_2d("ijData_2d", 2, 2);
1174  ijData_2d(0, 0) = 1.0; ijData_2d(0, 1) = 0.0;
1175  ijData_2d(1, 0) = 0.0; ijData_2d(1, 1) = 1.0;
1176 
1177  Kokkos::realloc(outData2,2,2);
1178  art::crossProductDataData<double>(outData2, ijData_1c, ijData_2d);
1179 
1180  if( !(outData2(0,0)==0.0 ) ) {
1181  *outStream << "\n\nINCORRECT crossProductDataData (17): i x i != 0; ";
1182  errorFlag++;
1183  }
1184  if( !(outData2(0,1)==1.0 ) ) {
1185  *outStream << "\n\nINCORRECT crossProductDataData (18): i x j != 1; ";
1186  errorFlag++;
1187  }
1188  if( !(outData2(1,0)==-1.0 ) ) {
1189  *outStream << "\n\nINCORRECT crossProductDataData (19): j x i != -1; ";
1190  errorFlag++;
1191  }
1192  if( !(outData2(1,1)==0.0 ) ) {
1193  *outStream << "\n\nINCORRECT crossProductDataData (20): j x j != 0; ";
1194  errorFlag++;
1195  }
1196 
1197 
1198  *outStream \
1199  << "\n"
1200  << "===============================================================================\n"\
1201  << "| TEST 3.a: outerProductDataField operations: (C,P,D) and (C,F,P,D) |\n"\
1202  << "===============================================================================\n";
1203  /*
1204  * ijkData_1a ijkFields_1a Expected result in (C,F,P,D,D) array:
1205  * i,i i,i j,j, k,k (0,0) (0,1) (0,2)
1206  * j,j i,i j,j, k,k (1,0) (1,1) (1,2)
1207  * k,k i,i j,j, k,k (2,0) (2,1) (2,2)
1208  * Indicates the only non-zero element of (C,F,P,*,*), specifically,
1209  * element with row = cell and col = field should equal 1; all other should equal 0
1210  */
1211 
1212  Kokkos::View<double*****> outFields5("outFields5",3, 3, 2, 3, 3);
1213  art::outerProductDataField<double>(outFields5, ijkData_1a, ijkFields_1a);
1214 
1215  for(unsigned int cell = 0; cell < ijkData_1a.dimension(0); cell++){
1216  for(unsigned int field = 0; field < ijkFields_1a.dimension(1); field++){
1217  for(unsigned int point = 0; point < ijkData_1a.dimension(1); point++){
1218  for(unsigned int row = 0; row < ijkData_1a.dimension(2); row++){
1219  for(unsigned int col = 0; col < ijkData_1a.dimension(2); col++){
1220 
1221  // element with row = cell and col = field should equal 1; all other should equal 0
1222  if( (row == cell && col == field) ){
1223  if(outFields5(cell, field, point, row, col) != 1.0) {
1224  *outStream << "\n\nINCORRECT outerProductDataField (1): computed value is "
1225  << outFields5(cell, field, point, row, col) << " whereas correct value is 1.0";
1226  errorFlag++;
1227  }
1228  }
1229  else {
1230  if(outFields5(cell, field, point, row, col) != 0.0) {
1231  *outStream << "\n\nINCORRECT outerProductDataField (2): computed value is "
1232  << outFields5(cell, field, point, row, col) << " whereas correct value is 0.0";
1233  errorFlag++;
1234  }
1235  } // if
1236  }// col
1237  }// row
1238  }// point
1239  }// field
1240  }// cell
1241 
1242  *outStream \
1243  << "\n"
1244  << "===============================================================================\n"\
1245  << "| TEST 3.b: outerProductDataField operations: (C,P,D) and (F,P,D) |\n"\
1246  << "===============================================================================\n";
1247  /*
1248  * ijkData_1b ijkFields_1b expected result in (C,F,P,D,D) array
1249  * i,i,i i,j,k (0,0) (0,1) (0,2)
1250  * j,j,j (1,0) (1,1) (1,2)
1251  * k,k,k (2,0) (2,1) (2,2)
1252  * Indicates the only non-zero element of (C,F,P,*,*), specifically,
1253  * element with row = cell and col = point should equal 1; all other should equal 0
1254  */
1255 
1256  Kokkos::realloc(outFields5,3, 1, 3, 3, 3);
1257  art::outerProductDataField<double>(outFields5, ijkData_1b, ijkFields_1b);
1258 
1259  for(unsigned int cell = 0; cell < ijkData_1b.dimension(0); cell++){
1260  for(unsigned int field = 0; field < ijkFields_1b.dimension(0); field++){
1261  for(unsigned int point = 0; point < ijkData_1b.dimension(1); point++){
1262  for(unsigned int row = 0; row < ijkData_1b.dimension(2); row++){
1263  for(unsigned int col = 0; col < ijkData_1b.dimension(2); col++){
1264 
1265  // element with row = cell and col = point should equal 1; all other should equal 0
1266  if( (row == cell && col == point) ){
1267  if(outFields5(cell, field, point, row, col) != 1.0) {
1268  *outStream << "\n\nINCORRECT outerProductDataField (3): computed value is "
1269  << outFields5(cell, field, point, row, col) << " whereas correct value is 1.0";
1270  errorFlag++;
1271 
1272  }
1273  }
1274  else {
1275  if(outFields5(cell, field, point, row, col) != 0.0) {
1276  *outStream << "\n\nINCORRECT outerProductDataField (4): computed value is "
1277  << outFields5(cell, field, point, row, col) << " whereas correct value is 0.0";
1278  errorFlag++;
1279  }
1280  } // if
1281  }// col
1282  }// row
1283  }// point
1284  }// field
1285  }// cell
1286  *outStream \
1287  << "\n"
1288  << "===============================================================================\n"\
1289  << "| TEST 4.a: outerProductDataData operations: (C,P,D) and (C,P,D) |\n"\
1290  << "===============================================================================\n";
1291  /*
1292  * ijkData_1a jkiData_2a kijData_2a
1293  * i,i j,j k,k
1294  * j,j k,k i,i
1295  * k,k i,i j,j
1296  *
1297  * Expected results are stated with each test case.
1298  */
1299  Kokkos::View<double****>outData4("outData4",3, 2, 3, 3);
1300  art::outerProductDataData<double>(outData4, ijkData_1a, ijkData_1a);
1301  for(unsigned int cell = 0; cell < ijkData_1a.dimension(0); cell++){
1302  for(unsigned int point = 0; point < ijkData_1a.dimension(1); point++){
1303  for(unsigned int row = 0; row < ijkData_1a.dimension(2); row++){
1304  for(unsigned int col = 0; col < ijkData_1a.dimension(2); col++){
1305 
1306  // element with row = cell and col = cell should equal 1; all other should equal 0
1307  if( (row == cell && col == cell) ){
1308  if(outData4(cell, point, row, col) != 1.0) {
1309  *outStream << "\n\nINCORRECT outerProductDataData (1): computed value is "
1310  << outData4(cell, point, row, col) << " whereas correct value is 1.0";
1311  errorFlag++;
1312  }
1313  }
1314  else {
1315  if(outData4(cell, point, row, col) != 0.0) {
1316  *outStream << "\n\nINCORRECT outerProductDataData (2): computed value is "
1317  << outData4(cell, point, row, col) << " whereas correct value is 0.0";
1318  errorFlag++;
1319  }
1320  } // if
1321  }// col
1322  }// row
1323  }// point
1324  }// cell
1325 
1326  Kokkos::deep_copy(outData4,0.0);
1327  art::outerProductDataData<double>(outData4, ijkData_1a, jkiData_2a);
1328  for(unsigned int cell = 0; cell < ijkData_1a.dimension(0); cell++){
1329  for(unsigned int point = 0; point < ijkData_1a.dimension(1); point++){
1330  for(unsigned int row = 0; row < ijkData_1a.dimension(2); row++){
1331  for(unsigned int col = 0; col < ijkData_1a.dimension(2); col++){
1332 
1333  // element with row = cell and col = cell + 1 (mod 3) should equal 1; all other should equal 0
1334  if( (row == cell && col == (cell + 1) % 3) ){
1335  if(outData4(cell, point, row, col) != 1.0) {
1336  *outStream << "\n\nINCORRECT outerProductDataData (3): computed value is "
1337  << outData4(cell, point, row, col) << " whereas correct value is 1.0";
1338  errorFlag++;
1339  }
1340  }
1341  else {
1342  if(outData4(cell, point, row, col) != 0.0) {
1343  *outStream << "\n\nINCORRECT outerProductDataData (4): computed value is "
1344  << outData4(cell, point, row, col) << " whereas correct value is 0.0";
1345  errorFlag++;
1346  }
1347  } // if
1348  }// col
1349  }// row
1350  }// point
1351  }// cell
1352 
1353 
1354  Kokkos::deep_copy(outData4,0.0);
1355  art::outerProductDataData<double>(outData4, ijkData_1a, kijData_2a);
1356  for(unsigned int cell = 0; cell < ijkData_1a.dimension(0); cell++){
1357  for(unsigned int point = 0; point < ijkData_1a.dimension(1); point++){
1358  for(unsigned int row = 0; row < ijkData_1a.dimension(2); row++){
1359  for(unsigned int col = 0; col < ijkData_1a.dimension(2); col++){
1360 
1361  // element with row = cell and col = cell + 2 (mod 3) should equal 1; all other should equal 0
1362  if( (row == cell && col == (cell + 2) % 3) ){
1363  if(outData4(cell, point, row, col) != 1.0) {
1364  *outStream << "\n\nINCORRECT outerProductDataData (5): computed value is "
1365  << outData4(cell, point, row, col) << " whereas correct value is 1.0";
1366  errorFlag++;
1367  }
1368  }
1369  else {
1370  if(outData4(cell, point, row, col) != 0.0) {
1371  *outStream << "\n\nINCORRECT outerProductDataData (6): computed value is "
1372  << outData4(cell, point, row, col) << " whereas correct value is 0.0";
1373  errorFlag++;
1374  }
1375  } // if
1376  }// col
1377  }// row
1378  }// point
1379  }// cell
1380 
1381 
1382  *outStream \
1383  << "\n"
1384  << "===============================================================================\n"\
1385  << "| TEST 4.b: outerProductDataData operations: (C,P,D) and (P,D) |\n"\
1386  << "===============================================================================\n";
1387  /*
1388  * ijkData_1b ijkData_2b expected result in (C,P,D,D) array
1389  * i,i,i i,j,k (0,0) (0,1) (0,2)
1390  * j,j,j (1,0) (1,1) (1,2)
1391  * k,k,k (2,0) (2,1) (2,2)
1392  * Indicates the only non-zero element of (C,P,*,*), specifically,
1393  * element with row = cell and col = point should equal 1; all other should equal 0
1394  */
1395  Kokkos::realloc(outData4,3, 3, 3, 3);
1396  art::outerProductDataData<double>(outData4, ijkData_1b, ijkData_2b);
1397  for(unsigned int cell = 0; cell < ijkData_1b.dimension(0); cell++){
1398  for(unsigned int point = 0; point < ijkData_1b.dimension(1); point++){
1399  for(unsigned int row = 0; row < ijkData_1b.dimension(2); row++){
1400  for(unsigned int col = 0; col < ijkData_1b.dimension(2); col++){
1401 
1402  // element with row = cell and col = cell + 2 (mod 3) should equal 1; all other should equal 0
1403  if( (row == cell && col == point) ){
1404  if(outData4(cell, point, row, col) != 1.0) {
1405  *outStream << "\n\nINCORRECT outerProductDataData (7): computed value is "
1406  << outData4(cell, point, row, col) << " whereas correct value is 1.0";
1407  errorFlag++;
1408  }
1409  }
1410  else {
1411  if(outData4(cell, point, row, col) != 0.0) {
1412  *outStream << "\n\nINCORRECT outerProductDataData (8): computed value is "
1413  << outData4(cell, point, row, col) << " whereas correct value is 0.0";
1414  errorFlag++;
1415  }
1416  } // if
1417  }// col
1418  }// row
1419  }// point
1420  }// cell
1421 
1422  *outStream \
1423  << "\n"
1424  << "===============================================================================\n"\
1425  << "| TEST 5.a: matvecProductDataField operations: (C,P,D,D) and (C,F,P,D) |\n"\
1426  << "===============================================================================\n";
1427  /*
1428  * inputMat inputVecFields outFields
1429  * 1 1 1 0 0 0 0 1 -1 -1 0 3 0 0
1430  * -1 2 -1 -1 -2 -3 0 1 -1 1 0 0 6 2
1431  * 1 2 3 -2 6 -4 0 1 -1 -1 0 6 0 12
1432  */
1433 
1434  // (C,P,D,D)
1435  Kokkos::View<double****> inputMat("inputMat",2,1,3,3);
1436  // cell 0
1437  inputMat(0,0,0,0) = 1.0; inputMat(0,0,0,1) = 1.0; inputMat(0,0,0,2) = 1.0;
1438  inputMat(0,0,1,0) =-1.0; inputMat(0,0,1,1) = 2.0; inputMat(0,0,1,2) =-1.0;
1439  inputMat(0,0,2,0) = 1.0; inputMat(0,0,2,1) = 2.0; inputMat(0,0,2,2) = 3.0;
1440  // cell 1
1441  inputMat(1,0,0,0) = 0.0; inputMat(1,0,0,1) = 0.0; inputMat(1,0,0,2) = 0.0;
1442  inputMat(1,0,1,0) =-1.0; inputMat(1,0,1,1) =-2.0; inputMat(1,0,1,2) =-3.0;
1443  inputMat(1,0,2,0) =-2.0; inputMat(1,0,2,1) = 6.0; inputMat(1,0,2,2) =-4.0;
1444 
1445  // (C,F,P,D)
1446  Kokkos::View<double****> inputVecFields4("inputVecFields4",2,2,1,3);
1447  // cell 0; fields 0,1
1448  inputVecFields4(0,0,0,0) = 0.0; inputVecFields4(0,0,0,1) = 0.0; inputVecFields4(0,0,0,2) = 0.0;
1449  inputVecFields4(0,1,0,0) = 1.0; inputVecFields4(0,1,0,1) = 1.0; inputVecFields4(0,1,0,2) = 1.0;
1450  // cell 1; fields 0,1
1451  inputVecFields4(1,0,0,0) =-1.0; inputVecFields4(1,0,0,1) =-1.0; inputVecFields4(1,0,0,2) =-1.0;
1452  inputVecFields4(1,1,0,0) =-1.0; inputVecFields4(1,1,0,1) = 1.0; inputVecFields4(1,1,0,2) =-1.0;
1453 
1454  // (C,F,P,D) - true
1455  Kokkos::View<double****> outFieldsCorrect("outFieldsCorrect",2,2,1,3);
1456  // cell 0; fields 0,1
1457  outFieldsCorrect(0,0,0,0) = 0.0; outFieldsCorrect(0,0,0,1) = 0.0; outFieldsCorrect(0,0,0,2) = 0.0;
1458  outFieldsCorrect(0,1,0,0) = 3.0; outFieldsCorrect(0,1,0,1) = 0.0; outFieldsCorrect(0,1,0,2) = 6.0;
1459  // cell 1; fields 0,1
1460  outFieldsCorrect(1,0,0,0) = 0.0; outFieldsCorrect(1,0,0,1) = 6.0; outFieldsCorrect(1,0,0,2) = 0.0;
1461  outFieldsCorrect(1,1,0,0) = 0.0; outFieldsCorrect(1,1,0,1) = 2.0; outFieldsCorrect(1,1,0,2) = 12.0;
1462 
1463  // (C,F,P,D)
1464  Kokkos::realloc(outFields4,2,2,1,3);
1465  art::matvecProductDataField<double>(outFields4, inputMat, inputVecFields4);
1466 
1467  // test loop
1468  for(unsigned int cell = 0; cell < outFields4.dimension(0); cell++){
1469  for(unsigned int field = 0; field < outFields4.dimension(1); field++){
1470  for(unsigned int point = 0; point < outFields4.dimension(2); point++){
1471  for(unsigned int row = 0; row < outFields4.dimension(3); row++){
1472  if(outFields4(cell, field, point, row) != outFieldsCorrect(cell, field, point, row)) {
1473  *outStream << "\n\nINCORRECT matvecProductDataField (1): \n value at multi-index ("
1474  << cell << "," << field << "," << point << "," << row << ") = "
1475  << outFields4(cell, field, point, row) << " but correct value is "
1476  << outFieldsCorrect(cell, field, point, row) <<"\n";
1477  errorFlag++;
1478  }
1479  }//row
1480  }// point
1481  }// field
1482  }// cell
1483 
1484 
1485  *outStream \
1486  << "\n"
1487  << "===============================================================================\n"\
1488  << "| TEST 5.b: matvecProductDataField operations: (C,P,D,D) and (F,P,D) |\n"\
1489  << "===============================================================================\n";
1490  /*
1491  * inputMat inputVecFields outFields
1492  * 1 1 1 0 0 0 0 1 -1 -1 0 3 -3 -1 0 0 0 0
1493  * -1 2 -1 -1 -2 -3 0 1 -1 1 0 0 0 4 0 -6 6 2
1494  * 1 2 3 -2 6 -4 0 1 -1 -1 0 6 -6 -2 0 0 0 12
1495  * Use the same 4 vector fields as above but formatted as (F,P,D) array
1496  */
1497  // (C,F,P,D)
1498  Kokkos::View<double***> inputVecFields3("inputVecFields3",4,1,3);
1499 
1500  // fields 0,1,2,3
1501  inputVecFields3(0,0,0) = 0.0; inputVecFields3(0,0,1) = 0.0; inputVecFields3(0,0,2) = 0.0;
1502  inputVecFields3(1,0,0) = 1.0; inputVecFields3(1,0,1) = 1.0; inputVecFields3(1,0,2) = 1.0;
1503  inputVecFields3(2,0,0) =-1.0; inputVecFields3(2,0,1) =-1.0; inputVecFields3(2,0,2) =-1.0;
1504  inputVecFields3(3,0,0) =-1.0; inputVecFields3(3,0,1) = 1.0; inputVecFields3(3,0,2) =-1.0;
1505 
1506  // (C,F,P,D) - true
1507  Kokkos::realloc(outFieldsCorrect,2,4,1,3);
1508  // cell 0; fields 0,1,2,3
1509  outFieldsCorrect(0,0,0,0) = 0.0; outFieldsCorrect(0,0,0,1) = 0.0; outFieldsCorrect(0,0,0,2) = 0.0;
1510  outFieldsCorrect(0,1,0,0) = 3.0; outFieldsCorrect(0,1,0,1) = 0.0; outFieldsCorrect(0,1,0,2) = 6.0;
1511  outFieldsCorrect(0,2,0,0) =-3.0; outFieldsCorrect(0,2,0,1) = 0.0; outFieldsCorrect(0,2,0,2) =-6.0;
1512  outFieldsCorrect(0,3,0,0) =-1.0; outFieldsCorrect(0,3,0,1) = 4.0; outFieldsCorrect(0,3,0,2) =-2.0;
1513  // cell 1; fields 0,1,2,3
1514  outFieldsCorrect(1,0,0,0) = 0.0; outFieldsCorrect(1,0,0,1) = 0.0; outFieldsCorrect(1,0,0,2) = 0.0;
1515  outFieldsCorrect(1,1,0,0) = 0.0; outFieldsCorrect(1,1,0,1) =-6.0; outFieldsCorrect(1,1,0,2) = 0.0;
1516  outFieldsCorrect(1,2,0,0) = 0.0; outFieldsCorrect(1,2,0,1) = 6.0; outFieldsCorrect(1,2,0,2) = 0.0;
1517  outFieldsCorrect(1,3,0,0) = 0.0; outFieldsCorrect(1,3,0,1) = 2.0; outFieldsCorrect(1,3,0,2) =12.0;
1518 
1519  // (C,F,P,D)
1520 
1521  Kokkos::realloc(outFields4,2,4,1,3);
1522  art::matvecProductDataField<double>(outFields4, inputMat, inputVecFields3);
1523 
1524  // test loop
1525  for(unsigned int cell = 0; cell < outFields4.dimension(0); cell++){
1526  for(unsigned int field = 0; field < outFields4.dimension(1); field++){
1527  for(unsigned int point = 0; point < outFields4.dimension(2); point++){
1528  for(unsigned int row = 0; row < outFields4.dimension(3); row++){
1529  if(outFields4(cell, field, point, row) != outFieldsCorrect(cell, field, point, row)) {
1530  *outStream << "\n\nINCORRECT matvecProductDataField (2): \n value at multi-index ("
1531  << cell << "," << field << "," << point << "," << row << ") = "
1532  << outFields4(cell, field, point, row) << " but correct value is "
1533  << outFieldsCorrect(cell, field, point, row) <<"\n";
1534  errorFlag++;
1535  }
1536  }//row
1537  }// point
1538  }// field
1539  }// cell
1540 
1541 
1542  *outStream \
1543  << "\n"
1544  << "===============================================================================\n"\
1545  << "| TEST 5.c: matvecProductDataField random tests: branch inputFields(C,F,P,D) |\n"\
1546  << "===============================================================================\n";
1547  /*
1548  * d1 is spatial dimension and should be 1, 2 or 3. If d1>3, the RealSpaceTools function 'inverse' will fail
1549  */
1550  {// test 5.c scope
1551  int c=5, p=9, f=7, d1=3;
1552  double zero = INTREPID_TOL*10000.0;
1553 
1554  Kokkos::View<double****> in_c_f_p_d("in_c_f_p_d",c, f, p, d1);
1555  Kokkos::View<double****> out_c_f_p_d("out_c_f_p_d",c, f, p, d1);
1556  Kokkos::View<double****> outi_c_f_p_d("outi_c_f_p_d",c, f, p, d1);
1557 
1558  Kokkos::View<double**> data_c_p("data_c_p",c, p);
1559  Kokkos::View<double**> datainv_c_p("datainv_c_p",c, p);
1560  Kokkos::View<double**> data_c_1("data_c_1",c, 1);
1561  Kokkos::View<double**> datainv_c_1("datainv_c_1",c, 1);
1562  Kokkos::View<double***> data_c_p_d("data_c_p_d",c, p, d1);
1563  Kokkos::View<double***> datainv_c_p_d("datainv_c_p_d",c, p, d1);
1564  Kokkos::View<double***> data_c_1_d("data_c_1_d",c, 1, d1);
1565  Kokkos::View<double***> datainv_c_1_d("datainv_c_1_d",c, 1, d1);
1566  Kokkos::View<double****> data_c_p_d_d("data_c_p_d_d",c, p, d1, d1);
1567  Kokkos::View<double****> datainv_c_p_d_d("datainv_c_p_d_d",c, p, d1, d1);
1568  Kokkos::View<double****> datainvtrn_c_p_d_d("datainvtrn_c_p_d_d",c, p, d1, d1);
1569  Kokkos::View<double****> data_c_1_d_d("data_c_1_d_d",c, 1, d1, d1);
1570  Kokkos::View<double****> datainv_c_1_d_d("datainv_c_1_d_d",c, 1, d1, d1);
1571  Kokkos::View<double****> datainvtrn_c_1_d_d("datainvtrn_c_1_d_d",c, 1, d1, d1);
1572  /***********************************************************************************************
1573  * Constant diagonal tensor: inputData(C,P) *
1574  **********************************************************************************************/
1575  for (unsigned int i=0; i<in_c_f_p_d.dimension(0); i++)
1576  for (unsigned int j=0; j<in_c_f_p_d.dimension(1); j++)
1577  for (unsigned int k=0; k<in_c_f_p_d.dimension(2); k++)
1578  for (unsigned int l=0; l<in_c_f_p_d.dimension(3); l++)
1579  in_c_f_p_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
1580 
1581  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
1582  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
1583  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
1584  datainv_c_p(i,j) = 1.0 / data_c_p(i,j);
1585  }
1586 
1587  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
1588  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
1589  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
1590  datainv_c_1(i,j) = 1.0 / data_c_1(i,j);
1591  }
1592 
1593  // Tensor values vary by point:
1594  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p, in_c_f_p_d);
1595  art::matvecProductDataField<double>(out_c_f_p_d, datainv_c_p, out_c_f_p_d);
1596  rst::subtract(out_c_f_p_d, in_c_f_p_d);
1597  if (rst::vectorNorm(out_c_f_p_d, NORM_ONE) > zero) {
1598  *outStream << "\n\nINCORRECT matvecProductDataField (3): check scalar inverse property\n\n";
1599  errorFlag = -1000;
1600  }
1601  // Tensor values do not vary by point:
1602  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1, in_c_f_p_d);
1603  art::matvecProductDataField<double>(out_c_f_p_d, datainv_c_1, out_c_f_p_d);
1604  rst::subtract(out_c_f_p_d, in_c_f_p_d);
1605  if (rst::vectorNorm(out_c_f_p_d, NORM_ONE) > zero) {
1606  *outStream << "\n\nINCORRECT matvecProductDataField (4): check scalar inverse property\n\n";
1607  errorFlag = -1000;
1608  }
1609  /***********************************************************************************************
1610  * Non-onstant diagonal tensor: inputData(C,P,D) *
1611  **********************************************************************************************/
1612  for (unsigned int i=0; i<in_c_f_p_d.dimension(0); i++)
1613  for (unsigned int j=0; j<in_c_f_p_d.dimension(1); j++)
1614  for (unsigned int k=0; k<in_c_f_p_d.dimension(2); k++)
1615  for (unsigned int l=0; l<in_c_f_p_d.dimension(3); l++)
1616  in_c_f_p_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
1617 
1618  for (unsigned int i=0; i<data_c_p_d.dimension(0); i++)
1619  for (unsigned int j=0; j<data_c_p_d.dimension(1); j++)
1620  for (unsigned int k=0; k<data_c_p_d.dimension(2); k++){
1621  data_c_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
1622  datainv_c_p_d(i,j,k) = 1.0 / data_c_p_d(i,j,k);
1623  }
1624 
1625  for (unsigned int i=0; i<data_c_1_d.dimension(0); i++)
1626  for (unsigned int j=0; j<data_c_1_d.dimension(1); j++)
1627  for (unsigned int k=0; k<data_c_1_d.dimension(2); k++){
1628  data_c_1_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
1629  datainv_c_1_d(i,j,k) = 1.0 / data_c_1_d(i,j,k);
1630  }
1631 
1632 
1633  // Tensor values vary by point:
1634  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d, in_c_f_p_d);
1635  art::matvecProductDataField<double>(out_c_f_p_d, datainv_c_p_d, out_c_f_p_d);
1636  rst::subtract(out_c_f_p_d, in_c_f_p_d);
1637  if (rst::vectorNorm(out_c_f_p_d, NORM_ONE) > zero) {
1638  *outStream << "\n\nINCORRECT matvecProductDataField (5): check scalar inverse property\n\n";
1639  errorFlag = -1000;
1640  }
1641  // Tensor values do not vary by point
1642  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d, in_c_f_p_d);
1643  art::matvecProductDataField<double>(out_c_f_p_d, datainv_c_1_d, out_c_f_p_d);
1644  rst::subtract(out_c_f_p_d, in_c_f_p_d);
1645  if (rst::vectorNorm(out_c_f_p_d, NORM_ONE) > zero) {
1646  *outStream << "\n\nINCORRECT matvecProductDataField (6): check scalar inverse property\n\n";
1647  errorFlag = -1000;
1648  }
1649  /***********************************************************************************************
1650  * Full tensor: inputData(C,P,D,D) *
1651  **********************************************************************************************/
1652 
1653  for (unsigned int i=0; i<in_c_f_p_d.dimension(0); i++)
1654  for (unsigned int j=0; j<in_c_f_p_d.dimension(1); j++)
1655  for (unsigned int k=0; k<in_c_f_p_d.dimension(2); k++)
1656  for (unsigned int l=0; l<in_c_f_p_d.dimension(3); l++)
1657  in_c_f_p_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
1658 
1659  for (unsigned int i=0; i<data_c_p_d_d.dimension(0); i++)
1660  for (unsigned int j=0; j<data_c_p_d_d.dimension(1); j++)
1661  for (unsigned int k=0; k<data_c_p_d_d.dimension(2); k++)
1662  for (unsigned int l=0; l<data_c_p_d_d.dimension(3); l++)
1663  data_c_p_d_d(i, j, k, l) = Teuchos::ScalarTraits<double>::random();
1664 
1665  RealSpaceTools<double>::inverse(datainv_c_p_d_d, data_c_p_d_d);
1666 
1667  for (int ic=0; ic < c; ic++) {
1668  for (int ip=0; ip < p; ip++) {
1669  for (int i=0; i<d1; i++) {
1670  for (int j=0; j<d1; j++) {
1671  data_c_p_d_d(ic, ip, i, j) = Teuchos::ScalarTraits<double>::random();
1672  }
1673  }
1674  }
1675  }
1676 
1677  RealSpaceTools<double>::inverse(datainv_c_p_d_d, data_c_p_d_d);
1678 
1679  for (int ic=0; ic < c; ic++) {
1680  for (int ip=0; ip < 1; ip++) {
1681  for (int i=0; i<d1; i++) {
1682  for (int j=0; j<d1; j++) {
1683  data_c_1_d_d(ic, 0, i, j) = Teuchos::ScalarTraits<double>::random();
1684  }
1685  }
1686  }
1687  }
1688  RealSpaceTools<double>::inverse(datainv_c_1_d_d, data_c_1_d_d);
1689 
1690 
1691  // Tensor values vary by point: test "N" and "T" options (no transpose/transpose)
1692  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d_d, in_c_f_p_d);
1693  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_p_d_d, out_c_f_p_d);
1694  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
1695  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
1696  *outStream << "\n\nINCORRECT matvecProductDataField (7): check matrix inverse property\n\n";
1697  errorFlag = -1000;
1698  }
1699  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d_d, in_c_f_p_d, 't');
1700  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_p_d_d, out_c_f_p_d, 't');
1701  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
1702  if (rst::vectorNorm(outi_c_f_p_d,NORM_ONE) > zero) {
1703  *outStream << "\n\nINCORRECT matvecProductDataField (8): check matrix inverse property, w/ double transpose\n\n";
1704  errorFlag = -1000;
1705  }
1706  // Tensor values do not vary by point: test "N" and "T" options (no transpose/transpose)
1707  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d_d, in_c_f_p_d);
1708  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_1_d_d, out_c_f_p_d);
1709  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
1710  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
1711  *outStream << "\n\nINCORRECT matvecProductDataField (9): check matrix inverse property\n\n";
1712  errorFlag = -1000;
1713  }
1714  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d_d, in_c_f_p_d, 't');
1715  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_1_d_d, out_c_f_p_d, 't');
1716  rst::subtract(outi_c_f_p_d,in_c_f_p_d);
1717  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
1718  *outStream << "\n\nINCORRECT matvecProductDataField (10): check matrix inverse property, w/ double transpose\n\n";
1719  errorFlag = -1000;
1720  }
1721  /***********************************************************************************************
1722  * Full tensor: inputData(C,P,D,D) test inverse transpose *
1723  **********************************************************************************************/
1724 
1725  for (unsigned int i=0; i<in_c_f_p_d.dimension(0); i++)
1726  for (unsigned int j=0; j<in_c_f_p_d.dimension(1); j++)
1727  for (unsigned int k=0; k<in_c_f_p_d.dimension(2); k++)
1728  for (unsigned int l=0; l<in_c_f_p_d.dimension(3); l++)
1729  in_c_f_p_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
1730 
1731 
1732  for (int ic=0; ic < c; ic++) {
1733  for (int ip=0; ip < p; ip++) {
1734  for (int i=0; i<d1; i++) {
1735  for (int j=0; j<d1; j++) {
1736  data_c_p_d_d(ic, ip, i, j) = Teuchos::ScalarTraits<double>::random();
1737  }
1738  }
1739 
1740  }
1741  }
1742  RealSpaceTools<double>::inverse(datainv_c_p_d_d, data_c_p_d_d);
1743  RealSpaceTools<double>::transpose(datainvtrn_c_p_d_d, datainv_c_p_d_d);
1744  for (int ic=0; ic < c; ic++) {
1745  for (int ip=0; ip < 1; ip++) {
1746  for (int i=0; i<d1; i++) {
1747  for (int j=0; j<d1; j++) {
1748  data_c_1_d_d(ic, 0, i, j) = Teuchos::ScalarTraits<double>::random();
1749  }
1750  }
1751  }
1752 }
1753  RealSpaceTools<double>::inverse(datainv_c_1_d_d, data_c_1_d_d);
1754  RealSpaceTools<double>::transpose(datainvtrn_c_1_d_d, datainv_c_1_d_d);
1755 
1756  // Tensor values vary by point:
1757  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d_d, in_c_f_p_d);
1758  art::matvecProductDataField<double>(outi_c_f_p_d, datainvtrn_c_p_d_d, out_c_f_p_d, 't');
1759  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
1760  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
1761  *outStream << "\n\nINCORRECT matvecProductDataField (11): check matrix inverse transpose property\n\n";
1762  errorFlag = -1000;
1763  }
1764  // Tensor values do not vary by point
1765  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d_d, in_c_f_p_d);
1766  art::matvecProductDataField<double>(outi_c_f_p_d, datainvtrn_c_1_d_d, out_c_f_p_d, 't');
1767  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
1768  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
1769  *outStream << "\n\nINCORRECT matvecProductDataField (12): check matrix inverse transpose property\n\n";
1770  errorFlag = -1000;
1771  }
1772  }// test 5.c scope
1773 
1774 
1775  *outStream \
1776  << "\n"
1777  << "===============================================================================\n"\
1778  << "| TEST 5.d: matvecProductDataField random tests: branch inputFields(F,P,D) |\n"\
1779  << "===============================================================================\n";
1780  /*
1781  * d1 is the spatial dimension and should be 1, 2 or 3. If d1>3, the RealSpaceTools function 'inverse' will fail
1782  */
1783  {// test 5.d scope
1784  int c=5, p=9, f=7, d1=3;
1785  double zero = INTREPID_TOL*10000.0;
1786 
1787  Kokkos::View<double***> in_f_p_d("in_f_p_d",f, p, d1);
1788  Kokkos::View<double****> in_c_f_p_d("in_c_f_p_d",c, f, p, d1);
1789  Kokkos::View<double**> data_c_p("data_c_p",c, p);
1790  Kokkos::View<double**> datainv_c_p("datainv_c_p",c, p);
1791  Kokkos::View<double**> data_c_1("data_c_1",c, 1);
1792  Kokkos::View<double**> datainv_c_1("datainv_c_1",c, 1);
1793  Kokkos::View<double***> data_c_p_d("data_c_p_d",c, p, d1);
1794  Kokkos::View<double***> datainv_c_p_d("datainv_c_p_d",c, p, d1);
1795  Kokkos::View<double***> data_c_1_d("data_c_1_d",c, 1, d1);
1796  Kokkos::View<double***> datainv_c_1_d("datainv_c_1_d",c, 1, d1);
1797  Kokkos::View<double****> data_c_p_d_d("data_c_p_d_d",c, p, d1, d1);
1798  Kokkos::View<double****> datainv_c_p_d_d("datainv_c_p_d_d",c, p, d1, d1);
1799  Kokkos::View<double****> datainvtrn_c_p_d_d("datainvtrn_c_p_d_d",c, p, d1, d1);
1800  Kokkos::View<double****> data_c_1_d_d("data_c_1_d_d",c, 1, d1, d1);
1801  Kokkos::View<double****> datainv_c_1_d_d("datainv_c_1_d_d",c, 1, d1, d1);
1802  Kokkos::View<double****> datainvtrn_c_1_d_d("datainvtrn_c_1_d_d",c, 1, d1, d1);
1803  Kokkos::View<double**> data_c_p_one("data_c_p_one",c, p);
1804  Kokkos::View<double**> data_c_1_one("data_c_1_one",c, 1);
1805  Kokkos::View<double****> out_c_f_p_d("out_c_f_p_d",c, f, p, d1);
1806  Kokkos::View<double****> outi_c_f_p_d("outi_c_f_p_d",c, f, p, d1);
1807  /***********************************************************************************************
1808  * Constant diagonal tensor: inputData(C,P) *
1809  **********************************************************************************************/
1810 
1811  for (unsigned int i=0; i<in_f_p_d.dimension(0); i++)
1812  for (unsigned int j=0; j<in_f_p_d.dimension(1); j++)
1813  for (unsigned int k=0; k<in_f_p_d.dimension(2); k++)
1814  in_f_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
1815 
1816  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
1817  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
1818  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
1819  datainv_c_p(i,j) = 1.0 / data_c_p(i,j);
1820  data_c_p_one(i,j) = 1.0;
1821  }
1822 
1823  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
1824  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
1825  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
1826  datainv_c_1(i,j) = 1.0 / data_c_1(i,j);
1827  }
1828 
1829 
1830  // Tensor values vary by point
1831  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1832  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p, in_f_p_d);
1833  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_p, out_c_f_p_d);
1834  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
1835  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
1836  *outStream << "\n\nINCORRECT matvecProductDataField (13): check scalar inverse property\n\n";
1837  errorFlag = -1000;
1838  }
1839  // Tensor values do not vary by point
1840  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1841  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1, in_f_p_d);
1842  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_1, out_c_f_p_d);
1843  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
1844  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
1845  *outStream << "\n\nINCORRECT matvecProductDataField (14): check scalar inverse property\n\n";
1846  errorFlag = -1000;
1847  }
1848  /***********************************************************************************************
1849  * Non-constant diagonal tensor: inputData(C,P,D) *
1850  **********************************************************************************************/
1851  for (unsigned int i=0; i<in_f_p_d.dimension(0); i++)
1852  for (unsigned int j=0; j<in_f_p_d.dimension(1); j++)
1853  for (unsigned int k=0; k<in_f_p_d.dimension(2); k++){
1854  in_f_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
1855  }
1856 
1857  for (unsigned int i=0; i<data_c_p_d.dimension(0); i++)
1858  for (unsigned int j=0; j<data_c_p_d.dimension(1); j++)
1859  for (unsigned int k=0; k<data_c_p_d.dimension(2); k++){
1860  data_c_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
1861  datainv_c_p_d(i,j,k) = 1.0 / data_c_p_d(i,j,k);
1862  }
1863 
1864  for (unsigned int i=0; i<data_c_1_d.dimension(0); i++)
1865  for (unsigned int j=0; j<data_c_1_d.dimension(1); j++)
1866  for (unsigned int k=0; k<data_c_1_d.dimension(2); k++){
1867  data_c_1_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
1868  datainv_c_1_d(i,j,k) = 1.0 / data_c_1_d(i,j,k);
1869  }
1870 
1871  // Tensor values vary by point:
1872  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1873  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d, in_f_p_d);
1874  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_p_d, out_c_f_p_d);
1875  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
1876  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
1877  *outStream << "\n\nINCORRECT matvecProductDataField (15): check scalar inverse property\n\n";
1878  errorFlag = -1000;
1879  }
1880  // Tensor values do not vary by point:
1881  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1882  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d, in_f_p_d);
1883  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_1_d, out_c_f_p_d);
1884  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
1885  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
1886  *outStream << "\n\nINCORRECT matvecProductDataField (16): check scalar inverse property\n\n";
1887  errorFlag = -1000;
1888  }
1889  /***********************************************************************************************
1890  * Full tensor: inputData(C,P,D,D) *
1891  **********************************************************************************************/
1892  for (unsigned int i=0; i<in_f_p_d.dimension(0); i++)
1893  for (unsigned int j=0; j<in_f_p_d.dimension(1); j++)
1894  for (unsigned int k=0; k<in_f_p_d.dimension(2); k++){
1895  in_f_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
1896  }
1897 
1898  for (int ic=0; ic < c; ic++) {
1899  for (int ip=0; ip < p; ip++) {
1900  for (int i=0; i<d1; i++) {
1901  for (int j=0; j<d1; j++) {
1902  data_c_p_d_d(ic, ip, i, j) = Teuchos::ScalarTraits<double>::random();
1903  }
1904  }
1905  }
1906  }
1907  RealSpaceTools<double>::inverse(datainv_c_p_d_d, data_c_p_d_d);
1908  for (int ic=0; ic < c; ic++) {
1909  for (int ip=0; ip < 1; ip++) {
1910  for (int i=0; i<d1; i++) {
1911  for (int j=0; j<d1; j++) {
1912  data_c_1_d_d(ic, 0, i, j) = Teuchos::ScalarTraits<double>::random();
1913  }
1914  }
1915  }
1916  }
1917  RealSpaceTools<double>::inverse(datainv_c_1_d_d, data_c_1_d_d);
1918  // Tensor values vary by point: test "N" and "T" (no-transpose/transpose) options
1919  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1920  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d_d, in_f_p_d);
1921  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_p_d_d, out_c_f_p_d);
1922  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
1923  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
1924  *outStream << "\n\nINCORRECT matvecProductDataField (17): check matrix inverse property\n\n";
1925  errorFlag = -1000;
1926  }
1927  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1928  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d_d, in_f_p_d, 't');
1929  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_p_d_d, out_c_f_p_d, 't');
1930  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
1931  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
1932  *outStream << "\n\nINCORRECT matvecProductDataField (18): check matrix inverse property, w/ double transpose\n\n";
1933  errorFlag = -1000;
1934  }
1935  // Tensor values do not vary by point: test "N" and "T" (no-transpose/transpose) options
1936  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1937  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d_d, in_f_p_d);
1938  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_1_d_d, out_c_f_p_d);
1939  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
1940  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
1941  *outStream << "\n\nINCORRECT matvecProductDataField (19): check matrix inverse property\n\n";
1942  errorFlag = -1000;
1943  }
1944  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1945  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d_d, in_f_p_d, 't');
1946  art::matvecProductDataField<double>(outi_c_f_p_d, datainv_c_1_d_d, out_c_f_p_d, 't');
1947  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
1948  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
1949  *outStream << "\n\nINCORRECT matvecProductDataField (20): check matrix inverse property, w/ double transpose\n\n";
1950  errorFlag = -1000;
1951  }
1952  /***********************************************************************************************
1953  * Full tensor: inputData(C,P,D,D) test inverse transpose *
1954  **********************************************************************************************/
1955  for (unsigned int i=0; i<in_f_p_d.dimension(0); i++)
1956  for (unsigned int j=0; j<in_f_p_d.dimension(1); j++)
1957  for (unsigned int k=0; k<in_f_p_d.dimension(2); k++){
1958  in_f_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
1959  }
1960 
1961  for (int ic=0; ic < c; ic++) {
1962  for (int ip=0; ip < p; ip++) {
1963  for (int i=0; i<d1; i++) {
1964  for (int j=0; j<d1; j++) {
1965  data_c_p_d_d(ic, ip, i, j) = Teuchos::ScalarTraits<double>::random();
1966  }
1967  }
1968  }
1969  }
1970 RealSpaceTools<double>::inverse(datainv_c_p_d_d, data_c_p_d_d);
1971 RealSpaceTools<double>::transpose(datainvtrn_c_p_d_d, datainv_c_p_d_d);
1972  for (int ic=0; ic < c; ic++) {
1973  for (int ip=0; ip < 1; ip++) {
1974  for (int i=0; i<d1; i++) {
1975  for (int j=0; j<d1; j++) {
1976  data_c_1_d_d(ic, 0, i, j) = Teuchos::ScalarTraits<double>::random();
1977  }
1978  }
1979  }
1980  }
1981  RealSpaceTools<double>::inverse(datainv_c_1_d_d, data_c_1_d_d);
1982  RealSpaceTools<double>::transpose(datainvtrn_c_1_d_d, datainv_c_1_d_d);
1983 
1984  // Tensor values vary by point:
1985  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1986  art::matvecProductDataField<double>(out_c_f_p_d, data_c_p_d_d, in_f_p_d);
1987  art::matvecProductDataField<double>(outi_c_f_p_d, datainvtrn_c_p_d_d, out_c_f_p_d, 't');
1988  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
1989  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
1990  *outStream << "\n\nINCORRECT matvecProductDataField (21): check matrix inverse transpose property\n\n";
1991  errorFlag = -1000;
1992  }
1993  // Tensor values do not vary by point:
1994  art::matvecProductDataField<double>(in_c_f_p_d, data_c_p_one, in_f_p_d);
1995  art::matvecProductDataField<double>(out_c_f_p_d, data_c_1_d_d, in_f_p_d);
1996  art::matvecProductDataField<double>(outi_c_f_p_d, datainvtrn_c_1_d_d, out_c_f_p_d, 't');
1997  rst::subtract(outi_c_f_p_d, in_c_f_p_d);
1998  if (rst::vectorNorm(outi_c_f_p_d, NORM_ONE) > zero) {
1999  *outStream << "\n\nINCORRECT matvecProductDataField (22): check matrix inverse transpose property\n\n";
2000  errorFlag = -1000;
2001  }
2002  }// test 5.d scope
2003 
2004 
2005 
2006  *outStream \
2007  << "\n"
2008  << "===============================================================================\n"\
2009  << "| TEST 6.a: matvecProductDataData operations: (C,P,D,D) and (C,P,D) |\n"\
2010  << "===============================================================================\n";
2011  /*
2012  * inputMat inputVecFields outFields
2013  * 1 1 1 0 0 0 1 1 1 0 0 0 0 1 -1 -1 0 0 -3 0
2014  * -1 2 -1 -1 -2 -3 -1 2 -1 -1 -2 -3 0 1 -1 1 0 -6 0 2
2015  * 1 2 3 -2 6 -4 1 2 3 -2 6 -4 0 1 -1 -1 0 0 -6 12
2016  */
2017 
2018  // (C,P,D,D)
2019 
2020  Kokkos::realloc(inputMat,4,1,3,3);
2021  // cell 0
2022  inputMat(0,0,0,0) = 1.0; inputMat(0,0,0,1) = 1.0; inputMat(0,0,0,2) = 1.0;
2023  inputMat(0,0,1,0) =-1.0; inputMat(0,0,1,1) = 2.0; inputMat(0,0,1,2) =-1.0;
2024  inputMat(0,0,2,0) = 1.0; inputMat(0,0,2,1) = 2.0; inputMat(0,0,2,2) = 3.0;
2025  // cell 1
2026  inputMat(1,0,0,0) = 0.0; inputMat(1,0,0,1) = 0.0; inputMat(1,0,0,2) = 0.0;
2027  inputMat(1,0,1,0) =-1.0; inputMat(1,0,1,1) =-2.0; inputMat(1,0,1,2) =-3.0;
2028  inputMat(1,0,2,0) =-2.0; inputMat(1,0,2,1) = 6.0; inputMat(1,0,2,2) =-4.0;
2029  // cell 2
2030  inputMat(2,0,0,0) = 1.0; inputMat(2,0,0,1) = 1.0; inputMat(2,0,0,2) = 1.0;
2031  inputMat(2,0,1,0) =-1.0; inputMat(2,0,1,1) = 2.0; inputMat(2,0,1,2) =-1.0;
2032  inputMat(2,0,2,0) = 1.0; inputMat(2,0,2,1) = 2.0; inputMat(2,0,2,2) = 3.0;
2033  // cell 3
2034  inputMat(3,0,0,0) = 0.0; inputMat(3,0,0,1) = 0.0; inputMat(3,0,0,2) = 0.0;
2035  inputMat(3,0,1,0) =-1.0; inputMat(3,0,1,1) =-2.0; inputMat(3,0,1,2) =-3.0;
2036  inputMat(3,0,2,0) =-2.0; inputMat(3,0,2,1) = 6.0; inputMat(3,0,2,2) =-4.0;
2037 
2038  // (C,P,D)
2039  Kokkos::realloc(inputVecFields3,4,1,3);
2040  inputVecFields3(0,0,0) = 0.0; inputVecFields3(0,0,1) = 0.0; inputVecFields3(0,0,2) = 0.0;
2041  inputVecFields3(1,0,0) = 1.0; inputVecFields3(1,0,1) = 1.0; inputVecFields3(1,0,2) = 1.0;
2042  inputVecFields3(2,0,0) =-1.0; inputVecFields3(2,0,1) =-1.0; inputVecFields3(2,0,2) =-1.0;
2043  inputVecFields3(3,0,0) =-1.0; inputVecFields3(3,0,1) = 1.0; inputVecFields3(3,0,2) =-1.0;
2044 
2045  // (C,P,D) - true
2046  Kokkos::View<double***>outFieldsCorrect3("outFieldsCorrect3",4,1,3);
2047  outFieldsCorrect3(0,0,0) = 0.0; outFieldsCorrect3(0,0,1) = 0.0; outFieldsCorrect3(0,0,2) = 0.0;
2048  outFieldsCorrect3(1,0,0) = 0.0; outFieldsCorrect3(1,0,1) =-6.0; outFieldsCorrect3(1,0,2) = 0.0;
2049  outFieldsCorrect3(2,0,0) =-3.0; outFieldsCorrect3(2,0,1) = 0.0; outFieldsCorrect3(2,0,2) =-6.0;
2050  outFieldsCorrect3(3,0,0) = 0.0; outFieldsCorrect3(3,0,1) = 2.0; outFieldsCorrect3(3,0,2) = 12.0;
2051 
2052  // (C,P,D)
2053  Kokkos::realloc(outFields3,4,1,3);
2054  art::matvecProductDataData<double>(outFields3, inputMat, inputVecFields3);
2055  // test loop
2056  for(unsigned int cell = 0; cell < outFields3.dimension(0); cell++){
2057  for(unsigned int point = 0; point < outFields3.dimension(1); point++){
2058  for(unsigned int row = 0; row < outFields3.dimension(2); row++){
2059  if(outFields3(cell, point, row) != outFieldsCorrect3(cell, point, row)) {
2060  *outStream << "\n\nINCORRECT matvecProductDataData (1): \n value at multi-index ("
2061  << cell << "," << point << "," << row << ") = "
2062  << outFields3(cell, point, row) << " but correct value is "
2063  << outFieldsCorrect3(cell, point, row) <<"\n";
2064  errorFlag++;
2065  }
2066  }//row
2067  }// point
2068  }// cell
2069 
2070 
2071  *outStream \
2072  << "\n"
2073  << "===============================================================================\n"\
2074  << "| TEST 6.b: matvecProductDataData operations: (C,P,D,D) and (P,D) |\n"\
2075  << "===============================================================================\n";
2076  /*
2077  * inputMat inputVecFields outFields
2078  * 1 1 1 0 0 0 1 1 1 0 0 0 0 1 -1 -1 0 0 -3 0
2079  * -1 2 -1 -1 -2 -3 -1 2 -1 -1 -2 -3 0 1 -1 1 0 -6 0 2
2080  * 1 2 3 -2 6 -4 1 2 3 -2 6 -4 0 1 -1 -1 0 0 -6 12
2081  */
2082  // (C,P,D,D)
2083  Kokkos::realloc(inputMat,1,4,3,3);
2084  // point 0
2085  inputMat(0,0,0,0) = 1.0; inputMat(0,0,0,1) = 1.0; inputMat(0,0,0,2) = 1.0;
2086  inputMat(0,0,1,0) =-1.0; inputMat(0,0,1,1) = 2.0; inputMat(0,0,1,2) =-1.0;
2087  inputMat(0,0,2,0) = 1.0; inputMat(0,0,2,1) = 2.0; inputMat(0,0,2,2) = 3.0;
2088  // point 1
2089  inputMat(0,1,0,0) = 0.0; inputMat(0,1,0,1) = 0.0; inputMat(0,1,0,2) = 0.0;
2090  inputMat(0,1,1,0) =-1.0; inputMat(0,1,1,1) =-2.0; inputMat(0,1,1,2) =-3.0;
2091  inputMat(0,1,2,0) =-2.0; inputMat(0,1,2,1) = 6.0; inputMat(0,1,2,2) =-4.0;
2092  // point 2
2093  inputMat(0,2,0,0) = 1.0; inputMat(0,2,0,1) = 1.0; inputMat(0,2,0,2) = 1.0;
2094  inputMat(0,2,1,0) =-1.0; inputMat(0,2,1,1) = 2.0; inputMat(0,2,1,2) =-1.0;
2095  inputMat(0,2,2,0) = 1.0; inputMat(0,2,2,1) = 2.0; inputMat(0,2,2,2) = 3.0;
2096  // point 3
2097  inputMat(0,3,0,0) = 0.0; inputMat(0,3,0,1) = 0.0; inputMat(0,3,0,2) = 0.0;
2098  inputMat(0,3,1,0) =-1.0; inputMat(0,3,1,1) =-2.0; inputMat(0,3,1,2) =-3.0;
2099  inputMat(0,3,2,0) =-2.0; inputMat(0,3,2,1) = 6.0; inputMat(0,3,2,2) =-4.0;
2100 
2101  // (P,D)
2102  Kokkos::View<double**>inputVecFields2("inputVecFields2",4,3);
2103  //
2104  inputVecFields2(0,0) = 0.0; inputVecFields2(0,1) = 0.0; inputVecFields2(0,2) = 0.0;
2105  inputVecFields2(1,0) = 1.0; inputVecFields2(1,1) = 1.0; inputVecFields2(1,2) = 1.0;
2106  inputVecFields2(2,0) =-1.0; inputVecFields2(2,1) =-1.0; inputVecFields2(2,2) =-1.0;
2107  inputVecFields2(3,0) =-1.0; inputVecFields2(3,1) = 1.0; inputVecFields2(3,2) =-1.0;
2108 
2109  // (C,P,D) - true
2110  Kokkos::realloc(outFieldsCorrect3,1,4,3);
2111  outFieldsCorrect3(0,0,0) = 0.0; outFieldsCorrect3(0,0,1) = 0.0; outFieldsCorrect3(0,0,2) = 0.0;
2112  outFieldsCorrect3(0,1,0) = 0.0; outFieldsCorrect3(0,1,1) =-6.0; outFieldsCorrect3(0,1,2) = 0.0;
2113  outFieldsCorrect3(0,2,0) =-3.0; outFieldsCorrect3(0,2,1) = 0.0; outFieldsCorrect3(0,2,2) =-6.0;
2114  outFieldsCorrect3(0,3,0) = 0.0; outFieldsCorrect3(0,3,1) = 2.0; outFieldsCorrect3(0,3,2) = 12.0;
2115 
2116  // (C,P,D)
2117  Kokkos::realloc(outFields3,1,4,3);
2118  art::matvecProductDataData<double>(outFields3, inputMat, inputVecFields2);
2119  // test loop
2120  for(unsigned int cell = 0; cell < outFields3.dimension(0); cell++){
2121  for(unsigned int point = 0; point < outFields3.dimension(1); point++){
2122  for(unsigned int row = 0; row < outFields3.dimension(2); row++){
2123  if(outFields3(cell, point, row) != outFieldsCorrect3(cell, point, row)) {
2124  *outStream << "\n\nINCORRECT matvecProductDataData (2): \n value at multi-index ("
2125  << cell << "," << point << "," << row << ") = "
2126  << outFields3(cell, point, row) << " but correct value is "
2127  << outFieldsCorrect3(cell, point, row) <<"\n";
2128  errorFlag++;
2129  }
2130  }//row
2131  }// point
2132  }// cell
2133 
2134 
2135  *outStream \
2136  << "\n"
2137  << "===============================================================================\n"\
2138  << "| TEST 6.c: matvecProductDataData random tests: branch inputDataRight(C,P,D) |\n"\
2139  << "===============================================================================\n";
2140  /*
2141  * Test derived from Test 5.c
2142  */
2143  {// test 6.c scope
2144  int c=5, p=9, d1=3;
2145  double zero = INTREPID_TOL*10000.0;
2146 
2147  Kokkos::View<double***> in_c_p_d("in_c_p_d",c, p, d1);
2148  Kokkos::View<double***> out_c_p_d("out_c_p_d",c, p, d1);
2149  Kokkos::View<double***> outi_c_p_d("outi_c_p_d",c, p, d1);
2150 
2151  Kokkos::View<double**> data_c_p("data_c_p",c, p);
2152  Kokkos::View<double**> datainv_c_p("datainv_c_p",c, p);
2153  Kokkos::View<double**> data_c_1("data_c_1",c, 1);
2154  Kokkos::View<double**> datainv_c_1("datainv_c_1",c, 1);
2155  Kokkos::View<double***> data_c_p_d("data_c_p_d",c, p, d1);
2156  Kokkos::View<double***> datainv_c_p_d("datainv_c_p_d",c, p, d1);
2157  Kokkos::View<double***> data_c_1_d("data_c_1_d",c, 1, d1);
2158  Kokkos::View<double***> datainv_c_1_d("datainv_c_1_d",c, 1, d1);
2159  Kokkos::View<double****> data_c_p_d_d("data_c_p_d_d",c, p, d1, d1);
2160  Kokkos::View<double****> datainv_c_p_d_d("datainv_c_p_d_d",c, p, d1, d1);
2161  Kokkos::View<double****> datainvtrn_c_p_d_d("datainvtrn_c_p_d_d",c, p, d1, d1);
2162  Kokkos::View<double****> data_c_1_d_d("data_c_1_d_d",c, 1, d1, d1);
2163  Kokkos::View<double****> datainv_c_1_d_d("datainv_c_1_d_d",c, 1, d1, d1);
2164  Kokkos::View<double****> datainvtrn_c_1_d_d("datainvtrn_c_1_d_d",c, 1, d1, d1);
2165  /***********************************************************************************************
2166  * Constant diagonal tensor: inputDataLeft(C,P) *
2167  **********************************************************************************************/
2168  for (unsigned int i=0; i<in_c_p_d.dimension(0); i++)
2169  for (unsigned int j=0; j<in_c_p_d.dimension(1); j++)
2170  for (unsigned int k=0; k<in_c_p_d.dimension(2); k++){
2171  in_c_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
2172  }
2173 
2174  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
2175  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
2176  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
2177  datainv_c_p(i,j) = 1.0 / data_c_p(i,j);
2178  }
2179 
2180  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
2181  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
2182  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
2183  datainv_c_1(i,j) = 1.0 / data_c_1(i,j);
2184  }
2185  // Tensor values vary by point:
2186  art::matvecProductDataData<double>(out_c_p_d, data_c_p, in_c_p_d);
2187  art::matvecProductDataData<double>(out_c_p_d, datainv_c_p, out_c_p_d);
2188  rst::subtract(out_c_p_d, in_c_p_d);
2189  if (rst::vectorNorm(out_c_p_d, NORM_ONE) > zero) {
2190  *outStream << "\n\nINCORRECT matvecProductDataData (3): check scalar inverse property\n\n";
2191  errorFlag = -1000;
2192  }
2193  // Tensor values do not vary by point:
2194  art::matvecProductDataData<double>(out_c_p_d, data_c_1, in_c_p_d);
2195  art::matvecProductDataData<double>(out_c_p_d, datainv_c_1, out_c_p_d);
2196  rst::subtract(out_c_p_d, in_c_p_d);
2197  if (rst::vectorNorm(out_c_p_d, NORM_ONE) > zero) {
2198  *outStream << "\n\nINCORRECT matvecProductDataData (4): check scalar inverse property\n\n";
2199  errorFlag = -1000;
2200  }
2201  /***********************************************************************************************
2202  * Non-onstant diagonal tensor: inputDataLeft(C,P,D) *
2203  **********************************************************************************************/
2204  for (unsigned int i=0; i<in_c_p_d.dimension(0); i++)
2205  for (unsigned int j=0; j<in_c_p_d.dimension(1); j++)
2206  for (unsigned int k=0; k<in_c_p_d.dimension(2); k++){
2207  in_c_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
2208  }
2209 
2210  for (unsigned int i=0; i<data_c_p_d.dimension(0); i++)
2211  for (unsigned int j=0; j<data_c_p_d.dimension(1); j++)
2212  for (unsigned int k=0; k<data_c_p_d.dimension(2); k++){
2213  data_c_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
2214  datainv_c_p_d(i,j,k) = 1.0 / data_c_p_d(i,j,k);
2215  }
2216 
2217  for (unsigned int i=0; i<data_c_1_d.dimension(0); i++)
2218  for (unsigned int j=0; j<data_c_1_d.dimension(1); j++)
2219  for (unsigned int k=0; k<data_c_1_d.dimension(2); k++){
2220  data_c_1_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
2221  datainv_c_1_d(i,j,k) = 1.0 / data_c_1_d(i,j,k);
2222  }
2223  // Tensor values vary by point:
2224  art::matvecProductDataData<double>(out_c_p_d, data_c_p_d, in_c_p_d);
2225  art::matvecProductDataData<double>(out_c_p_d, datainv_c_p_d, out_c_p_d);
2226  rst::subtract(out_c_p_d, in_c_p_d);
2227  if (rst::vectorNorm(out_c_p_d, NORM_ONE) > zero) {
2228  *outStream << "\n\nINCORRECT matvecProductDataData (5): check scalar inverse property\n\n";
2229  errorFlag = -1000;
2230  }
2231  // Tensor values do not vary by point
2232  art::matvecProductDataData<double>(out_c_p_d, data_c_1_d, in_c_p_d);
2233  art::matvecProductDataData<double>(out_c_p_d, datainv_c_1_d, out_c_p_d);
2234  rst::subtract(out_c_p_d, in_c_p_d);
2235  if (rst::vectorNorm(out_c_p_d, NORM_ONE) > zero) {
2236  *outStream << "\n\nINCORRECT matvecProductDataData (6): check scalar inverse property\n\n";
2237  errorFlag = -1000;
2238  }
2239  /***********************************************************************************************
2240  * Full tensor: inputDataLeft(C,P,D,D) *
2241  **********************************************************************************************/
2242  for (unsigned int i=0; i<in_c_p_d.dimension(0); i++)
2243  for (unsigned int j=0; j<in_c_p_d.dimension(1); j++)
2244  for (unsigned int k=0; k<in_c_p_d.dimension(2); k++){
2245  in_c_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
2246  }
2247  for (int ic=0; ic < c; ic++) {
2248  for (int ip=0; ip < p; ip++) {
2249  for (int i=0; i<d1; i++) {
2250  for (int j=0; j<d1; j++) {
2251  data_c_p_d_d(ic, ip, i, j) = Teuchos::ScalarTraits<double>::random();
2252  }
2253  }
2254  }
2255  }
2256 
2257  RealSpaceTools<double>::inverse(datainv_c_p_d_d, data_c_p_d_d);
2258 
2259  for (int ic=0; ic < c; ic++) {
2260  for (int ip=0; ip < 1; ip++) {
2261  for (int i=0; i<d1; i++) {
2262  for (int j=0; j<d1; j++) {
2263  data_c_1_d_d(ic, 0, i, j) = Teuchos::ScalarTraits<double>::random();
2264  }
2265  }
2266  }
2267  }
2268 
2269  RealSpaceTools<double>::inverse(datainv_c_1_d_d, data_c_1_d_d);
2270  // Tensor values vary by point: test "N" and "T" options (no transpose/transpose)
2271  art::matvecProductDataData<double>(out_c_p_d, data_c_p_d_d, in_c_p_d);
2272  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_p_d_d, out_c_p_d);
2273  rst::subtract(outi_c_p_d, in_c_p_d);
2274  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
2275  *outStream << "\n\nINCORRECT matvecProductDataData (7): check matrix inverse property\n\n";
2276  errorFlag = -1000;
2277  }
2278  art::matvecProductDataData<double>(out_c_p_d, data_c_p_d_d, in_c_p_d, 't');
2279  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_p_d_d, out_c_p_d, 't');
2280  rst::subtract(outi_c_p_d, in_c_p_d);
2281  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
2282  *outStream << "\n\nINCORRECT matvecProductDataData (8): check matrix inverse property, w/ double transpose\n\n";
2283  errorFlag = -1000;
2284  }
2285  // Tensor values do not vary by point: test "N" and "T" options (no transpose/transpose)
2286  art::matvecProductDataData<double>(out_c_p_d, data_c_1_d_d, in_c_p_d);
2287  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_1_d_d, out_c_p_d);
2288  rst::subtract(outi_c_p_d, in_c_p_d);
2289  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
2290  *outStream << "\n\nINCORRECT matvecProductDataData (9): check matrix inverse property\n\n";
2291  errorFlag = -1000;
2292  }
2293  art::matvecProductDataData<double>(out_c_p_d, data_c_1_d_d, in_c_p_d, 't');
2294  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_1_d_d, out_c_p_d, 't');
2295  rst::subtract(outi_c_p_d, in_c_p_d);
2296  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
2297  *outStream << "\n\nINCORRECT matvecProductDataData (10): check matrix inverse property, w/ double transpose\n\n";
2298  errorFlag = -1000;
2299  }
2300 
2301  /***********************************************************************************************
2302  * Full tensor: inputDataLeft(C,P,D,D) test inverse transpose *
2303  **********************************************************************************************/
2304  for (unsigned int i=0; i<in_c_p_d.dimension(0); i++)
2305  for (unsigned int j=0; j<in_c_p_d.dimension(1); j++)
2306  for (unsigned int k=0; k<in_c_p_d.dimension(2); k++){
2307  in_c_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
2308  }
2309 
2310  for (int ic=0; ic < c; ic++) {
2311  for (int ip=0; ip < p; ip++) {
2312  for (int i=0; i<d1; i++) {
2313  for (int j=0; j<d1; j++) {
2314  data_c_p_d_d(ic, ip, i, j) = Teuchos::ScalarTraits<double>::random();
2315  }
2316  }
2317  }
2318  }
2319  RealSpaceTools<double>::inverse(datainv_c_p_d_d, data_c_p_d_d);
2320  RealSpaceTools<double>::transpose(datainvtrn_c_p_d_d, datainv_c_p_d_d);
2321 
2322  for (int ic=0; ic < c; ic++) {
2323  for (int ip=0; ip < 1; ip++) {
2324  for (int i=0; i<d1; i++) {
2325  for (int j=0; j<d1; j++) {
2326  data_c_1_d_d(ic, 0, i, j) = Teuchos::ScalarTraits<double>::random();
2327  }
2328  }
2329  }
2330  }
2331  RealSpaceTools<double>::inverse(datainv_c_1_d_d, data_c_1_d_d);
2332  RealSpaceTools<double>::transpose(datainvtrn_c_1_d_d, datainv_c_1_d_d);
2333  // Tensor values vary by point:
2334  art::matvecProductDataData<double>(out_c_p_d, data_c_p_d_d, in_c_p_d);
2335  art::matvecProductDataData<double>(outi_c_p_d, datainvtrn_c_p_d_d, out_c_p_d, 't');
2336  rst::subtract(outi_c_p_d, in_c_p_d);
2337  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
2338  *outStream << "\n\nINCORRECT matvecProductDataData (11): check matrix inverse transpose property\n\n";
2339  errorFlag = -1000;
2340  }
2341  // Tensor values do not vary by point
2342  art::matvecProductDataData<double>(out_c_p_d, data_c_1_d_d, in_c_p_d);
2343  art::matvecProductDataData<double>(outi_c_p_d, datainvtrn_c_1_d_d, out_c_p_d, 't');
2344  rst::subtract(outi_c_p_d, in_c_p_d);
2345  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
2346  *outStream << "\n\nINCORRECT matvecProductDataData (12): check matrix inverse transpose property\n\n";
2347  errorFlag = -1000;
2348  }
2349  }// test 6.c scope
2350 
2351 
2352  *outStream \
2353  << "\n"
2354  << "===============================================================================\n"\
2355  << "| TEST 6.d: matvecProductDataData random tests: branch inputDataRight(P,D) |\n"\
2356  << "===============================================================================\n";
2357  /*
2358  * Test derived from Test 5.d
2359  */
2360  {// test 6.d scope
2361  int c=5, p=9, d1=3;
2362  double zero = INTREPID_TOL*10000.0;
2363 
2364  Kokkos::View<double**> in_p_d("in_p_d",p, d1);
2365  Kokkos::View<double***> in_c_p_d("in_c_p_d",c, p, d1);
2366  Kokkos::View<double**> data_c_p("data_c_p",c, p);
2367  Kokkos::View<double**> datainv_c_p("datainv_c_p",c, p);
2368  Kokkos::View<double**> data_c_1("data_c_1",c, 1);
2369  Kokkos::View<double**> datainv_c_1("datainv_c_1",c, 1);
2370  Kokkos::View<double***> data_c_p_d("data_c_p_d",c, p, d1);
2371  Kokkos::View<double***> datainv_c_p_d("datainv_c_p_d",c, p, d1);
2372  Kokkos::View<double***> data_c_1_d("data_c_1_d",c, 1, d1);
2373  Kokkos::View<double***> datainv_c_1_d("datainv_c_1_d",c, 1, d1);
2374  Kokkos::View<double****> data_c_p_d_d("data_c_p_d_d",c, p, d1, d1);
2375  Kokkos::View<double****> datainv_c_p_d_d("datainv_c_p_d_d",c, p, d1, d1);
2376  Kokkos::View<double****> datainvtrn_c_p_d_d("datainvtrn_c_p_d_d",c, p, d1, d1);
2377  Kokkos::View<double****> data_c_1_d_d("data_c_1_d_d",c, 1, d1, d1);
2378  Kokkos::View<double****> datainv_c_1_d_d("datainv_c_1_d_d",c, 1, d1, d1);
2379  Kokkos::View<double****> datainvtrn_c_1_d_d("datainvtrn_c_1_d_d",c, 1, d1, d1);
2380  Kokkos::View<double**> data_c_p_one("data_c_p_one",c, p);
2381  Kokkos::View<double**> data_c_1_one("data_c_1_one",c, 1);
2382  Kokkos::View<double***> out_c_p_d("out_c_p_d",c, p, d1);
2383  Kokkos::View<double***> outi_c_p_d("outi_c_p_d",c, p, d1);
2384  /***********************************************************************************************
2385  * Constant diagonal tensor: inputData(C,P) *
2386  **********************************************************************************************/
2387  for (unsigned int i=0; i<in_p_d.dimension(0); i++)
2388  for (unsigned int j=0; j<in_p_d.dimension(1); j++){
2389  in_p_d(i,j) = Teuchos::ScalarTraits<double>::random();
2390  }
2391 
2392  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
2393  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
2394  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
2395  datainv_c_p(i,j) = 1.0 / data_c_p(i,j);
2396  data_c_p_one(i,j) = 1.0;
2397  }
2398 
2399  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
2400  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
2401  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
2402  datainv_c_1(i,j) = 1.0 / data_c_1(i,j);
2403  }
2404  // Tensor values vary by point
2405  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2406  art::matvecProductDataData<double>(out_c_p_d, data_c_p, in_p_d);
2407  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_p, out_c_p_d);
2408  rst::subtract(outi_c_p_d, in_c_p_d);
2409  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
2410  *outStream << "\n\nINCORRECT matvecProductDataData (13): check scalar inverse property\n\n";
2411  errorFlag = -1000;
2412  }
2413  // Tensor values do not vary by point
2414  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2415  art::matvecProductDataData<double>(out_c_p_d, data_c_1, in_p_d);
2416  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_1, out_c_p_d);
2417  rst::subtract(outi_c_p_d, in_c_p_d);
2418  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
2419  *outStream << "\n\nINCORRECT matvecProductDataData (14): check scalar inverse property\n\n";
2420  errorFlag = -1000;
2421  }
2422  /***********************************************************************************************
2423  * Non-constant diagonal tensor: inputData(C,P,D) *
2424  **********************************************************************************************/
2425  for (unsigned int i=0; i<in_p_d.dimension(0); i++)
2426  for (unsigned int j=0; j<in_p_d.dimension(1); j++){
2427  in_p_d(i,j) = Teuchos::ScalarTraits<double>::random();
2428  }
2429 
2430  for (unsigned int i=0; i<data_c_p_d.dimension(0); i++)
2431  for (unsigned int j=0; j<data_c_p_d.dimension(1); j++)
2432  for (unsigned int k=0; k<data_c_p_d.dimension(2); k++){
2433  data_c_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
2434  datainv_c_p_d(i,j,k) = 1.0 / data_c_p_d(i,j,k);
2435  }
2436 
2437  for (unsigned int i=0; i<data_c_1_d.dimension(0); i++)
2438  for (unsigned int j=0; j<data_c_1_d.dimension(1); j++)
2439  for (unsigned int k=0; k<data_c_1_d.dimension(2); k++){
2440  data_c_1_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
2441  datainv_c_1_d(i,j,k) = 1.0 / data_c_1_d(i,j,k);
2442  }
2443  // Tensor values vary by point:
2444  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2445  art::matvecProductDataData<double>(out_c_p_d, data_c_p_d, in_p_d);
2446  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_p_d, out_c_p_d);
2447  rst::subtract(outi_c_p_d, in_c_p_d);
2448  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
2449  *outStream << "\n\nINCORRECT matvecProductDataData (15): check scalar inverse property\n\n";
2450  errorFlag = -1000;
2451  }
2452  // Tensor values do not vary by point:
2453  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2454  art::matvecProductDataData<double>(out_c_p_d, data_c_1_d, in_p_d);
2455  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_1_d, out_c_p_d);
2456  rst::subtract(outi_c_p_d, in_c_p_d);
2457  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
2458  *outStream << "\n\nINCORRECT matvecProductDataData (16): check scalar inverse property\n\n";
2459  errorFlag = -1000;
2460  }
2461  /***********************************************************************************************
2462  * Full tensor: inputData(C,P,D,D) *
2463  **********************************************************************************************/
2464  for (unsigned int i=0; i<in_p_d.dimension(0); i++)
2465  for (unsigned int j=0; j<in_p_d.dimension(1); j++){
2466  in_p_d(i,j) = Teuchos::ScalarTraits<double>::random();
2467  }
2468 
2469  for (int ic=0; ic < c; ic++) {
2470  for (int ip=0; ip < p; ip++) {
2471  for (int i=0; i<d1; i++) {
2472  for (int j=0; j<d1; j++) {
2473  data_c_p_d_d(ic, ip, i, j) = Teuchos::ScalarTraits<double>::random();
2474  }
2475  }
2476  }
2477  }
2478  RealSpaceTools<double>::inverse(datainv_c_p_d_d, data_c_p_d_d);
2479  for (int ic=0; ic < c; ic++) {
2480  for (int ip=0; ip < 1; ip++) {
2481  for (int i=0; i<d1; i++) {
2482  for (int j=0; j<d1; j++) {
2483  data_c_1_d_d(ic, 0, i, j) = Teuchos::ScalarTraits<double>::random();
2484  }
2485  }
2486  }
2487  }
2488 RealSpaceTools<double>::inverse(datainv_c_1_d_d, data_c_1_d_d);
2489 
2490  // Tensor values vary by point: test "N" and "T" (no-transpose/transpose) options
2491  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2492  art::matvecProductDataData<double>(out_c_p_d, data_c_p_d_d, in_p_d);
2493  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_p_d_d, out_c_p_d);
2494  rst::subtract(outi_c_p_d, in_c_p_d);
2495  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
2496  *outStream << "\n\nINCORRECT matvecProductDataData (17): check matrix inverse property\n\n";
2497  errorFlag = -1000;
2498  }
2499  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2500  art::matvecProductDataData<double>(out_c_p_d, data_c_p_d_d, in_p_d, 't');
2501  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_p_d_d, out_c_p_d, 't');
2502  rst::subtract(outi_c_p_d, in_c_p_d);
2503  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
2504  *outStream << "\n\nINCORRECT matvecProductDataData (18): check matrix inverse property, w/ double transpose\n\n";
2505  errorFlag = -1000;
2506  }
2507  // Tensor values do not vary by point: test "N" and "T" (no-transpose/transpose) options
2508  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2509  art::matvecProductDataData<double>(out_c_p_d, data_c_1_d_d, in_p_d);
2510  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_1_d_d, out_c_p_d);
2511  rst::subtract(outi_c_p_d, in_c_p_d);
2512  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
2513  *outStream << "\n\nINCORRECT matvecProductDataData (19): check matrix inverse property\n\n";
2514  errorFlag = -1000;
2515  }
2516  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2517  art::matvecProductDataData<double>(out_c_p_d, data_c_1_d_d, in_p_d, 't');
2518  art::matvecProductDataData<double>(outi_c_p_d, datainv_c_1_d_d, out_c_p_d, 't');
2519  rst::subtract(outi_c_p_d, in_c_p_d);
2520  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
2521  *outStream << "\n\nINCORRECT matvecProductDataData (20): check matrix inverse property, w/ double transpose\n\n";
2522  errorFlag = -1000;
2523  }
2524  /***********************************************************************************************
2525  * Full tensor: inputData(C,P,D,D) test inverse transpose *
2526  **********************************************************************************************/
2527  for (unsigned int i=0; i<in_p_d.dimension(0); i++)
2528  for (unsigned int j=0; j<in_p_d.dimension(1); j++){
2529  in_p_d(i,j) = Teuchos::ScalarTraits<double>::random();
2530  }
2531 
2532  for (int ic=0; ic < c; ic++) {
2533  for (int ip=0; ip < p; ip++) {
2534  for (int i=0; i<d1; i++) {
2535  for (int j=0; j<d1; j++) {
2536  data_c_p_d_d(ic, ip, i, j) = Teuchos::ScalarTraits<double>::random();
2537  }
2538  }
2539  }
2540  }
2541  RealSpaceTools<double>::inverse(datainv_c_p_d_d, data_c_p_d_d);
2542  RealSpaceTools<double>::transpose(datainvtrn_c_p_d_d, datainv_c_p_d_d);
2543 
2544  for (int ic=0; ic < c; ic++) {
2545  for (int ip=0; ip < 1; ip++) {
2546  for (int i=0; i<d1; i++) {
2547  for (int j=0; j<d1; j++) {
2548  data_c_1_d_d(ic, 0, i, j) = Teuchos::ScalarTraits<double>::random();
2549  }
2550  }
2551  }
2552  }
2553  RealSpaceTools<double>::inverse(datainv_c_1_d_d, data_c_1_d_d);
2554  RealSpaceTools<double>::transpose(datainvtrn_c_1_d_d, datainv_c_1_d_d);
2555 
2556  // Tensor values vary by point:
2557  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2558  art::matvecProductDataData<double>(out_c_p_d, data_c_p_d_d, in_p_d);
2559  art::matvecProductDataData<double>(outi_c_p_d, datainvtrn_c_p_d_d, out_c_p_d, 't');
2560  rst::subtract(outi_c_p_d, in_c_p_d);
2561  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
2562  *outStream << "\n\nINCORRECT matvecProductDataData (21): check matrix inverse transpose property\n\n";
2563  errorFlag = -1000;
2564  }
2565  // Tensor values do not vary by point:
2566  art::matvecProductDataData<double>(in_c_p_d, data_c_p_one, in_p_d);
2567  art::matvecProductDataData<double>(out_c_p_d, data_c_1_d_d, in_p_d);
2568  art::matvecProductDataData<double>(outi_c_p_d, datainvtrn_c_1_d_d, out_c_p_d, 't');
2569  rst::subtract(outi_c_p_d, in_c_p_d);
2570  if (rst::vectorNorm(outi_c_p_d, NORM_ONE) > zero) {
2571  *outStream << "\n\nINCORRECT matvecProductDataData (22): check matrix inverse transpose property\n\n";
2572  errorFlag = -1000;
2573  }
2574  }// test 6.d scope
2575 
2576 
2577  *outStream \
2578  << "\n"
2579  << "===============================================================================\n"\
2580  << "| TEST 7.a: matmatProductDataField random tests: branch inputFields(C,F,P,D,D)|\n"\
2581  << "===============================================================================\n";
2582  {// Test 7.a scope
2583  int c=5, p=9, f=7, d1=3;
2584  double zero = INTREPID_TOL*10000.0;
2585 
2586  Kokkos::View<double*****> in_c_f_p_d_d("in_c_f_p_d_d",c, f, p, d1, d1);
2587  Kokkos::View<double**> data_c_p("data_c_p",c, p);
2588  Kokkos::View<double**> datainv_c_p("datainv_c_p",c, p);
2589  Kokkos::View<double**> data_c_1("data_c_1",c, 1);
2590  Kokkos::View<double**> datainv_c_1("datainv_c_1",c, 1);
2591  Kokkos::View<double***> data_c_p_d("data_c_p_d",c, p, d1);
2592  Kokkos::View<double***> datainv_c_p_d("datainv_c_p_d",c, p, d1);
2593  Kokkos::View<double***> data_c_1_d("data_c_1_d",c, 1, d1);
2594  Kokkos::View<double***> datainv_c_1_d("datainv_c_1_d",c, 1, d1);
2595  Kokkos::View<double****> data_c_p_d_d("data_c_p_d_d",c, p, d1, d1);
2596  Kokkos::View<double****> datainv_c_p_d_d("datainv_c_p_d_d",c, p, d1, d1);
2597  Kokkos::View<double****> datainvtrn_c_p_d_d("datainvtrn_c_p_d_d",c, p, d1, d1);
2598  Kokkos::View<double****> data_c_1_d_d("data_c_1_d_d",c, 1, d1, d1);
2599  Kokkos::View<double****> datainv_c_1_d_d("datainv_c_1_d_d",c, 1, d1, d1);
2600  Kokkos::View<double****> datainvtrn_c_1_d_d("datainvtrn_c_1_d_d",c, 1, d1, d1);
2601  Kokkos::View<double*****> out_c_f_p_d_d("out_c_f_p_d_d",c, f, p, d1, d1);
2602  Kokkos::View<double*****> outi_c_f_p_d_d("outi_c_f_p_d_d",c, f, p, d1, d1);
2603  /***********************************************************************************************
2604  * Constant diagonal tensor: inputData(C,P) *
2605  **********************************************************************************************/
2606  for (unsigned int i=0; i<in_c_f_p_d_d.dimension(0); i++)
2607  for (unsigned int j=0; j<in_c_f_p_d_d.dimension(1); j++)
2608  for (unsigned int k=0; k<in_c_f_p_d_d.dimension(2); k++)
2609  for (unsigned int l=0; l<in_c_f_p_d_d.dimension(3); l++)
2610  for (unsigned int m=0; m<in_c_f_p_d_d.dimension(4); m++){
2611  in_c_f_p_d_d(i,j,k,l,m) = Teuchos::ScalarTraits<double>::random();
2612  }
2613 
2614  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
2615  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
2616  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
2617  datainv_c_p(i,j) = 1.0 / data_c_p(i,j);
2618  }
2619 
2620  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
2621  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
2622  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
2623  datainv_c_1(i,j) = 1.0 / data_c_1(i,j);
2624  }
2625 
2626  // Tensor values vary by point:
2627  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p, in_c_f_p_d_d);
2628  art::matmatProductDataField<double>(out_c_f_p_d_d, datainv_c_p, out_c_f_p_d_d);
2629  rst::subtract(out_c_f_p_d_d, in_c_f_p_d_d);
2630  if (rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) > zero) {
2631  *outStream << "\n\nINCORRECT matmatProductDataField (1): check scalar inverse property\n\n";
2632  errorFlag = -1000;
2633  }
2634  // Tensor values do not vary by point:
2635  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1, in_c_f_p_d_d);
2636  art::matmatProductDataField<double>(out_c_f_p_d_d, datainv_c_1, out_c_f_p_d_d);
2637  rst::subtract(out_c_f_p_d_d, in_c_f_p_d_d);
2638  if (rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) > zero) {
2639  *outStream << "\n\nINCORRECT matmatProductDataField (2): check scalar inverse property\n\n";
2640  errorFlag = -1000;
2641  }
2642  /***********************************************************************************************
2643  * Non-onstant diagonal tensor: inputData(C,P,D) *
2644  **********************************************************************************************/
2645  for (unsigned int i=0; i<in_c_f_p_d_d.dimension(0); i++)
2646  for (unsigned int j=0; j<in_c_f_p_d_d.dimension(1); j++)
2647  for (unsigned int k=0; k<in_c_f_p_d_d.dimension(2); k++)
2648  for (unsigned int l=0; l<in_c_f_p_d_d.dimension(3); l++)
2649  for (unsigned int m=0; m<in_c_f_p_d_d.dimension(4); m++){
2650  in_c_f_p_d_d(i,j,k,l,m) = Teuchos::ScalarTraits<double>::random();
2651  }
2652 
2653  for (unsigned int i=0; i<data_c_p_d.dimension(0); i++)
2654  for (unsigned int j=0; j<data_c_p_d.dimension(1); j++)
2655  for (unsigned int k=0; k<data_c_p_d.dimension(2); k++){
2656  data_c_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
2657  datainv_c_p_d(i,j,k) = 1.0 / data_c_p_d(i,j,k);
2658  }
2659 
2660  for (unsigned int i=0; i<data_c_1_d.dimension(0); i++)
2661  for (unsigned int j=0; j<data_c_1_d.dimension(1); j++)
2662  for (unsigned int k=0; k<data_c_1_d.dimension(2); k++){
2663  data_c_1_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
2664  datainv_c_1_d(i,j,k) = 1.0 / data_c_1_d(i,j,k);
2665  }
2666 
2667  // Tensor values vary by point:
2668  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d, in_c_f_p_d_d);
2669  art::matmatProductDataField<double>(out_c_f_p_d_d, datainv_c_p_d, out_c_f_p_d_d);
2670  rst::subtract(out_c_f_p_d_d, in_c_f_p_d_d);
2671  if (rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) > zero) {
2672  *outStream << "\n\nINCORRECT matmatProductDataField (3): check scalar inverse property\n\n";
2673  errorFlag = -1000;
2674  }
2675  // Tensor values do not vary by point
2676  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d, in_c_f_p_d_d);
2677  art::matmatProductDataField<double>(out_c_f_p_d_d, datainv_c_1_d, out_c_f_p_d_d);
2678  rst::subtract(out_c_f_p_d_d, in_c_f_p_d_d);
2679  if (rst::vectorNorm(out_c_f_p_d_d, NORM_ONE) > zero) {
2680  *outStream << "\n\nINCORRECT matmatProductDataField (4): check scalar inverse property\n\n";
2681  errorFlag = -1000;
2682  }
2683  /***********************************************************************************************
2684  * Full tensor: inputData(C,P,D,D) *
2685  **********************************************************************************************/
2686  for (unsigned int i=0; i<in_c_f_p_d_d.dimension(0); i++)
2687  for (unsigned int j=0; j<in_c_f_p_d_d.dimension(1); j++)
2688  for (unsigned int k=0; k<in_c_f_p_d_d.dimension(2); k++)
2689  for (unsigned int l=0; l<in_c_f_p_d_d.dimension(3); l++)
2690  for (unsigned int m=0; m<in_c_f_p_d_d.dimension(4); m++){
2691  in_c_f_p_d_d(i,j,k,l,m) = Teuchos::ScalarTraits<double>::random();
2692  }
2693 
2694  for (int ic=0; ic < c; ic++) {
2695  for (int ip=0; ip < p; ip++) {
2696  for (int i=0; i<d1; i++) {
2697  for (int j=0; j<d1; j++) {
2698  data_c_p_d_d(ic, ip, i, j) = Teuchos::ScalarTraits<double>::random();
2699  }
2700  }
2701  }
2702  }
2703  RealSpaceTools<double>::inverse(datainv_c_p_d_d, data_c_p_d_d);
2704 
2705  for (int ic=0; ic < c; ic++) {
2706  for (int ip=0; ip < 1; ip++) {
2707  for (int i=0; i<d1; i++) {
2708  for (int j=0; j<d1; j++) {
2709  data_c_1_d_d(ic, 0, i, j) = Teuchos::ScalarTraits<double>::random();
2710  }
2711  }
2712  }
2713  }
2714  RealSpaceTools<double>::inverse(datainv_c_1_d_d, data_c_1_d_d);
2715 
2716 
2717  // Tensor values vary by point: test "N" and "T" options (no transpose/transpose)
2718  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d_d, in_c_f_p_d_d);
2719  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_p_d_d, out_c_f_p_d_d);
2720  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
2721  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
2722  *outStream << "\n\nINCORRECT matmatProductDataField (5): check matrix inverse property\n\n";
2723  errorFlag = -1000;
2724  }
2725  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d_d, in_c_f_p_d_d, 't');
2726  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_p_d_d, out_c_f_p_d_d, 't');
2727  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
2728  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
2729  *outStream << "\n\nINCORRECT matmatProductDataField (6): check matrix inverse property, w/ double transpose\n\n";
2730  errorFlag = -1000;
2731  }
2732  // Tensor values do not vary by point: test "N" and "T" options (no transpose/transpose)
2733  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d_d, in_c_f_p_d_d);
2734  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_1_d_d, out_c_f_p_d_d);
2735  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
2736  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
2737  *outStream << "\n\nINCORRECT matmatProductDataField (7): check matrix inverse property\n\n";
2738  errorFlag = -1000;
2739  }
2740  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d_d, in_c_f_p_d_d,'t');
2741  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_1_d_d, out_c_f_p_d_d, 't');
2742  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
2743  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
2744  *outStream << "\n\nINCORRECT matmatProductDataField (8): check matrix inverse property, w/ double transpose\n\n";
2745  errorFlag = -1000;
2746  }
2747  /***********************************************************************************************
2748  * Full tensor: inputData(C,P,D,D) test inverse transpose *
2749  **********************************************************************************************/
2750  for (unsigned int i=0; i<in_c_f_p_d_d.dimension(0); i++)
2751  for (unsigned int j=0; j<in_c_f_p_d_d.dimension(1); j++)
2752  for (unsigned int k=0; k<in_c_f_p_d_d.dimension(2); k++)
2753  for (unsigned int l=0; l<in_c_f_p_d_d.dimension(3); l++)
2754  for (unsigned int m=0; m<in_c_f_p_d_d.dimension(4); m++){
2755  in_c_f_p_d_d(i,j,k,l,m) = Teuchos::ScalarTraits<double>::random();
2756  }
2757  for (int ic=0; ic < c; ic++) {
2758  for (int ip=0; ip < p; ip++) {
2759  for (int i=0; i<d1; i++) {
2760  for (int j=0; j<d1; j++) {
2761  data_c_p_d_d(ic, ip, i, j) = Teuchos::ScalarTraits<double>::random();
2762  }
2763  }
2764  }
2765  }
2766 RealSpaceTools<double>::inverse(datainv_c_p_d_d, data_c_p_d_d);
2767 RealSpaceTools<double>::transpose(datainvtrn_c_p_d_d, datainv_c_p_d_d);
2768 
2769  for (int ic=0; ic < c; ic++) {
2770  for (int ip=0; ip < 1; ip++) {
2771  for (int i=0; i<d1; i++) {
2772  for (int j=0; j<d1; j++) {
2773  data_c_1_d_d(ic, 0, i, j) = Teuchos::ScalarTraits<double>::random();
2774  }
2775  }
2776  }
2777  }
2778  RealSpaceTools<double>::inverse(datainv_c_1_d_d, data_c_1_d_d);
2779  RealSpaceTools<double>::transpose(datainvtrn_c_1_d_d, datainv_c_1_d_d);
2780 
2781  // Tensor values vary by point:
2782  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d_d, in_c_f_p_d_d);
2783  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainvtrn_c_p_d_d, out_c_f_p_d_d, 't');
2784  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
2785  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
2786  *outStream << "\n\nINCORRECT matmatProductDataField (9): check matrix inverse transpose property\n\n";
2787  errorFlag = -1000;
2788  }
2789  // Tensor values do not vary by point
2790  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d_d, in_c_f_p_d_d);
2791  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainvtrn_c_1_d_d, out_c_f_p_d_d, 't');
2792  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
2793  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
2794  *outStream << "\n\nINCORRECT matmatProductDataField (10): check matrix inverse transpose property\n\n";
2795  errorFlag = -1000;
2796  }
2797  }// test 7.a scope
2798 
2799 
2800 
2801  *outStream \
2802  << "\n"
2803  << "===============================================================================\n"\
2804  << "| TEST 7.b: matmatProductDataField random tests: branch inputFields(F,P,D,D) |\n"\
2805  << "===============================================================================\n";
2806  {// Test 7.b scope
2807  int c=5, p=9, f=7, d1=3;
2808  double zero = INTREPID_TOL*10000.0;
2809 
2810  Kokkos::View<double****> in_f_p_d_d("in_f_p_d_d",f, p, d1, d1);
2811  Kokkos::View<double*****> in_c_f_p_d_d("in_c_f_p_d_d",c, f, p, d1, d1);
2812  Kokkos::View<double**> data_c_p("data_c_p",c, p);
2813  Kokkos::View<double**> datainv_c_p("datainv_c_p",c, p);
2814  Kokkos::View<double**> data_c_1("data_c_1",c, 1);
2815  Kokkos::View<double**> datainv_c_1("datainv_c_1",c, 1);
2816  Kokkos::View<double***> data_c_p_d("data_c_p_d",c, p, d1);
2817  Kokkos::View<double***> datainv_c_p_d("datainv_c_p_d",c, p, d1);
2818  Kokkos::View<double***> data_c_1_d("data_c_1_d",c, 1, d1);
2819  Kokkos::View<double***> datainv_c_1_d("datainv_c_1_d",c, 1, d1);
2820  Kokkos::View<double****> data_c_p_d_d("data_c_p_d_d",c, p, d1, d1);
2821  Kokkos::View<double****> datainv_c_p_d_d("datainv_c_p_d_d",c, p, d1, d1);
2822  Kokkos::View<double****> datainvtrn_c_p_d_d("datainvtrn_c_p_d_d",c, p, d1, d1);
2823  Kokkos::View<double****> data_c_1_d_d("data_c_1_d_d",c, 1, d1, d1);
2824  Kokkos::View<double****> datainv_c_1_d_d("datainv_c_1_d_d",c, 1, d1, d1);
2825  Kokkos::View<double****> datainvtrn_c_1_d_d("datainvtrn_c_1_d_d",c, 1, d1, d1);
2826  Kokkos::View<double**> data_c_p_one("data_c_p_one",c, p);
2827  Kokkos::View<double**> data_c_1_one("data_c_1_one",c, 1);
2828  Kokkos::View<double*****> out_c_f_p_d_d("out_c_f_p_d_d",c, f, p, d1, d1);
2829  Kokkos::View<double*****> outi_c_f_p_d_d("outi_c_f_p_d_d",c, f, p, d1, d1);
2830  /***********************************************************************************************
2831  * Constant diagonal tensor: inputData(C,P) *
2832  **********************************************************************************************/
2833  for (unsigned int i=0; i<in_f_p_d_d.dimension(0); i++)
2834  for (unsigned int j=0; j<in_f_p_d_d.dimension(1); j++)
2835  for (unsigned int k=0; k<in_f_p_d_d.dimension(2); k++)
2836  for (unsigned int l=0; l<in_f_p_d_d.dimension(3); l++){
2837  in_f_p_d_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
2838 
2839  }
2840 
2841  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
2842  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
2843  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
2844  datainv_c_p(i,j) = 1.0 / data_c_p(i,j);
2845  data_c_p_one(i,j) = 1.0;
2846  }
2847 
2848  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
2849  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
2850  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
2851  datainv_c_1(i,j) = 1.0 / data_c_1(i,j);
2852  }
2853 
2854 
2855  // Tensor values vary by point
2856  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
2857  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p, in_f_p_d_d);
2858  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_p, out_c_f_p_d_d);
2859  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
2860  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
2861  *outStream << "\n\nINCORRECT matmatProductDataField (11): check scalar inverse property\n\n";
2862  errorFlag = -1000;
2863  }
2864  // Tensor values do not vary by point
2865  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
2866  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1, in_f_p_d_d);
2867  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_1, out_c_f_p_d_d);
2868  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
2869  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
2870  *outStream << "\n\nINCORRECT matmatProductDataField (12): check scalar inverse property\n\n";
2871  errorFlag = -1000;
2872  }
2873  /***********************************************************************************************
2874  * Non-constant diagonal tensor: inputData(C,P,D) *
2875  **********************************************************************************************/
2876  for (unsigned int i=0; i<in_f_p_d_d.dimension(0); i++)
2877  for (unsigned int j=0; j<in_f_p_d_d.dimension(1); j++)
2878  for (unsigned int k=0; k<in_f_p_d_d.dimension(2); k++)
2879  for (unsigned int l=0; l<in_f_p_d_d.dimension(3); l++){
2880  in_f_p_d_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
2881  }
2882  for (unsigned int i=0; i<data_c_p_d.dimension(0); i++)
2883  for (unsigned int j=0; j<data_c_p_d.dimension(1); j++)
2884  for (unsigned int k=0; k<data_c_p_d.dimension(2); k++){
2885  data_c_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
2886  datainv_c_p_d(i,j,k) = 1.0 / data_c_p_d(i,j,k);
2887  }
2888  for (unsigned int i=0; i<data_c_1_d.dimension(0); i++)
2889  for (unsigned int j=0; j<data_c_1_d.dimension(1); j++)
2890  for (unsigned int k=0; k<data_c_1_d.dimension(2); k++){
2891  data_c_1_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
2892  datainv_c_1_d(i,j,k) = 1.0 / data_c_1_d(i,j,k);
2893  }
2894  // Tensor values vary by point:
2895  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
2896  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d, in_f_p_d_d);
2897  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_p_d, out_c_f_p_d_d);
2898  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
2899  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
2900  *outStream << "\n\nINCORRECT matmatProductDataField (13): check scalar inverse property\n\n";
2901  errorFlag = -1000;
2902  }
2903  // Tensor values do not vary by point:
2904  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
2905  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d, in_f_p_d_d);
2906  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_1_d, out_c_f_p_d_d);
2907  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
2908  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
2909  *outStream << "\n\nINCORRECT matmatProductDataField (14): check scalar inverse property\n\n";
2910  errorFlag = -1000;
2911  }
2912  /***********************************************************************************************
2913  * Full tensor: inputData(C,P,D,D) *
2914  **********************************************************************************************/
2915  for (unsigned int i=0; i<in_f_p_d_d.dimension(0); i++)
2916  for (unsigned int j=0; j<in_f_p_d_d.dimension(1); j++)
2917  for (unsigned int k=0; k<in_f_p_d_d.dimension(2); k++)
2918  for (unsigned int l=0; l<in_f_p_d_d.dimension(3); l++){
2919  in_f_p_d_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
2920  }
2921 
2922 
2923  for (int ic=0; ic < c; ic++) {
2924  for (int ip=0; ip < p; ip++) {
2925  for (int i=0; i<d1; i++) {
2926  for (int j=0; j<d1; j++) {
2927  data_c_p_d_d(ic, ip, i, j) = Teuchos::ScalarTraits<double>::random();
2928  }
2929  }
2930  }
2931  }
2932 RealSpaceTools<double>::inverse(datainv_c_p_d_d, data_c_p_d_d);
2933 
2934  for (int ic=0; ic < c; ic++) {
2935  for (int ip=0; ip < 1; ip++) {
2936  for (int i=0; i<d1; i++) {
2937  for (int j=0; j<d1; j++) {
2938  data_c_1_d_d(ic, 0, i, j) = Teuchos::ScalarTraits<double>::random();
2939  }
2940  }
2941  }
2942  }
2943  RealSpaceTools<double>::inverse(datainv_c_1_d_d, data_c_1_d_d);
2944  // Tensor values vary by point: test "N" and "T" (no-transpose/transpose) options
2945  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
2946  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d_d, in_f_p_d_d);
2947  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_p_d_d, out_c_f_p_d_d);
2948  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
2949  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
2950  *outStream << "\n\nINCORRECT matmatProductDataField (15): check matrix inverse property\n\n";
2951  errorFlag = -1000;
2952  }
2953  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
2954  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d_d, in_f_p_d_d, 't');
2955  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_p_d_d, out_c_f_p_d_d, 't');
2956  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
2957  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
2958  *outStream << "\n\nINCORRECT matmatProductDataField (16): check matrix inverse property, w/ double transpose\n\n";
2959  errorFlag = -1000;
2960  }
2961  // Tensor values do not vary by point: test "N" and "T" (no-transpose/transpose) options
2962  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
2963  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d_d, in_f_p_d_d);
2964  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_1_d_d, out_c_f_p_d_d);
2965  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
2966  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
2967  *outStream << "\n\nINCORRECT matmatProductDataField (17): check matrix inverse property\n\n";
2968  errorFlag = -1000;
2969  }
2970  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
2971  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d_d, in_f_p_d_d, 't');
2972  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainv_c_1_d_d, out_c_f_p_d_d, 't');
2973  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
2974  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
2975  *outStream << "\n\nINCORRECT matmatProductDataField (18): check matrix inverse property, w/ double transpose\n\n";
2976  errorFlag = -1000;
2977  }
2978  /***********************************************************************************************
2979  * Full tensor: inputData(C,P,D,D) test inverse transpose *
2980  **********************************************************************************************/
2981  for (unsigned int i=0; i<in_f_p_d_d.dimension(0); i++)
2982  for (unsigned int j=0; j<in_f_p_d_d.dimension(1); j++)
2983  for (unsigned int k=0; k<in_f_p_d_d.dimension(2); k++)
2984  for (unsigned int l=0; l<in_f_p_d_d.dimension(3); l++){
2985  in_f_p_d_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
2986  }
2987  for (int ic=0; ic < c; ic++) {
2988  for (int ip=0; ip < p; ip++) {
2989  for (int i=0; i<d1; i++) {
2990  for (int j=0; j<d1; j++) {
2991  data_c_p_d_d(ic, ip, i, j) = Teuchos::ScalarTraits<double>::random();
2992  }
2993  }
2994  }
2995  }
2996 RealSpaceTools<double>::inverse(datainv_c_p_d_d, data_c_p_d_d);
2997 RealSpaceTools<double>::transpose(datainvtrn_c_p_d_d, datainv_c_p_d_d);
2998 
2999  for (int ic=0; ic < c; ic++) {
3000  for (int ip=0; ip < 1; ip++) {
3001  for (int i=0; i<d1; i++) {
3002  for (int j=0; j<d1; j++) {
3003  data_c_1_d_d(ic, 0, i, j) = Teuchos::ScalarTraits<double>::random();
3004  }
3005  }
3006  }
3007  }
3008  RealSpaceTools<double>::inverse(datainv_c_1_d_d, data_c_1_d_d);
3009  RealSpaceTools<double>::transpose(datainvtrn_c_1_d_d, datainv_c_1_d_d);
3010  // Tensor values vary by point:
3011  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
3012  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_p_d_d, in_f_p_d_d);
3013  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainvtrn_c_p_d_d, out_c_f_p_d_d, 't');
3014  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
3015  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
3016  *outStream << "\n\nINCORRECT matmatProductDataField (19): check matrix inverse transpose property\n\n";
3017  errorFlag = -1000;
3018  }
3019  // Tensor values do not vary by point:
3020  art::matmatProductDataField<double>(in_c_f_p_d_d, data_c_p_one, in_f_p_d_d);
3021  art::matmatProductDataField<double>(out_c_f_p_d_d, data_c_1_d_d, in_f_p_d_d);
3022  art::matmatProductDataField<double>(outi_c_f_p_d_d, datainvtrn_c_1_d_d, out_c_f_p_d_d, 't');
3023  rst::subtract(outi_c_f_p_d_d, in_c_f_p_d_d);
3024  if (rst::vectorNorm(outi_c_f_p_d_d, NORM_ONE) > zero) {
3025  *outStream << "\n\nINCORRECT matmatProductDataField (20): check matrix inverse transpose property\n\n";
3026  errorFlag = -1000;
3027  }
3028  }// test 7.b scope
3029 
3030 
3031 
3032  *outStream \
3033  << "\n"
3034  << "===============================================================================\n"\
3035  << "| TEST 8.a: matmatProductDataData random tests: branch inputDataRight(C,P,D,D)|\n"\
3036  << "===============================================================================\n";
3037  /*
3038  * Test derived from Test 7.a
3039  */
3040  {// test 8.a scope
3041  int c=5, p=9, d1=3;
3042  double zero = INTREPID_TOL*10000.0;
3043 
3044  Kokkos::View<double****> in_c_p_d_d("in_c_p_d_d",c, p, d1, d1);
3045  Kokkos::View<double****> out_c_p_d_d("out_c_p_d_d",c, p, d1, d1);
3046  Kokkos::View<double****> outi_c_p_d_d("outi_c_p_d_d",c, p, d1, d1);
3047 
3048  Kokkos::View<double**> data_c_p("data_c_p",c, p);
3049  Kokkos::View<double**> datainv_c_p("datainv_c_p",c, p);
3050  Kokkos::View<double**> data_c_1("data_c_1",c, 1);
3051  Kokkos::View<double**> datainv_c_1("datainv_c_1",c, 1);
3052  Kokkos::View<double***> data_c_p_d("data_c_p_d",c, p, d1);
3053  Kokkos::View<double***> datainv_c_p_d("datainv_c_p_d",c, p, d1);
3054  Kokkos::View<double***> data_c_1_d("data_c_1_d",c, 1, d1);
3055  Kokkos::View<double***> datainv_c_1_d("datainv_c_1_d",c, 1, d1);
3056  Kokkos::View<double****> data_c_p_d_d("data_c_p_d_d",c, p, d1, d1);
3057  Kokkos::View<double****> datainv_c_p_d_d("datainv_c_p_d_d",c, p, d1, d1);
3058  Kokkos::View<double****> datainvtrn_c_p_d_d("datainvtrn_c_p_d_d",c, p, d1, d1);
3059  Kokkos::View<double****> data_c_1_d_d("data_c_1_d_d",c, 1, d1, d1);
3060  Kokkos::View<double****> datainv_c_1_d_d("datainv_c_1_d_d",c, 1, d1, d1);
3061  Kokkos::View<double****> datainvtrn_c_1_d_d("datainvtrn_c_1_d_d",c, 1, d1, d1);
3062  /***********************************************************************************************
3063  * Constant diagonal tensor: inputDataLeft(C,P) *
3064  **********************************************************************************************/
3065  for (unsigned int i=0; i<in_c_p_d_d.dimension(0); i++)
3066  for (unsigned int j=0; j<in_c_p_d_d.dimension(1); j++)
3067  for (unsigned int k=0; k<in_c_p_d_d.dimension(2); k++)
3068  for (unsigned int l=0; l<in_c_p_d_d.dimension(3); l++){
3069  in_c_p_d_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
3070  }
3071  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
3072  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
3073  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
3074  datainv_c_p(i,j) = 1.0 / data_c_p(i,j);
3075  }
3076  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
3077  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
3078  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
3079  datainv_c_1(i,j) = 1.0 / data_c_1(i,j);
3080  }
3081 
3082  // Tensor values vary by point:
3083  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p, in_c_p_d_d);
3084  art::matmatProductDataData<double>(out_c_p_d_d, datainv_c_p, out_c_p_d_d);
3085  rst::subtract(out_c_p_d_d, in_c_p_d_d);
3086  if (rst::vectorNorm(out_c_p_d_d, NORM_ONE) > zero) {
3087  *outStream << "\n\nINCORRECT matmatProductDataData (1): check scalar inverse property\n\n";
3088  errorFlag = -1000;
3089  }
3090  // Tensor values do not vary by point:
3091  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1, in_c_p_d_d);
3092  art::matmatProductDataData<double>(out_c_p_d_d, datainv_c_1, out_c_p_d_d);
3093  rst::subtract(out_c_p_d_d, in_c_p_d_d);
3094  if (rst::vectorNorm(out_c_p_d_d, NORM_ONE) > zero) {
3095  *outStream << "\n\nINCORRECT matmatProductDataData (2): check scalar inverse property\n\n";
3096  errorFlag = -1000;
3097  }
3098  /***********************************************************************************************
3099  * Non-onstant diagonal tensor: inputDataLeft(C,P,D) *
3100  **********************************************************************************************/
3101 for (unsigned int i=0; i<in_c_p_d_d.dimension(0); i++)
3102  for (unsigned int j=0; j<in_c_p_d_d.dimension(1); j++)
3103  for (unsigned int k=0; k<in_c_p_d_d.dimension(2); k++)
3104  for (unsigned int l=0; l<in_c_p_d_d.dimension(3); l++){
3105  in_c_p_d_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
3106  }
3107  for (unsigned int i=0; i<data_c_p_d.dimension(0); i++)
3108  for (unsigned int j=0; j<data_c_p_d.dimension(1); j++)
3109  for (unsigned int k=0; k<data_c_p_d.dimension(2); k++){
3110  data_c_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
3111  datainv_c_p_d(i,j,k) = 1.0 / data_c_p_d(i,j,k);
3112  }
3113  for (unsigned int i=0; i<data_c_1_d.dimension(0); i++)
3114  for (unsigned int j=0; j<data_c_1_d.dimension(1); j++)
3115  for (unsigned int k=0; k<data_c_1_d.dimension(2); k++){
3116  data_c_1_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
3117  datainv_c_1_d(i,j,k) = 1.0 / data_c_1_d(i,j,k);
3118  }
3119  // Tensor values vary by point:
3120  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d, in_c_p_d_d);
3121  art::matmatProductDataData<double>(out_c_p_d_d, datainv_c_p_d, out_c_p_d_d);
3122  rst::subtract(out_c_p_d_d, in_c_p_d_d);
3123  if (rst::vectorNorm(out_c_p_d_d,NORM_ONE) > zero) {
3124  *outStream << "\n\nINCORRECT matmatProductDataData (3): check scalar inverse property\n\n";
3125  errorFlag = -1000;
3126  }
3127  // Tensor values do not vary by point
3128  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d, in_c_p_d_d);
3129  art::matmatProductDataData<double>(out_c_p_d_d, datainv_c_1_d, out_c_p_d_d);
3130  rst::subtract(out_c_p_d_d, in_c_p_d_d);
3131  if (rst::vectorNorm(out_c_p_d_d, NORM_ONE) > zero) {
3132  *outStream << "\n\nINCORRECT matmatProductDataData (4): check scalar inverse property\n\n";
3133  errorFlag = -1000;
3134  }
3135  /***********************************************************************************************
3136  * Full tensor: inputDataLeft(C,P,D,D) *
3137  **********************************************************************************************/
3138 for (unsigned int i=0; i<in_c_p_d_d.dimension(0); i++)
3139  for (unsigned int j=0; j<in_c_p_d_d.dimension(1); j++)
3140  for (unsigned int k=0; k<in_c_p_d_d.dimension(2); k++)
3141  for (unsigned int l=0; l<in_c_p_d_d.dimension(3); l++){
3142  in_c_p_d_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
3143  }
3144 for (int ic=0; ic < c; ic++) {
3145  for (int ip=0; ip < p; ip++) {
3146  for (int i=0; i<d1; i++) {
3147  for (int j=0; j<d1; j++) {
3148  data_c_p_d_d(ic, ip, i, j) = Teuchos::ScalarTraits<double>::random();
3149  }
3150  }
3151  }
3152  }
3153 RealSpaceTools<double>::inverse(datainv_c_p_d_d, data_c_p_d_d);
3154 
3155  for (int ic=0; ic < c; ic++) {
3156  for (int ip=0; ip < 1; ip++) {
3157  for (int i=0; i<d1; i++) {
3158  for (int j=0; j<d1; j++) {
3159  data_c_1_d_d(ic, 0, i, j) = Teuchos::ScalarTraits<double>::random();
3160  }
3161  }
3162  }
3163  }
3164  RealSpaceTools<double>::inverse(datainv_c_1_d_d, data_c_1_d_d);
3165 
3166  // Tensor values vary by point: test "N" and "T" options (no transpose/transpose)
3167  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d_d, in_c_p_d_d);
3168  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_p_d_d, out_c_p_d_d);
3169  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
3170  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
3171  *outStream << "\n\nINCORRECT matmatProductDataData (5): check matrix inverse property\n\n";
3172  errorFlag = -1000;
3173  }
3174  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d_d, in_c_p_d_d, 't');
3175  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_p_d_d, out_c_p_d_d, 't');
3176  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
3177  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
3178  *outStream << "\n\nINCORRECT matmatProductDataData (6): check matrix inverse property, w/ double transpose\n\n";
3179  errorFlag = -1000;
3180  }
3181  // Tensor values do not vary by point: test "N" and "T" options (no transpose/transpose)
3182  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d_d, in_c_p_d_d);
3183  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_1_d_d, out_c_p_d_d);
3184  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
3185  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
3186  *outStream << "\n\nINCORRECT matmatProductDataData (7): check matrix inverse property\n\n";
3187  errorFlag = -1000;
3188  }
3189  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d_d, in_c_p_d_d, 't');
3190  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_1_d_d, out_c_p_d_d, 't');
3191  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
3192  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
3193  *outStream << "\n\nINCORRECT matmatProductDataData (8): check matrix inverse property, w/ double transpose\n\n";
3194  errorFlag = -1000;
3195  }
3196  /***********************************************************************************************
3197  * Full tensor: inputDataLeft(C,P,D,D) test inverse transpose *
3198  **********************************************************************************************/
3199 for (unsigned int i=0; i<in_c_p_d_d.dimension(0); i++)
3200  for (unsigned int j=0; j<in_c_p_d_d.dimension(1); j++)
3201  for (unsigned int k=0; k<in_c_p_d_d.dimension(2); k++)
3202  for (unsigned int l=0; l<in_c_p_d_d.dimension(3); l++){
3203  in_c_p_d_d(i,j,k,l) = Teuchos::ScalarTraits<double>::random();
3204  }
3205  for (int ic=0; ic < c; ic++) {
3206  for (int ip=0; ip < p; ip++) {
3207  for (int i=0; i<d1; i++) {
3208  for (int j=0; j<d1; j++) {
3209  data_c_p_d_d(ic, ip, i, j) = Teuchos::ScalarTraits<double>::random();
3210  }
3211  }
3212  }
3213  }
3214 RealSpaceTools<double>::inverse(datainv_c_p_d_d, data_c_p_d_d);
3215 RealSpaceTools<double>::transpose(datainvtrn_c_p_d_d, datainv_c_p_d_d);
3216 
3217  for (int ic=0; ic < c; ic++) {
3218  for (int ip=0; ip < 1; ip++) {
3219  for (int i=0; i<d1; i++) {
3220  for (int j=0; j<d1; j++) {
3221  data_c_1_d_d(ic, 0, i, j) = Teuchos::ScalarTraits<double>::random();
3222  }
3223  }
3224  }
3225  }
3226  RealSpaceTools<double>::inverse(datainv_c_1_d_d, data_c_1_d_d);
3227  RealSpaceTools<double>::transpose(datainvtrn_c_1_d_d, datainv_c_1_d_d);
3228  // Tensor values vary by point:
3229  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d_d, in_c_p_d_d);
3230  art::matmatProductDataData<double>(outi_c_p_d_d, datainvtrn_c_p_d_d, out_c_p_d_d, 't');
3231  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
3232  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
3233  *outStream << "\n\nINCORRECT matmatProductDataData (9): check matrix inverse transpose property\n\n";
3234  errorFlag = -1000;
3235  }
3236  // Tensor values do not vary by point
3237  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d_d, in_c_p_d_d);
3238  art::matmatProductDataData<double>(outi_c_p_d_d, datainvtrn_c_1_d_d, out_c_p_d_d, 't');
3239  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
3240  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
3241  *outStream << "\n\nINCORRECT matmatProductDataData (10): check matrix inverse transpose property\n\n";
3242  errorFlag = -1000;
3243  }
3244  }// test 8.a scope
3245  *outStream \
3246  << "\n"
3247  << "===============================================================================\n"\
3248  << "| TEST 8.b: matmatProductDataData random tests: branch inputDataRight(P,D,D) |\n"\
3249  << "===============================================================================\n";
3250  /*
3251  * Test derived from Test 7.b
3252  */
3253  {// test 8.b scope
3254  int c=5, p=9, d1=3;
3255  double zero = INTREPID_TOL*10000.0;
3256 
3257  Kokkos::View<double***> in_p_d_d("in_p_d_d",p, d1, d1);
3258  Kokkos::View<double****> in_c_p_d_d("in_c_p_d_d",c, p, d1, d1);
3259  Kokkos::View<double****> out_c_p_d_d("out_c_p_d_d",c, p, d1, d1);
3260  Kokkos::View<double****> outi_c_p_d_d("outi_c_p_d_d",c, p, d1, d1);
3261 
3262  Kokkos::View<double**> data_c_p("data_c_p",c, p);
3263  Kokkos::View<double**> datainv_c_p("datainv_c_p",c, p);
3264  Kokkos::View<double**> data_c_1("data_c_1",c, 1);
3265  Kokkos::View<double**> datainv_c_1("datainv_c_1",c, 1);
3266  Kokkos::View<double***> data_c_p_d("data_c_p_d",c, p, d1);
3267  Kokkos::View<double***> datainv_c_p_d("datainv_c_p_d",c, p, d1);
3268  Kokkos::View<double***> data_c_1_d("data_c_1_d",c, 1, d1);
3269  Kokkos::View<double***> datainv_c_1_d("datainv_c_1_d",c, 1, d1);
3270  Kokkos::View<double****> data_c_p_d_d("data_c_p_d_d",c, p, d1, d1);
3271  Kokkos::View<double****> datainv_c_p_d_d("datainv_c_p_d_d",c, p, d1, d1);
3272  Kokkos::View<double****> datainvtrn_c_p_d_d("datainvtrn_c_p_d_d",c, p, d1, d1);
3273  Kokkos::View<double****> data_c_1_d_d("data_c_1_d_d",c, 1, d1, d1);
3274  Kokkos::View<double****> datainv_c_1_d_d("datainv_c_1_d_d",c, 1, d1, d1);
3275  Kokkos::View<double****> datainvtrn_c_1_d_d("datainvtrn_c_1_d_d",c, 1, d1, d1);
3276  Kokkos::View<double**> data_c_p_one("data_c_p_one",c, p);
3277  Kokkos::View<double**> data_c_1_one("data_c_1_one",c, 1);
3278  /***********************************************************************************************
3279  * Constant diagonal tensor: inputData(C,P) *
3280  **********************************************************************************************/
3281 for (unsigned int i=0; i<in_p_d_d.dimension(0); i++)
3282  for (unsigned int j=0; j<in_p_d_d.dimension(1); j++)
3283  for (unsigned int k=0; k<in_p_d_d.dimension(2); k++){
3284  in_p_d_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
3285  }
3286  for (unsigned int i=0; i<data_c_p.dimension(0); i++)
3287  for (unsigned int j=0; j<data_c_p.dimension(1); j++){
3288  data_c_p(i,j) = Teuchos::ScalarTraits<double>::random();
3289  datainv_c_p(i,j) = 1.0 / data_c_p(i,j);
3290  data_c_p_one(i,j) = 1.0;
3291  }
3292 
3293  for (unsigned int i=0; i<data_c_1.dimension(0); i++)
3294  for (unsigned int j=0; j<data_c_1.dimension(1); j++){
3295  data_c_1(i,j) = Teuchos::ScalarTraits<double>::random();
3296  datainv_c_1(i,j) = 1.0 / data_c_1(i,j);
3297  }
3298  // Tensor values vary by point
3299  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3300  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p, in_p_d_d);
3301  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_p, out_c_p_d_d);
3302  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
3303  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
3304  *outStream << "\n\nINCORRECT matmatProductDataData (11): check scalar inverse property\n\n";
3305  errorFlag = -1000;
3306  }
3307  // Tensor values do not vary by point
3308  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3309  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1, in_p_d_d);
3310  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_1, out_c_p_d_d);
3311  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
3312  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
3313  *outStream << "\n\nINCORRECT matmatProductDataData (12): check scalar inverse property\n\n";
3314  errorFlag = -1000;
3315  }
3316  /***********************************************************************************************
3317  * Non-constant diagonal tensor: inputData(C,P,D) *
3318  **********************************************************************************************/
3319 for (unsigned int i=0; i<in_p_d_d.dimension(0); i++)
3320  for (unsigned int j=0; j<in_p_d_d.dimension(1); j++)
3321  for (unsigned int k=0; k<in_p_d_d.dimension(2); k++){
3322  in_p_d_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
3323  }
3324  for (unsigned int i=0; i<data_c_p_d.dimension(0); i++)
3325  for (unsigned int j=0; j<data_c_p_d.dimension(1); j++)
3326  for (unsigned int k=0; k<data_c_p_d.dimension(2); k++){
3327  data_c_p_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
3328  datainv_c_p_d(i,j,k) = 1.0 / data_c_p_d(i,j,k);
3329  }
3330  for (unsigned int i=0; i<data_c_1_d.dimension(0); i++)
3331  for (unsigned int j=0; j<data_c_1_d.dimension(1); j++)
3332  for (unsigned int k=0; k<data_c_1_d.dimension(2); k++){
3333  data_c_1_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
3334  datainv_c_1_d(i,j,k) = 1.0 / data_c_1_d(i,j,k);
3335  }
3336 
3337  // Tensor values vary by point:
3338  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3339  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d, in_p_d_d);
3340  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_p_d, out_c_p_d_d);
3341  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
3342  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
3343  *outStream << "\n\nINCORRECT matmatProductDataData (13): check scalar inverse property\n\n";
3344  errorFlag = -1000;
3345  }
3346  // Tensor values do not vary by point:
3347  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3348  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d, in_p_d_d);
3349  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_1_d, out_c_p_d_d);
3350  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
3351  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
3352  *outStream << "\n\nINCORRECT matmatProductDataData (14): check scalar inverse property\n\n";
3353  errorFlag = -1000;
3354  }
3355  /***********************************************************************************************
3356  * Full tensor: inputData(C,P,D,D) *
3357  **********************************************************************************************/
3358 for (unsigned int i=0; i<in_p_d_d.dimension(0); i++)
3359  for (unsigned int j=0; j<in_p_d_d.dimension(1); j++)
3360  for (unsigned int k=0; k<in_p_d_d.dimension(2); k++){
3361  in_p_d_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
3362  }
3363  for (int ic=0; ic < c; ic++) {
3364  for (int ip=0; ip < p; ip++) {
3365  for (int i=0; i<d1; i++) {
3366  for (int j=0; j<d1; j++) {
3367  data_c_p_d_d(ic, ip, i, j) = Teuchos::ScalarTraits<double>::random();
3368  }
3369  }
3370  }
3371  }
3372 RealSpaceTools<double>::inverse(datainv_c_p_d_d, data_c_p_d_d);
3373 
3374  for (int ic=0; ic < c; ic++) {
3375  for (int ip=0; ip < 1; ip++) {
3376  for (int i=0; i<d1; i++) {
3377  for (int j=0; j<d1; j++) {
3378  data_c_1_d_d(ic, 0, i, j) = Teuchos::ScalarTraits<double>::random();
3379  }
3380  }
3381  }
3382  }
3383  RealSpaceTools<double>::inverse(datainv_c_1_d_d, data_c_1_d_d);
3384 
3385  // Tensor values vary by point: test "N" and "T" (no-transpose/transpose) options
3386  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3387  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d_d, in_p_d_d);
3388  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_p_d_d, out_c_p_d_d);
3389  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
3390  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
3391  *outStream << "\n\nINCORRECT matmatProductDataData (15): check matrix inverse property\n\n";
3392  errorFlag = -1000;
3393  }
3394  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3395  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d_d, in_p_d_d, 't');
3396  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_p_d_d, out_c_p_d_d, 't');
3397  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
3398  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
3399  *outStream << "\n\nINCORRECT matmatProductDataData (16): check matrix inverse property, w/ double transpose\n\n";
3400  errorFlag = -1000;
3401  }
3402  // Tensor values do not vary by point: test "N" and "T" (no-transpose/transpose) options
3403  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3404  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d_d, in_p_d_d);
3405  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_1_d_d, out_c_p_d_d);
3406  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
3407  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
3408  *outStream << "\n\nINCORRECT matmatProductDataData (17): check matrix inverse property\n\n";
3409  errorFlag = -1000;
3410  }
3411  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3412  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d_d, in_p_d_d, 't');
3413  art::matmatProductDataData<double>(outi_c_p_d_d, datainv_c_1_d_d, out_c_p_d_d, 't');
3414  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
3415  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
3416  *outStream << "\n\nINCORRECT matmatProductDataData (18): check matrix inverse property, w/ double transpose\n\n";
3417  errorFlag = -1000;
3418  }
3419  /***********************************************************************************************
3420  * Full tensor: inputData(C,P,D,D) test inverse transpose *
3421  **********************************************************************************************/
3422 for (unsigned int i=0; i<in_p_d_d.dimension(0); i++)
3423  for (unsigned int j=0; j<in_p_d_d.dimension(1); j++)
3424  for (unsigned int k=0; k<in_p_d_d.dimension(2); k++){
3425  in_p_d_d(i,j,k) = Teuchos::ScalarTraits<double>::random();
3426  }
3427  for (int ic=0; ic < c; ic++) {
3428  for (int ip=0; ip < p; ip++) {
3429  for (int i=0; i<d1; i++) {
3430  for (int j=0; j<d1; j++) {
3431  data_c_p_d_d(ic, ip, i, j) = Teuchos::ScalarTraits<double>::random();
3432  }
3433  }
3434  }
3435  }
3436 RealSpaceTools<double>::inverse(datainv_c_p_d_d, data_c_p_d_d);
3437 RealSpaceTools<double>::transpose(datainvtrn_c_p_d_d, datainv_c_p_d_d);
3438 
3439  for (int ic=0; ic < c; ic++) {
3440  for (int ip=0; ip < 1; ip++) {
3441  for (int i=0; i<d1; i++) {
3442  for (int j=0; j<d1; j++) {
3443  data_c_1_d_d(ic, 0, i, j) = Teuchos::ScalarTraits<double>::random();
3444  }
3445  }
3446  }
3447  }
3448  RealSpaceTools<double>::inverse(datainv_c_1_d_d, data_c_1_d_d);
3449  RealSpaceTools<double>::transpose(datainvtrn_c_1_d_d, datainv_c_1_d_d);
3450  // Tensor values vary by point:
3451  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3452  art::matmatProductDataData<double>(out_c_p_d_d, data_c_p_d_d, in_p_d_d);
3453  art::matmatProductDataData<double>(outi_c_p_d_d, datainvtrn_c_p_d_d, out_c_p_d_d, 't');
3454  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
3455  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
3456  *outStream << "\n\nINCORRECT matmatProductDataData (19): check matrix inverse transpose property\n\n";
3457  errorFlag = -1000;
3458  }
3459  // Tensor values do not vary by point:
3460  art::matmatProductDataData<double>(in_c_p_d_d, data_c_p_one, in_p_d_d);
3461  art::matmatProductDataData<double>(out_c_p_d_d, data_c_1_d_d, in_p_d_d);
3462  art::matmatProductDataData<double>(outi_c_p_d_d, datainvtrn_c_1_d_d, out_c_p_d_d, 't');
3463  rst::subtract(outi_c_p_d_d, in_c_p_d_d);
3464  if (rst::vectorNorm(outi_c_p_d_d, NORM_ONE) > zero) {
3465  *outStream << "\n\nINCORRECT matmatProductDataData (20): check matrix inverse transpose property\n\n";
3466  errorFlag = -1000;
3467  }
3468  }// test 8.b scope
3469  }//try
3470 
3471  /*************************************************************************************************
3472  * Finish test *
3473  ************************************************************************************************/
3474 
3475  catch (std::logic_error err) {
3476  *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
3477  *outStream << err.what() << '\n';
3478  *outStream << "-------------------------------------------------------------------------------" << "\n\n";
3479  errorFlag = -1000;
3480  };
3481 
3482 
3483  if (errorFlag != 0)
3484  std::cout << "End Result: TEST FAILED\n";
3485  else
3486  std::cout << "End Result: TEST PASSED\n";
3487 
3488  // reset format state of std::cout
3489  std::cout.copyfmt(oldFormatState);
3490  Kokkos::finalize();
3491 
3492  return errorFlag;
3493 }
Implementation of basic linear algebra functionality in Euclidean space.
static void matvecProductDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const char transpose='N')
There are two use cases: (1) matrix-vector product of a rank-3 container inputDataRight with dimensio...
static void matvecProductDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields, const char transpose='N')
There are two use cases: (1) matrix-vector product of a rank-4 container inputFields with dimensions ...
Header file for utility class to provide multidimensional containers.
Header file for utility class to provide array tools, such as tensor contractions, etc.
static void matmatProductDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight, const char transpose='N')
There are two use cases: (1) matrix-matrix product of a rank-4 container inputDataRight with dimensio...
static void crossProductDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields)
There are two use cases: (1) cross product of a rank-4 container inputFields with dimensions (C...
static void matmatProductDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields, const char transpose='N')
There are two use cases: (1) matrix-matrix product of a rank-5 container inputFields with dimensions ...
Intrepid utilities.
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D...
static void crossProductDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight)
There are two use cases: (1) cross product of a rank-3 container inputDataRight with dimensions (C...
static void outerProductDataData(ArrayOutData &outputData, const ArrayInDataLeft &inputDataLeft, const ArrayInDataRight &inputDataRight)
There are two use cases: (1) outer product of a rank-3 container inputDataRight with dimensions (C...
static void outerProductDataField(ArrayOutFields &outputFields, const ArrayInData &inputData, const ArrayInFields &inputFields)
There are two use cases: (1) outer product of a rank-4 container inputFields with dimensions (C...