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"
50 #include "Intrepid_FieldContainer.hpp"
51 #include "Intrepid_RealSpaceTools.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 }
Intrepid::RealSpaceTools
Implementation of basic linear algebra functionality in Euclidean space.
Definition: Intrepid_RealSpaceTools.hpp:68
Intrepid::ArrayTools::matmatProductDataData
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...
Definition: Intrepid_ArrayToolsDefTensor.hpp:1976
Intrepid_FieldContainer.hpp
Header file for utility class to provide multidimensional containers.
Intrepid_ArrayTools.hpp
Header file for utility class to provide array tools, such as tensor contractions, etc.
Intrepid::ArrayTools::crossProductDataField
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...
Definition: Intrepid_ArrayToolsDefTensor.hpp:55
Intrepid_Utils.hpp
Intrepid utilities.
Intrepid_RealSpaceTools.hpp
Header file for classes providing basic linear algebra functionality in 1D, 2D and 3D...
Intrepid::ArrayTools
Utility class that provides methods for higher-order algebraic manipulation of user-defined arrays...
Definition: Intrepid_ArrayTools.hpp:71
Intrepid::ArrayTools::crossProductDataData
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...
Definition: Intrepid_ArrayToolsDefTensor.hpp:238
Intrepid::ArrayTools::matmatProductDataField
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 ...
Definition: Intrepid_ArrayToolsDefTensor.hpp:1497
Intrepid::ArrayTools::outerProductDataData
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...
Definition: Intrepid_ArrayToolsDefTensor.hpp:532
Intrepid::ArrayTools::outerProductDataField
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...
Definition: Intrepid_ArrayToolsDefTensor.hpp:417
Intrepid::ArrayTools::matvecProductDataData
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...
Definition: Intrepid_ArrayToolsDefTensor.hpp:1085
Intrepid::ArrayTools::matvecProductDataField
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 ...
Definition: Intrepid_ArrayToolsDefTensor.hpp:642
Generated by   doxygen 1.8.5