|
Intrepid
|
00001 // @HEADER 00002 // ************************************************************************ 00003 // 00004 // Intrepid Package 00005 // Copyright (2007) Sandia Corporation 00006 // 00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 00008 // license for use of this work by or on behalf of the U.S. Government. 00009 // 00010 // This library is free software; you can redistribute it and/or modify 00011 // it under the terms of the GNU Lesser General Public License as 00012 // published by the Free Software Foundation; either version 2.1 of the 00013 // License, or (at your option) any later version. 00014 // 00015 // This library is distributed in the hope that it will be useful, but 00016 // WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00018 // Lesser General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Lesser General Public 00021 // License along with this library; if not, write to the Free Software 00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 00023 // USA 00024 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) or 00025 // Denis Ridzal (dridzal@sandia.gov). 00026 // 00027 // ************************************************************************ 00028 // @HEADER 00029 00035 namespace Intrepid { 00036 00037 template<class Scalar, 00038 class ArrayOutFields, 00039 class ArrayInData, 00040 class ArrayInFields> 00041 void ArrayTools::crossProductDataField(ArrayOutFields & outputFields, 00042 const ArrayInData & inputData, 00043 const ArrayInFields & inputFields){ 00044 00045 #ifdef HAVE_INTREPID_DEBUG 00046 std::string errmsg = ">>> ERROR (ArrayTools::crossProductDataField):"; 00047 /* 00048 * Check array rank and spatial dimension range (if applicable) 00049 * (1) inputData(C,P,D); 00050 * (2) inputFields(C,F,P,D) or (F,P,D); 00051 * (3) outputFields(C,F,P,D) in 3D, or (C,F,P) in 2D 00052 */ 00053 // (1) inputData is (C, P, D) and 2 <= D <= 3 is required 00054 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputData, 3, 3), 00055 std::invalid_argument, errmsg); 00056 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputData, 2, 2,3), 00057 std::invalid_argument, errmsg); 00058 // (2) inputFields is (C, F, P, D) or (F, P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required. 00059 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputFields, 3,4), std::invalid_argument, errmsg); 00060 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputFields, inputFields.rank()-1, 2,3), 00061 std::invalid_argument, errmsg); 00062 // (3) outputFields is (C,F,P,D) in 3D and (C,F,P) in 2D => rank = inputData.dimension(2) + 1 00063 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputFields, inputData.dimension(2)+1, inputData.dimension(2)+1), 00064 std::invalid_argument, errmsg); 00065 /* 00066 * Dimension cross-checks: 00067 * (1) inputData vs. inputFields 00068 * (2) outputFields vs. inputData 00069 * (3) outputFields vs. inputFields 00070 * 00071 * Cross-check (1): 00072 */ 00073 if( inputFields.rank() == 4) { 00074 // inputData(C,P,D) vs. inputFields(C,F,P,D): dimensions C, P, D must match 00075 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputData, 0,1,2, inputFields, 0,2,3), 00076 std::invalid_argument, errmsg); 00077 } 00078 else{ 00079 // inputData(C,P,D) vs. inputFields(F,P,D): dimensions P, D must match 00080 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputData, 1,2, inputFields, 1,2), 00081 std::invalid_argument, errmsg); 00082 } 00083 /* 00084 * Cross-check (2): 00085 */ 00086 if(inputData.dimension(2) == 2) { 00087 // in 2D: outputFields(C,F,P) vs. inputData(C,P,D): dimensions C,P must match 00088 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, 0,2, inputData, 0,1), 00089 std::invalid_argument, errmsg); 00090 } 00091 else{ 00092 // in 3D: outputFields(C,F,P,D) vs. inputData(C,P,D): dimensions C,P,D must match 00093 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, 0,2,3, inputData, 0,1,2), 00094 std::invalid_argument, errmsg); 00095 } 00096 /* 00097 * Cross-check (3): 00098 */ 00099 if(inputData.dimension(2) == 2) { 00100 // In 2D: 00101 if(inputFields.rank() == 4){ 00102 // and rank-4 inputFields: outputFields(C,F,P) vs. inputFields(C,F,P,D): dimensions C,F,P must match 00103 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, 0,1,2, inputFields, 0,1,2), 00104 std::invalid_argument, errmsg); 00105 } 00106 else{ 00107 // and rank-3 inputFields: outputFields(C,F,P) vs. inputFields(F,P,D): dimensions F,P must match 00108 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, 1,2, inputFields, 0,1), 00109 std::invalid_argument, errmsg); 00110 } 00111 } 00112 else{ 00113 // In 3D: 00114 if(inputFields.rank() == 4){ 00115 // and rank-4 inputFields: outputFields(C,F,P,D) vs. inputFields(C,F,P,D): all dimensions C,F,P,D must match 00116 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, inputFields), 00117 std::invalid_argument, errmsg); 00118 } 00119 else{ 00120 // and rank-3 inputFields: outputFields(C,F,P,D) vs. inputFields(F,P,D): dimensions F,P,D must match 00121 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, 1,2,3, inputFields, 0,1,2), 00122 std::invalid_argument, errmsg); 00123 } 00124 } 00125 #endif 00126 // 3D cross product 00127 if(inputData.dimension(2) == 3) { 00128 00129 // inputFields is (C,F,P,D) 00130 if(inputFields.rank() == 4){ 00131 00132 for(int cell = 0; cell < outputFields.dimension(0); cell++){ 00133 for(int field = 0; field < outputFields.dimension(1); field++){ 00134 for(int point = 0; point < outputFields.dimension(2); point++){ 00135 // 00136 outputFields(cell, field, point, 0) = \ 00137 inputData(cell, point, 1)*inputFields(cell, field, point, 2) - 00138 inputData(cell, point, 2)*inputFields(cell, field, point, 1); 00139 // 00140 outputFields(cell, field, point, 1) = \ 00141 inputData(cell, point, 2)*inputFields(cell, field, point, 0) - 00142 inputData(cell, point, 0)*inputFields(cell, field, point, 2); 00143 // 00144 outputFields(cell, field, point, 2) = \ 00145 inputData(cell, point, 0)*inputFields(cell, field, point, 1) - 00146 inputData(cell, point, 1)*inputFields(cell, field, point, 0); 00147 }// point 00148 }// field 00149 } // cell 00150 }// rank = 4 00151 // inputFields is (F,P,D) 00152 else if(inputFields.rank() == 3){ 00153 00154 for(int cell = 0; cell < outputFields.dimension(0); cell++){ 00155 for(int field = 0; field < outputFields.dimension(1); field++){ 00156 for(int point = 0; point < outputFields.dimension(2); point++){ 00157 // 00158 outputFields(cell, field, point, 0) = \ 00159 inputData(cell, point, 1)*inputFields(field, point, 2) - 00160 inputData(cell, point, 2)*inputFields(field, point, 1); 00161 // 00162 outputFields(cell, field, point, 1) = \ 00163 inputData(cell, point, 2)*inputFields(field, point, 0) - 00164 inputData(cell, point, 0)*inputFields(field, point, 2); 00165 // 00166 outputFields(cell, field, point, 2) = \ 00167 inputData(cell, point, 0)*inputFields(field, point, 1) - 00168 inputData(cell, point, 1)*inputFields(field, point, 0); 00169 }// point 00170 }// field 00171 } // cell 00172 }// rank = 3 00173 else{ 00174 TEST_FOR_EXCEPTION( true, std::invalid_argument, 00175 ">>> ERROR (ArrayTools::crossProductDataField): inputFields rank 3 or 4 required.") 00176 } 00177 } 00178 // 2D cross product 00179 else if(inputData.dimension(2) == 2){ 00180 00181 // inputFields is (C,F,P,D) 00182 if(inputFields.rank() == 4){ 00183 00184 for(int cell = 0; cell < outputFields.dimension(0); cell++){ 00185 for(int field = 0; field < outputFields.dimension(1); field++){ 00186 for(int point = 0; point < outputFields.dimension(2); point++){ 00187 outputFields(cell, field, point) = \ 00188 inputData(cell, point, 0)*inputFields(cell, field, point, 1) - 00189 inputData(cell, point, 1)*inputFields(cell, field, point, 0); 00190 }// point 00191 }// field 00192 } // cell 00193 }// rank = 4 00194 // inputFields is (F,P,D) 00195 else if(inputFields.rank() == 3) { 00196 00197 for(int cell = 0; cell < outputFields.dimension(0); cell++){ 00198 for(int field = 0; field < outputFields.dimension(1); field++){ 00199 for(int point = 0; point < outputFields.dimension(2); point++){ 00200 outputFields(cell, field, point) = \ 00201 inputData(cell, point, 0)*inputFields(field, point, 1) - 00202 inputData(cell, point, 1)*inputFields(field, point, 0); 00203 }// point 00204 }// field 00205 } // cell 00206 }// rank = 3 00207 } 00208 // Error: wrong dimension 00209 else { 00210 TEST_FOR_EXCEPTION( true, std::invalid_argument, 00211 ">>> ERROR (ArrayTools::crossProductDataField): spatial dimension 2 or 3 required.") 00212 } 00213 } 00214 00215 00216 00217 template<class Scalar, 00218 class ArrayOutData, 00219 class ArrayInDataLeft, 00220 class ArrayInDataRight> 00221 void ArrayTools::crossProductDataData(ArrayOutData & outputData, 00222 const ArrayInDataLeft & inputDataLeft, 00223 const ArrayInDataRight & inputDataRight){ 00224 #ifdef HAVE_INTREPID_DEBUG 00225 std::string errmsg = ">>> ERROR (ArrayTools::crossProductDataData):"; 00226 /* 00227 * Check array rank and spatial dimension range (if applicable) 00228 * (1) inputDataLeft(C,P,D); 00229 * (2) inputDataRight(C,P,D) or (P,D); 00230 * (3) outputData(C,P,D) in 3D, or (C,P) in 2D 00231 */ 00232 // (1) inputDataLeft is (C, P, D) and 2 <= D <= 3 is required 00233 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataLeft, 3,3), 00234 std::invalid_argument, errmsg); 00235 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataLeft, 2, 2,3), 00236 std::invalid_argument, errmsg); 00237 // (2) inputDataRight is (C, P, D) or (P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required. 00238 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataRight, 2,3), 00239 std::invalid_argument, errmsg); 00240 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataRight, inputDataRight.rank()-1, 2,3), 00241 std::invalid_argument, errmsg); 00242 // (3) outputData is (C,P,D) in 3D and (C,P) in 2D => rank = inputDataLeft.dimension(2) 00243 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputData, inputDataLeft.dimension(2), inputDataLeft.dimension(2)), 00244 std::invalid_argument, errmsg); 00245 /* 00246 * Dimension cross-checks: 00247 * (1) inputDataLeft vs. inputDataRight 00248 * (2) outputData vs. inputDataLeft 00249 * (3) outputData vs. inputDataRight 00250 * 00251 * Cross-check (1): 00252 */ 00253 if( inputDataRight.rank() == 3) { 00254 // inputDataLeft(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match 00255 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputDataLeft, inputDataRight), 00256 std::invalid_argument, errmsg); 00257 } 00258 // inputDataLeft(C, P,D) vs. inputDataRight(P,D): dimensions P, D must match 00259 else{ 00260 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputDataLeft, 1,2, inputDataRight, 0,1), 00261 std::invalid_argument, errmsg); 00262 } 00263 /* 00264 * Cross-check (2): 00265 */ 00266 if(inputDataLeft.dimension(2) == 2){ 00267 // in 2D: outputData(C,P) vs. inputDataLeft(C,P,D): dimensions C, P must match 00268 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputDataLeft, 0,1, outputData, 0,1), 00269 std::invalid_argument, errmsg); 00270 } 00271 else{ 00272 // in 3D: outputData(C,P,D) vs. inputDataLeft(C,P,D): all dimensions C, P, D must match 00273 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputDataLeft, outputData), 00274 std::invalid_argument, errmsg); 00275 } 00276 /* 00277 * Cross-check (3): 00278 */ 00279 if(inputDataLeft.dimension(2) == 2) { 00280 // In 2D: 00281 if(inputDataRight.rank() == 3){ 00282 // and rank-3 inputDataRight: outputData(C,P) vs. inputDataRight(C,P,D): dimensions C,P must match 00283 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputData, 0,1, inputDataRight, 0,1), 00284 std::invalid_argument, errmsg); 00285 } 00286 else{ 00287 // and rank-2 inputDataRight: outputData(C,P) vs. inputDataRight(P,D): dimension P must match 00288 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputData, 1, inputDataRight, 0), 00289 std::invalid_argument, errmsg); 00290 } 00291 } 00292 else{ 00293 // In 3D: 00294 if(inputDataRight.rank() == 3){ 00295 // and rank-3 inputDataRight: outputData(C,P,D) vs. inputDataRight(C,P,D): all dimensions C,P,D must match 00296 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputData, inputDataRight), 00297 std::invalid_argument, errmsg); 00298 } 00299 else{ 00300 // and rank-2 inputDataRight: outputData(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match 00301 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputData, 1,2, inputDataRight, 0,1), 00302 std::invalid_argument, errmsg); 00303 } 00304 } 00305 #endif 00306 // 3D cross product 00307 if(inputDataLeft.dimension(2) == 3) { 00308 00309 // inputDataRight is (C,P,D) 00310 if(inputDataRight.rank() == 3){ 00311 00312 for(int cell = 0; cell < inputDataLeft.dimension(0); cell++){ 00313 for(int point = 0; point < inputDataLeft.dimension(1); point++){ 00314 // 00315 outputData(cell, point, 0) = \ 00316 inputDataLeft(cell, point, 1)*inputDataRight(cell, point, 2) - 00317 inputDataLeft(cell, point, 2)*inputDataRight(cell, point, 1); 00318 // 00319 outputData(cell, point, 1) = \ 00320 inputDataLeft(cell, point, 2)*inputDataRight(cell, point, 0) - 00321 inputDataLeft(cell, point, 0)*inputDataRight(cell, point, 2); 00322 // 00323 outputData(cell, point, 2) = \ 00324 inputDataLeft(cell, point, 0)*inputDataRight(cell, point, 1) - 00325 inputDataLeft(cell, point, 1)*inputDataRight(cell, point, 0); 00326 }// point 00327 } // cell 00328 }// rank = 3 00329 // inputDataRight is (P,D) 00330 else if(inputDataRight.rank() == 2){ 00331 00332 for(int cell = 0; cell < inputDataLeft.dimension(0); cell++){ 00333 for(int point = 0; point < inputDataLeft.dimension(1); point++){ 00334 // 00335 outputData(cell, point, 0) = \ 00336 inputDataLeft(cell, point, 1)*inputDataRight(point, 2) - 00337 inputDataLeft(cell, point, 2)*inputDataRight(point, 1); 00338 // 00339 outputData(cell, point, 1) = \ 00340 inputDataLeft(cell, point, 2)*inputDataRight(point, 0) - 00341 inputDataLeft(cell, point, 0)*inputDataRight(point, 2); 00342 // 00343 outputData(cell, point, 2) = \ 00344 inputDataLeft(cell, point, 0)*inputDataRight(point, 1) - 00345 inputDataLeft(cell, point, 1)*inputDataRight(point, 0); 00346 }// point 00347 } // cell 00348 }// rank = 2 00349 else{ 00350 TEST_FOR_EXCEPTION( true, std::invalid_argument, 00351 ">>> ERROR (ArrayTools::crossProductDataData): inputDataRight rank 2 or 3 required.") 00352 } 00353 } 00354 // 2D cross product 00355 else if(inputDataLeft.dimension(2) == 2){ 00356 00357 // inputDataRight is (C,P,D) 00358 if(inputDataRight.rank() == 3){ 00359 00360 for(int cell = 0; cell < inputDataLeft.dimension(0); cell++){ 00361 for(int point = 0; point < inputDataLeft.dimension(1); point++){ 00362 outputData(cell, point) = \ 00363 inputDataLeft(cell, point, 0)*inputDataRight(cell, point, 1) - 00364 inputDataLeft(cell, point, 1)*inputDataRight(cell, point, 0); 00365 }// point 00366 } // cell 00367 }// rank = 3 00368 // inputDataRight is (P,D) 00369 else if(inputDataRight.rank() == 2) { 00370 00371 for(int cell = 0; cell < inputDataLeft.dimension(0); cell++){ 00372 for(int point = 0; point < inputDataLeft.dimension(1); point++){ 00373 outputData(cell, point) = \ 00374 inputDataLeft(cell, point, 0)*inputDataRight(point, 1) - 00375 inputDataLeft(cell, point, 1)*inputDataRight(point, 0); 00376 }// point 00377 } // cell 00378 }// rank = 2 00379 } 00380 // Error: wrong dimension 00381 else { 00382 TEST_FOR_EXCEPTION( true, std::invalid_argument, 00383 ">>> ERROR (ArrayTools::crossProductDataData): spatial dimension 2 or 3 required.") 00384 } 00385 } 00386 00387 00388 00389 template<class Scalar, 00390 class ArrayOutFields, 00391 class ArrayInData, 00392 class ArrayInFields> 00393 void ArrayTools::outerProductDataField(ArrayOutFields & outputFields, 00394 const ArrayInData & inputData, 00395 const ArrayInFields & inputFields){ 00396 #ifdef HAVE_INTREPID_DEBUG 00397 std::string errmsg = ">>> ERROR (ArrayTools::outerProductDataField):"; 00398 /* 00399 * Check array rank and spatial dimension range (if applicable) 00400 * (1) inputData(C,P,D); 00401 * (2) inputFields(C,F,P,D) or (F,P,D); 00402 * (3) outputFields(C,F,P,D,D) 00403 */ 00404 // (1) inputData is (C, P, D) and 2 <= D <= 3 is required 00405 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputData, 3,3), 00406 std::invalid_argument, errmsg); 00407 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputData, 2, 2,3), 00408 std::invalid_argument, errmsg); 00409 // (2) inputFields is (C, F, P, D) or (F, P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required. 00410 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputFields, 3,4), std::invalid_argument, errmsg); 00411 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputFields, inputFields.rank()-1, 2,3), 00412 std::invalid_argument, errmsg); 00413 // (3) outputFields is (C,F,P,D,D) 00414 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputFields, 5,5), std::invalid_argument, errmsg); 00415 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputFields, 3, 2,3), 00416 std::invalid_argument, errmsg); 00417 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputFields, 4, 2,3), 00418 std::invalid_argument, errmsg); 00419 /* 00420 * Dimension cross-checks: 00421 * (1) inputData vs. inputFields 00422 * (2) outputFields vs. inputData 00423 * (3) outputFields vs. inputFields 00424 * 00425 * Cross-check (2): outputFields(C,F,P,D,D) vs. inputData(C,P,D): dimensions C, P, D must match 00426 */ 00427 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00428 outputFields, 0,2,3,4, 00429 inputData, 0,1,2,2), 00430 std::invalid_argument, errmsg); 00431 /* 00432 * Cross-checks (1,3): 00433 */ 00434 if( inputFields.rank() == 4) { 00435 // Cross-check (1): inputData(C,P,D) vs. inputFields(C,F,P,D): dimensions C, P, D must match 00436 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00437 inputData, 0,1,2, 00438 inputFields, 0,2,3), 00439 std::invalid_argument, errmsg); 00440 // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(C,F,P,D): dimensions C, F, P, D must match 00441 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00442 outputFields, 0,1,2,3,4, 00443 inputFields, 0,1,2,3,3), 00444 std::invalid_argument, errmsg); 00445 } 00446 else{ 00447 // Cross-check (1): inputData(C,P,D) vs. inputFields(F,P,D): dimensions P, D must match 00448 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00449 inputData, 1,2, 00450 inputFields, 1,2), 00451 std::invalid_argument, errmsg); 00452 // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(F,P,D): dimensions F, P, D must match 00453 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00454 outputFields, 1,2,3,4, 00455 inputFields, 0,1,2,2), 00456 std::invalid_argument, errmsg); 00457 } 00458 #endif 00459 00460 // inputFields is (C,F,P,D) 00461 if(inputFields.rank() == 4){ 00462 00463 for(int cell = 0; cell < outputFields.dimension(0); cell++){ 00464 for(int field = 0; field < outputFields.dimension(1); field++){ 00465 for(int point = 0; point < outputFields.dimension(2); point++){ 00466 for(int row = 0; row < outputFields.dimension(3); row++){ 00467 for(int col = 0; col < outputFields.dimension(4); col++){ 00468 outputFields(cell, field, point, row, col) = \ 00469 inputData(cell, point, row)*inputFields(cell, field, point, col); 00470 }// col 00471 }// row 00472 }// point 00473 }// field 00474 } // cell 00475 }// rank = 4 00476 // inputFields is (F,P,D) 00477 else if(inputFields.rank() == 3){ 00478 00479 for(int cell = 0; cell < outputFields.dimension(0); cell++){ 00480 for(int field = 0; field < outputFields.dimension(1); field++){ 00481 for(int point = 0; point < outputFields.dimension(2); point++){ 00482 for(int row = 0; row < outputFields.dimension(3); row++){ 00483 for(int col = 0; col < outputFields.dimension(4); col++){ 00484 outputFields(cell, field, point, row, col) = \ 00485 inputData(cell, point, row)*inputFields(field, point, col); 00486 }// col 00487 }// row 00488 }// point 00489 }// field 00490 } // cell 00491 }// rank = 3 00492 else{ 00493 TEST_FOR_EXCEPTION( true, std::invalid_argument, 00494 ">>> ERROR (ArrayTools::outerProductDataField): inputFields rank 3 or 4 required.") 00495 } 00496 } 00497 00498 00499 00500 template<class Scalar, 00501 class ArrayOutData, 00502 class ArrayInDataLeft, 00503 class ArrayInDataRight> 00504 void ArrayTools::outerProductDataData(ArrayOutData & outputData, 00505 const ArrayInDataLeft & inputDataLeft, 00506 const ArrayInDataRight & inputDataRight){ 00507 #ifdef HAVE_INTREPID_DEBUG 00508 std::string errmsg = ">>> ERROR (ArrayTools::outerProductDataData):"; 00509 /* 00510 * Check array rank and spatial dimension range (if applicable) 00511 * (1) inputDataLeft(C,P,D); 00512 * (2) inputDataRight(C,P,D) or (P,D); 00513 * (3) outputData(C,P,D,D) 00514 */ 00515 // (1) inputDataLeft is (C, P, D) and 2 <= D <= 3 is required 00516 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataLeft, 3,3), 00517 std::invalid_argument, errmsg); 00518 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataLeft, 2, 2,3), 00519 std::invalid_argument, errmsg); 00520 // (2) inputDataRight is (C, P, D) or (P, D) and 2 <= (D=dimension(rank - 1)) <= 3 is required. 00521 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataRight, 2,3), 00522 std::invalid_argument, errmsg); 00523 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataRight, inputDataRight.rank()-1, 2,3), 00524 std::invalid_argument, errmsg); 00525 // (3) outputData is (C,P,D,D) 00526 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputData, 4, 4), std::invalid_argument, errmsg); 00527 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputData, 2, 2,3), 00528 std::invalid_argument, errmsg); 00529 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputData, 3, 2,3), 00530 std::invalid_argument, errmsg); 00531 /* 00532 * Dimension cross-checks: 00533 * (1) inputDataLeft vs. inputDataRight 00534 * (2) outputData vs. inputDataLeft 00535 * (3) outputData vs. inputDataRight 00536 * 00537 * Cross-check (2): outputData(C,P,D,D) vs. inputDataLeft(C,P,D): dimensions C, P, D must match 00538 */ 00539 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00540 outputData, 0,1,2,3, 00541 inputDataLeft, 0,1,2,2), 00542 std::invalid_argument, errmsg); 00543 /* 00544 * Cross-checks (1,3): 00545 */ 00546 if( inputDataRight.rank() == 3) { 00547 // Cross-check (1): inputDataLeft(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match 00548 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, inputDataLeft, inputDataRight), 00549 std::invalid_argument, errmsg); 00550 // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(C,P,D): dimensions C, P, D must match 00551 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00552 outputData, 0,1,2,3, 00553 inputDataRight, 0,1,2,2), 00554 std::invalid_argument, errmsg); 00555 } 00556 else{ 00557 // Cross-check (1): inputDataLeft(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match 00558 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00559 inputDataLeft, 1,2, 00560 inputDataRight, 0,1), 00561 std::invalid_argument, errmsg); 00562 // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(P,D): dimensions P, D must match 00563 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00564 outputData, 1,2,3, 00565 inputDataRight, 0,1,1), 00566 std::invalid_argument, errmsg); 00567 } 00568 #endif 00569 // inputDataRight is (C,P,D) 00570 if(inputDataRight.rank() == 3){ 00571 00572 for(int cell = 0; cell < inputDataLeft.dimension(0); cell++){ 00573 for(int point = 0; point < inputDataLeft.dimension(1); point++){ 00574 for(int row = 0; row < inputDataLeft.dimension(2); row++){ 00575 for(int col = 0; col < inputDataLeft.dimension(2); col++){ 00576 00577 outputData(cell, point, row, col) = \ 00578 inputDataLeft(cell, point, row)*inputDataRight(cell, point, col); 00579 }// col 00580 }// row 00581 }// point 00582 } // cell 00583 }// rank = 3 00584 // inputDataRight is (P,D) 00585 else if(inputDataRight.rank() == 2){ 00586 00587 for(int cell = 0; cell < inputDataLeft.dimension(0); cell++){ 00588 for(int point = 0; point < inputDataLeft.dimension(1); point++){ 00589 for(int row = 0; row < inputDataLeft.dimension(2); row++){ 00590 for(int col = 0; col < inputDataLeft.dimension(2); col++){ 00591 // 00592 outputData(cell, point, row, col) = \ 00593 inputDataLeft(cell, point, row)*inputDataRight(point, col); 00594 } // col 00595 } // row 00596 } // point 00597 } // cell 00598 }// rank = 2 00599 else{ 00600 TEST_FOR_EXCEPTION( true, std::invalid_argument, 00601 ">>> ERROR (ArrayTools::crossProductDataData): inputDataRight rank 2 or 3 required.") 00602 } 00603 } 00604 00605 00606 00607 template<class Scalar, 00608 class ArrayOutFields, 00609 class ArrayInData, 00610 class ArrayInFields> 00611 void ArrayTools::matvecProductDataField(ArrayOutFields & outputFields, 00612 const ArrayInData & inputData, 00613 const ArrayInFields & inputFields, 00614 const char transpose){ 00615 #ifdef HAVE_INTREPID_DEBUG 00616 std::string errmsg = ">>> ERROR (ArrayTools::matvecProductDataField):"; 00617 /* 00618 * Check array rank and spatial dimension range (if applicable) 00619 * (1) inputData(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. data 00620 * (2) inputFields(C,F,P,D) or (F,P,D); 00621 * (3) outputFields(C,F,P,D) 00622 */ 00623 // (1) inputData is (C,P), (C, P, D) or (C, P, D, D) and 1 <= D <= 3 is required 00624 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputData, 2,4), 00625 std::invalid_argument, errmsg); 00626 if(inputData.rank() > 2) { 00627 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputData, 2, 1,3), 00628 std::invalid_argument, errmsg); 00629 } 00630 if(inputData.rank() == 4) { 00631 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputData, 3, 1,3), 00632 std::invalid_argument, errmsg); 00633 } 00634 // (2) inputFields is (C, F, P, D) or (F, P, D) and 1 <= (D=dimension(rank - 1)) <= 3 is required. 00635 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputFields, 3,4), std::invalid_argument, errmsg); 00636 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputFields, inputFields.rank()-1, 1,3), 00637 std::invalid_argument, errmsg); 00638 // (3) outputFields is (C,F,P,D) 00639 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputFields, 4,4), std::invalid_argument, errmsg); 00640 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputFields, 3, 1,3), 00641 std::invalid_argument, errmsg); 00642 /* 00643 * Dimension cross-checks: 00644 * (1) inputData vs. inputFields 00645 * (2) outputFields vs. inputData 00646 * (3) outputFields vs. inputFields 00647 * 00648 * Cross-check (2): outputFields(C,F,P,D) vs. inputData(C,P), (C,P,D) or (C,P,D,D): 00649 * dimensions C, and D must match in all cases, dimension P must match only when non-constant 00650 * data is specified (P>1). Do not check P dimensions with constant data, i.e., when P=1 in 00651 * inputData(C,1,...) 00652 */ 00653 if(inputData.dimension(1) > 1){ // check P dimension if P>1 in inputData 00654 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00655 outputFields, 2, 00656 inputData, 1), 00657 std::invalid_argument, errmsg); 00658 } 00659 if(inputData.rank() == 2) { // inputData(C,P) -> C match 00660 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00661 outputFields, 0, 00662 inputData, 0), 00663 std::invalid_argument, errmsg); 00664 } 00665 if(inputData.rank() == 3){ // inputData(C,P,D) -> C, D match 00666 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00667 outputFields, 0,3, 00668 inputData, 0,2), 00669 std::invalid_argument, errmsg); 00670 00671 } 00672 if(inputData.rank() == 4){ // inputData(C,P,D,D) -> C, D, D match 00673 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00674 outputFields, 0,3,3, 00675 inputData, 0,2,3), 00676 std::invalid_argument, errmsg); 00677 } 00678 /* 00679 * Cross-checks (1,3): 00680 */ 00681 if(inputFields.rank() == 4) { 00682 // 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 00683 if(inputData.dimension(1) > 1){ // check P dimension if P>1 in inputData 00684 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00685 inputFields, 2, 00686 inputData, 1), 00687 std::invalid_argument, errmsg); 00688 } 00689 if(inputData.rank() == 2){ // inputData(C,P) -> C match 00690 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00691 inputData, 0, 00692 inputFields, 0), 00693 std::invalid_argument, errmsg); 00694 } 00695 if(inputData.rank() == 3){ // inputData(C,P,D) -> C, D match 00696 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00697 inputData, 0,2, 00698 inputFields, 0,3), 00699 std::invalid_argument, errmsg); 00700 } 00701 if(inputData.rank() == 4){ // inputData(C,P,D,D) -> C, D, D match 00702 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00703 inputData, 0,2,3, 00704 inputFields, 0,3,3), 00705 std::invalid_argument, errmsg); 00706 } 00707 // Cross-check (3): outputFields(C,F,P,D) vs. inputFields(C,F,P,D): all dimensions C, F, P, D must match 00708 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, inputFields), 00709 std::invalid_argument, errmsg); 00710 } 00711 else{ 00712 // Cross-check (1): inputData(C,P), (C,P,D) or (C,P,D,D) vs. inputFields(F,P,D): dimensions P, D must match 00713 if(inputData.dimension(1) > 1){ // check P if P>1 in inputData 00714 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00715 inputData, 1, 00716 inputFields, 1), 00717 std::invalid_argument, errmsg); 00718 } 00719 if(inputData.rank() == 3){ 00720 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00721 inputData, 2, 00722 inputFields, 2), 00723 std::invalid_argument, errmsg); 00724 } 00725 if(inputData.rank() == 4){ 00726 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00727 inputData, 2,3, 00728 inputFields, 2,2), 00729 std::invalid_argument, errmsg); 00730 } 00731 // Cross-check (3): outputFields(C,F,P,D) vs. inputFields(F,P,D): dimensions F, P, D must match 00732 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 00733 outputFields, 1,2,3, 00734 inputFields, 0,1,2), 00735 std::invalid_argument, errmsg); 00736 } 00737 #endif 00738 int dataRank = inputData.rank(); 00739 int numDataPts = inputData.dimension(1); 00740 int inRank = inputFields.rank(); 00741 int numCells = outputFields.dimension(0); 00742 int numFields = outputFields.dimension(1); 00743 int numPoints = outputFields.dimension(2); 00744 int matDim = outputFields.dimension(3); 00745 /********************************************************************************************* 00746 * inputFields is (C,F,P,D) * 00747 *********************************************************************************************/ 00748 if(inRank == 4){ 00749 if(numDataPts != 1){ // non-constant data 00750 00751 switch(dataRank){ 00752 case 2: 00753 for(int cell = 0; cell < numCells; cell++) { 00754 for(int field = 0; field < numFields; field++) { 00755 for(int point = 0; point < numPoints; point++) { 00756 for( int row = 0; row < matDim; row++) { 00757 outputFields(cell, field, point, row) = \ 00758 inputData(cell, point)*inputFields(cell, field, point, row); 00759 } // Row-loop 00760 } // P-loop 00761 } // F-loop 00762 }// C-loop 00763 break; 00764 00765 case 3: 00766 for(int cell = 0; cell < numCells; cell++) { 00767 for(int field = 0; field < numFields; field++) { 00768 for(int point = 0; point < numPoints; point++) { 00769 for( int row = 0; row < matDim; row++) { 00770 outputFields(cell, field, point, row) = \ 00771 inputData(cell, point, row)*inputFields(cell, field, point, row); 00772 } // Row-loop 00773 } // P-loop 00774 } // F-loop 00775 }// C-loop 00776 break; 00777 00778 case 4: 00779 if ((transpose == 'n') || (transpose == 'N')) { 00780 for(int cell = 0; cell < numCells; cell++){ 00781 for(int field = 0; field < numFields; field++){ 00782 for(int point = 0; point < numPoints; point++){ 00783 for(int row = 0; row < matDim; row++){ 00784 outputFields(cell, field, point, row) = 0.0; 00785 for(int col = 0; col < matDim; col++){ 00786 outputFields(cell, field, point, row) += \ 00787 inputData(cell, point, row, col)*inputFields(cell, field, point, col); 00788 }// col 00789 } //row 00790 }// point 00791 }// field 00792 }// cell 00793 } // no transpose 00794 else if ((transpose == 't') || (transpose == 'T')) { 00795 for(int cell = 0; cell < numCells; cell++){ 00796 for(int field = 0; field < numFields; field++){ 00797 for(int point = 0; point < numPoints; point++){ 00798 for(int row = 0; row < matDim; row++){ 00799 outputFields(cell, field, point, row) = 0.0; 00800 for(int col = 0; col < matDim; col++){ 00801 outputFields(cell, field, point, row) += \ 00802 inputData(cell, point, col, row)*inputFields(cell, field, point, col); 00803 }// col 00804 } //row 00805 }// point 00806 }// field 00807 }// cell 00808 } //transpose 00809 else { 00810 TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument, 00811 ">>> ERROR (ArrayTools::matvecProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'."); 00812 } 00813 break; 00814 00815 default: 00816 TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument, 00817 ">>> ERROR (ArrayTools::matvecProductDataField): inputData rank 2, 3 or 4 required.") 00818 } // switch inputData rank 00819 } 00820 else{ // constant data case 00821 switch(dataRank){ 00822 case 2: 00823 for(int cell = 0; cell < numCells; cell++) { 00824 for(int field = 0; field < numFields; field++) { 00825 for(int point = 0; point < numPoints; point++) { 00826 for( int row = 0; row < matDim; row++) { 00827 outputFields(cell, field, point, row) = \ 00828 inputData(cell, 0)*inputFields(cell, field, point, row); 00829 } // Row-loop 00830 } // P-loop 00831 } // F-loop 00832 }// C-loop 00833 break; 00834 00835 case 3: 00836 for(int cell = 0; cell < numCells; cell++) { 00837 for(int field = 0; field < numFields; field++) { 00838 for(int point = 0; point < numPoints; point++) { 00839 for( int row = 0; row < matDim; row++) { 00840 outputFields(cell, field, point, row) = \ 00841 inputData(cell, 0, row)*inputFields(cell, field, point, row); 00842 } // Row-loop 00843 } // P-loop 00844 } // F-loop 00845 }// C-loop 00846 break; 00847 00848 case 4: 00849 if ((transpose == 'n') || (transpose == 'N')) { 00850 for(int cell = 0; cell < numCells; cell++){ 00851 for(int field = 0; field < numFields; field++){ 00852 for(int point = 0; point < numPoints; point++){ 00853 for(int row = 0; row < matDim; row++){ 00854 outputFields(cell, field, point, row) = 0.0; 00855 for(int col = 0; col < matDim; col++){ 00856 outputFields(cell, field, point, row) += \ 00857 inputData(cell, 0, row, col)*inputFields(cell, field, point, col); 00858 }// col 00859 } //row 00860 }// point 00861 }// field 00862 }// cell 00863 } // no transpose 00864 else if ((transpose == 't') || (transpose == 'T')) { 00865 for(int cell = 0; cell < numCells; cell++){ 00866 for(int field = 0; field < numFields; field++){ 00867 for(int point = 0; point < numPoints; point++){ 00868 for(int row = 0; row < matDim; row++){ 00869 outputFields(cell, field, point, row) = 0.0; 00870 for(int col = 0; col < matDim; col++){ 00871 outputFields(cell, field, point, row) += \ 00872 inputData(cell, 0, col, row)*inputFields(cell, field, point, col); 00873 }// col 00874 } //row 00875 }// point 00876 }// field 00877 }// cell 00878 } //transpose 00879 else { 00880 TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument, 00881 ">>> ERROR (ArrayTools::matvecProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'."); 00882 } 00883 break; 00884 00885 default: 00886 TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument, 00887 ">>> ERROR (ArrayTools::matvecProductDataField): inputData rank 2, 3 or 4 required.") 00888 } // switch inputData rank 00889 } // end constant data case 00890 } // inputFields rank 4 00891 /********************************************************************************************* 00892 * inputFields is (F,P,D) * 00893 *********************************************************************************************/ 00894 else if(inRank == 3) { 00895 if(numDataPts != 1){ // non-constant data 00896 00897 switch(dataRank){ 00898 case 2: 00899 for(int cell = 0; cell < numCells; cell++) { 00900 for(int field = 0; field < numFields; field++) { 00901 for(int point = 0; point < numPoints; point++) { 00902 for( int row = 0; row < matDim; row++) { 00903 outputFields(cell, field, point, row) = \ 00904 inputData(cell, point)*inputFields(field, point, row); 00905 } // Row-loop 00906 } // P-loop 00907 } // F-loop 00908 }// C-loop 00909 break; 00910 00911 case 3: 00912 for(int cell = 0; cell < numCells; cell++) { 00913 for(int field = 0; field < numFields; field++) { 00914 for(int point = 0; point < numPoints; point++) { 00915 for( int row = 0; row < matDim; row++) { 00916 outputFields(cell, field, point, row) = \ 00917 inputData(cell, point, row)*inputFields(field, point, row); 00918 } // Row-loop 00919 } // P-loop 00920 } // F-loop 00921 }// C-loop 00922 break; 00923 00924 case 4: 00925 if ((transpose == 'n') || (transpose == 'N')) { 00926 for(int cell = 0; cell < numCells; cell++){ 00927 for(int field = 0; field < numFields; field++){ 00928 for(int point = 0; point < numPoints; point++){ 00929 for(int row = 0; row < matDim; row++){ 00930 outputFields(cell, field, point, row) = 0.0; 00931 for(int col = 0; col < matDim; col++){ 00932 outputFields(cell, field, point, row) += \ 00933 inputData(cell, point, row, col)*inputFields(field, point, col); 00934 }// col 00935 } //row 00936 }// point 00937 }// field 00938 }// cell 00939 } // no transpose 00940 else if ((transpose == 't') || (transpose == 'T')) { 00941 for(int cell = 0; cell < numCells; cell++){ 00942 for(int field = 0; field < numFields; field++){ 00943 for(int point = 0; point < numPoints; point++){ 00944 for(int row = 0; row < matDim; row++){ 00945 outputFields(cell, field, point, row) = 0.0; 00946 for(int col = 0; col < matDim; col++){ 00947 outputFields(cell, field, point, row) += \ 00948 inputData(cell, point, col, row)*inputFields(field, point, col); 00949 }// col 00950 } //row 00951 }// point 00952 }// field 00953 }// cell 00954 } //transpose 00955 else { 00956 TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument, 00957 ">>> ERROR (ArrayTools::matvecProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'."); 00958 } 00959 break; 00960 00961 default: 00962 TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument, 00963 ">>> ERROR (ArrayTools::matvecProductDataField): inputData rank 2, 3 or 4 required.") 00964 } // switch inputData rank 00965 } 00966 else{ // constant data case 00967 switch(dataRank){ 00968 case 2: 00969 for(int cell = 0; cell < numCells; cell++) { 00970 for(int field = 0; field < numFields; field++) { 00971 for(int point = 0; point < numPoints; point++) { 00972 for( int row = 0; row < matDim; row++) { 00973 outputFields(cell, field, point, row) = \ 00974 inputData(cell, 0)*inputFields(field, point, row); 00975 } // Row-loop 00976 } // P-loop 00977 } // F-loop 00978 }// C-loop 00979 break; 00980 00981 case 3: 00982 for(int cell = 0; cell < numCells; cell++) { 00983 for(int field = 0; field < numFields; field++) { 00984 for(int point = 0; point < numPoints; point++) { 00985 for( int row = 0; row < matDim; row++) { 00986 outputFields(cell, field, point, row) = \ 00987 inputData(cell, 0, row)*inputFields(field, point, row); 00988 } // Row-loop 00989 } // P-loop 00990 } // F-loop 00991 }// C-loop 00992 break; 00993 00994 case 4: 00995 if ((transpose == 'n') || (transpose == 'N')) { 00996 for(int cell = 0; cell < numCells; cell++){ 00997 for(int field = 0; field < numFields; field++){ 00998 for(int point = 0; point < numPoints; point++){ 00999 for(int row = 0; row < matDim; row++){ 01000 outputFields(cell, field, point, row) = 0.0; 01001 for(int col = 0; col < matDim; col++){ 01002 outputFields(cell, field, point, row) += \ 01003 inputData(cell, 0, row, col)*inputFields(field, point, col); 01004 }// col 01005 } //row 01006 }// point 01007 }// field 01008 }// cell 01009 } // no transpose 01010 else if ((transpose == 't') || (transpose == 'T')) { 01011 for(int cell = 0; cell < numCells; cell++){ 01012 for(int field = 0; field < numFields; field++){ 01013 for(int point = 0; point < numPoints; point++){ 01014 for(int row = 0; row < matDim; row++){ 01015 outputFields(cell, field, point, row) = 0.0; 01016 for(int col = 0; col < matDim; col++){ 01017 outputFields(cell, field, point, row) += \ 01018 inputData(cell, 0, col, row)*inputFields(field, point, col); 01019 }// col 01020 } //row 01021 }// point 01022 }// field 01023 }// cell 01024 } //transpose 01025 else { 01026 TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument, 01027 ">>> ERROR (ArrayTools::matvecProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'."); 01028 } 01029 break; 01030 01031 default: 01032 TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument, 01033 ">>> ERROR (ArrayTools::matvecProductDataField): inputData rank 2, 3 or 4 required.") 01034 } // switch inputData rank 01035 } // end constant data case 01036 } // inputFields rank 3 01037 else { 01038 TEST_FOR_EXCEPTION( true, std::invalid_argument, 01039 ">>> ERROR (ArrayTools::matvecProductDataField): inputFields rank 3 or 4 required.") 01040 }// rank error 01041 } 01042 01043 01044 01045 template<class Scalar, 01046 class ArrayOutData, 01047 class ArrayInDataLeft, 01048 class ArrayInDataRight> 01049 void ArrayTools::matvecProductDataData(ArrayOutData & outputData, 01050 const ArrayInDataLeft & inputDataLeft, 01051 const ArrayInDataRight & inputDataRight, 01052 const char transpose){ 01053 #ifdef HAVE_INTREPID_DEBUG 01054 std::string errmsg = ">>> ERROR (ArrayTools::matvecProductDataData):"; 01055 /* 01056 * Check array rank and spatial dimension range (if applicable) 01057 * (1) inputDataLeft(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. left data 01058 * (2) inputDataRight(C,P,D) or (P,D); 01059 * (3) outputData(C,P,D) 01060 */ 01061 // (1) inputDataLeft is (C,P), (C,P,D) or (C,P,D,D) and 1 <= D <= 3 is required 01062 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataLeft, 2,4), 01063 std::invalid_argument, errmsg); 01064 if(inputDataLeft.rank() > 2) { 01065 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataLeft, 2, 1,3), 01066 std::invalid_argument, errmsg); 01067 } 01068 if(inputDataLeft.rank() == 4) { 01069 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataLeft, 3, 1,3), 01070 std::invalid_argument, errmsg); 01071 } 01072 // (2) inputDataRight is (C, P, D) or (P, D) and 1 <= (D=dimension(rank - 1)) <= 3 is required. 01073 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataRight, 2,3), 01074 std::invalid_argument, errmsg); 01075 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataRight, inputDataRight.rank()-1, 1,3), 01076 std::invalid_argument, errmsg); 01077 // (3) outputData is (C,P,D) 01078 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputData, 3,3), std::invalid_argument, errmsg); 01079 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputData, 2, 1,3), 01080 std::invalid_argument, errmsg); 01081 /* 01082 * Dimension cross-checks: 01083 * (1) inputDataLeft vs. inputDataRight 01084 * (2) outputData vs. inputDataLeft 01085 * (3) outputData vs. inputDataRight 01086 * 01087 * Cross-check (2): outputData(C,P,D) vs. inputDataLeft(C,P), (C,P,D) or (C,P,D,D): 01088 * dimensions C, and D must match in all cases, dimension P must match only when non-constant 01089 * data is specified (P>1). Do not check P dimensions with constant left data, i.e., when P=1 in 01090 * inputDataLeft(C,1,...) 01091 */ 01092 if(inputDataLeft.dimension(1) > 1){ // check P dimension if P>1 in inputDataLeft 01093 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01094 outputData, 1, 01095 inputDataLeft, 1), 01096 std::invalid_argument, errmsg); 01097 } 01098 if(inputDataLeft.rank() == 2){ // inputDataLeft(C,P): check C 01099 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01100 outputData, 0, 01101 inputDataLeft, 0), 01102 std::invalid_argument, errmsg); 01103 } 01104 if(inputDataLeft.rank() == 3){ // inputDataLeft(C,P,D): check C and D 01105 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01106 outputData, 0,2, 01107 inputDataLeft, 0,2), 01108 std::invalid_argument, errmsg); 01109 } 01110 if(inputDataLeft.rank() == 4){ // inputDataLeft(C,P,D,D): check C and D 01111 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01112 outputData, 0,2,2, 01113 inputDataLeft, 0,2,3), 01114 std::invalid_argument, errmsg); 01115 } 01116 /* 01117 * Cross-checks (1,3): 01118 */ 01119 if( inputDataRight.rank() == 3) { 01120 // 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 01121 if(inputDataLeft.dimension(1) > 1){ // check P dimension if P>1 in inputDataLeft 01122 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01123 inputDataLeft, 1, 01124 inputDataRight, 1), 01125 std::invalid_argument, errmsg); 01126 } 01127 if(inputDataLeft.rank() == 2){ // inputDataLeft(C,P): check C 01128 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01129 inputDataLeft, 0, 01130 inputDataRight, 0), 01131 std::invalid_argument, errmsg); 01132 } 01133 if(inputDataLeft.rank() == 3){ // inputDataLeft(C,P,D): check C and D 01134 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01135 inputDataLeft, 0,2, 01136 inputDataRight, 0,2), 01137 std::invalid_argument, errmsg); 01138 } 01139 if(inputDataLeft.rank() == 4){ // inputDataLeft(C,P,D,D): check C and D 01140 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01141 inputDataLeft, 0,2,3, 01142 inputDataRight, 0,2,2), 01143 std::invalid_argument, errmsg); 01144 } 01145 01146 // Cross-check (3): outputData(C,P,D) vs. inputDataRight(C,P,D): all dimensions C, P, D must match 01147 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputData, inputDataRight), 01148 std::invalid_argument, errmsg); 01149 } 01150 else{ 01151 // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(P,D): dimensions P, D must match 01152 if(inputDataLeft.dimension(1) > 1){ // check P if P>1 in inputData 01153 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01154 inputDataLeft, 1, 01155 inputDataRight, 0), 01156 std::invalid_argument, errmsg); 01157 } 01158 if(inputDataLeft.rank() == 3){ // inputDataLeft(C,P,D): check D 01159 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01160 inputDataLeft, 2, 01161 inputDataRight, 1), 01162 std::invalid_argument, errmsg); 01163 } 01164 if(inputDataLeft.rank() == 4){ // inputDataLeft(C,P,D,D): check D 01165 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01166 inputDataLeft, 2,3, 01167 inputDataRight, 1,1), 01168 std::invalid_argument, errmsg); 01169 } 01170 // Cross-check (3): outputData(C,P,D) vs. inputDataRight(P,D): dimensions P, D must match 01171 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01172 outputData, 1,2, 01173 inputDataRight, 0,1), 01174 std::invalid_argument, errmsg); 01175 } 01176 #endif 01177 int dataLeftRank = inputDataLeft.rank(); 01178 int numDataLeftPts = inputDataLeft.dimension(1); 01179 int dataRightRank = inputDataRight.rank(); 01180 int numCells = outputData.dimension(0); 01181 int numPoints = outputData.dimension(1); 01182 int matDim = outputData.dimension(2); 01183 01184 /********************************************************************************************* 01185 * inputDataRight is (C,P,D) * 01186 *********************************************************************************************/ 01187 if(dataRightRank == 3){ 01188 if(numDataLeftPts != 1){ // non-constant left data 01189 01190 switch(dataLeftRank){ 01191 case 2: 01192 for(int cell = 0; cell < numCells; cell++) { 01193 for(int point = 0; point < numPoints; point++) { 01194 for( int row = 0; row < matDim; row++) { 01195 outputData(cell, point, row) = \ 01196 inputDataLeft(cell, point)*inputDataRight(cell, point, row); 01197 } // Row-loop 01198 } // P-loop 01199 }// C-loop 01200 break; 01201 01202 case 3: 01203 for(int cell = 0; cell < numCells; cell++) { 01204 for(int point = 0; point < numPoints; point++) { 01205 for( int row = 0; row < matDim; row++) { 01206 outputData(cell, point, row) = \ 01207 inputDataLeft(cell, point, row)*inputDataRight(cell, point, row); 01208 } // Row-loop 01209 } // P-loop 01210 }// C-loop 01211 break; 01212 01213 case 4: 01214 if ((transpose == 'n') || (transpose == 'N')) { 01215 for(int cell = 0; cell < numCells; cell++){ 01216 for(int point = 0; point < numPoints; point++){ 01217 for(int row = 0; row < matDim; row++){ 01218 outputData(cell, point, row) = 0.0; 01219 for(int col = 0; col < matDim; col++){ 01220 outputData(cell, point, row) += \ 01221 inputDataLeft(cell, point, row, col)*inputDataRight(cell, point, col); 01222 }// col 01223 } //row 01224 }// point 01225 }// cell 01226 } // no transpose 01227 else if ((transpose == 't') || (transpose == 'T')) { 01228 for(int cell = 0; cell < numCells; cell++){ 01229 for(int point = 0; point < numPoints; point++){ 01230 for(int row = 0; row < matDim; row++){ 01231 outputData(cell, point, row) = 0.0; 01232 for(int col = 0; col < matDim; col++){ 01233 outputData(cell, point, row) += \ 01234 inputDataLeft(cell, point, col, row)*inputDataRight(cell, point, col); 01235 }// col 01236 } //row 01237 }// point 01238 }// cell 01239 } //transpose 01240 else { 01241 TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument, 01242 ">>> ERROR (ArrayTools::matvecProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'."); 01243 } 01244 break; 01245 01246 default: 01247 TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument, 01248 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft rank 2, 3 or 4 required.") 01249 } // switch inputDataLeft rank 01250 } 01251 else{ // constant data case 01252 switch(dataLeftRank){ 01253 case 2: 01254 for(int cell = 0; cell < numCells; cell++) { 01255 for(int point = 0; point < numPoints; point++) { 01256 for( int row = 0; row < matDim; row++) { 01257 outputData(cell, point, row) = \ 01258 inputDataLeft(cell, 0)*inputDataRight(cell, point, row); 01259 } // Row-loop 01260 } // F-loop 01261 }// C-loop 01262 break; 01263 01264 case 3: 01265 for(int cell = 0; cell < numCells; cell++) { 01266 for(int point = 0; point < numPoints; point++) { 01267 for( int row = 0; row < matDim; row++) { 01268 outputData(cell, point, row) = \ 01269 inputDataLeft(cell, 0, row)*inputDataRight(cell, point, row); 01270 } // Row-loop 01271 } // P-loop 01272 }// C-loop 01273 break; 01274 01275 case 4: 01276 if ((transpose == 'n') || (transpose == 'N')) { 01277 for(int cell = 0; cell < numCells; cell++){ 01278 for(int point = 0; point < numPoints; point++){ 01279 for(int row = 0; row < matDim; row++){ 01280 outputData(cell, point, row) = 0.0; 01281 for(int col = 0; col < matDim; col++){ 01282 outputData(cell, point, row) += \ 01283 inputDataLeft(cell, 0, row, col)*inputDataRight(cell, point, col); 01284 }// col 01285 } //row 01286 }// point 01287 }// cell 01288 } // no transpose 01289 else if ((transpose == 't') || (transpose == 'T')) { 01290 for(int cell = 0; cell < numCells; cell++){ 01291 for(int point = 0; point < numPoints; point++){ 01292 for(int row = 0; row < matDim; row++){ 01293 outputData(cell, point, row) = 0.0; 01294 for(int col = 0; col < matDim; col++){ 01295 outputData(cell, point, row) += \ 01296 inputDataLeft(cell, 0, col, row)*inputDataRight(cell, point, col); 01297 }// col 01298 } //row 01299 }// point 01300 }// cell 01301 } //transpose 01302 else { 01303 TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument, 01304 ">>> ERROR (ArrayTools::matvecProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'."); 01305 } 01306 break; 01307 01308 default: 01309 TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument, 01310 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft rank 2, 3 or 4 required.") 01311 } // switch inputDataLeft rank 01312 } // end constant data case 01313 } // inputDataRight rank 4 01314 /********************************************************************************************* 01315 * inputDataRight is (P,D) * 01316 *********************************************************************************************/ 01317 else if(dataRightRank == 2) { 01318 if(numDataLeftPts != 1){ // non-constant data 01319 01320 switch(dataLeftRank){ 01321 case 2: 01322 for(int cell = 0; cell < numCells; cell++) { 01323 for(int point = 0; point < numPoints; point++) { 01324 for( int row = 0; row < matDim; row++) { 01325 outputData(cell, point, row) = \ 01326 inputDataLeft(cell, point)*inputDataRight(point, row); 01327 } // Row-loop 01328 } // P-loop 01329 }// C-loop 01330 break; 01331 01332 case 3: 01333 for(int cell = 0; cell < numCells; cell++) { 01334 for(int point = 0; point < numPoints; point++) { 01335 for( int row = 0; row < matDim; row++) { 01336 outputData(cell, point, row) = \ 01337 inputDataLeft(cell, point, row)*inputDataRight(point, row); 01338 } // Row-loop 01339 } // P-loop 01340 }// C-loop 01341 break; 01342 01343 case 4: 01344 if ((transpose == 'n') || (transpose == 'N')) { 01345 for(int cell = 0; cell < numCells; cell++){ 01346 for(int point = 0; point < numPoints; point++){ 01347 for(int row = 0; row < matDim; row++){ 01348 outputData(cell, point, row) = 0.0; 01349 for(int col = 0; col < matDim; col++){ 01350 outputData(cell, point, row) += \ 01351 inputDataLeft(cell, point, row, col)*inputDataRight(point, col); 01352 }// col 01353 } //row 01354 }// point 01355 }// cell 01356 } // no transpose 01357 else if ((transpose == 't') || (transpose == 'T')) { 01358 for(int cell = 0; cell < numCells; cell++){ 01359 for(int point = 0; point < numPoints; point++){ 01360 for(int row = 0; row < matDim; row++){ 01361 outputData(cell, point, row) = 0.0; 01362 for(int col = 0; col < matDim; col++){ 01363 outputData(cell, point, row) += \ 01364 inputDataLeft(cell, point, col, row)*inputDataRight(point, col); 01365 }// col 01366 } //row 01367 }// point 01368 }// cell 01369 } //transpose 01370 else { 01371 TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument, 01372 ">>> ERROR (ArrayTools::matvecProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'."); 01373 } 01374 break; 01375 01376 default: 01377 TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument, 01378 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft rank 2, 3 or 4 required.") 01379 } // switch inputDataLeft rank 01380 } 01381 else{ // constant data case 01382 switch(dataLeftRank){ 01383 case 2: 01384 for(int cell = 0; cell < numCells; cell++) { 01385 for(int point = 0; point < numPoints; point++) { 01386 for( int row = 0; row < matDim; row++) { 01387 outputData(cell, point, row) = \ 01388 inputDataLeft(cell, 0)*inputDataRight(point, row); 01389 } // Row-loop 01390 } // P-loop 01391 }// C-loop 01392 break; 01393 01394 case 3: 01395 for(int cell = 0; cell < numCells; cell++) { 01396 for(int point = 0; point < numPoints; point++) { 01397 for( int row = 0; row < matDim; row++) { 01398 outputData(cell, point, row) = \ 01399 inputDataLeft(cell, 0, row)*inputDataRight(point, row); 01400 } // Row-loop 01401 } // P-loop 01402 }// C-loop 01403 break; 01404 01405 case 4: 01406 if ((transpose == 'n') || (transpose == 'N')) { 01407 for(int cell = 0; cell < numCells; cell++){ 01408 for(int point = 0; point < numPoints; point++){ 01409 for(int row = 0; row < matDim; row++){ 01410 outputData(cell, point, row) = 0.0; 01411 for(int col = 0; col < matDim; col++){ 01412 outputData(cell, point, row) += \ 01413 inputDataLeft(cell, 0, row, col)*inputDataRight(point, col); 01414 }// col 01415 } //row 01416 }// point 01417 }// cell 01418 } // no transpose 01419 else if ((transpose == 't') || (transpose == 'T')) { 01420 for(int cell = 0; cell < numCells; cell++){ 01421 for(int point = 0; point < numPoints; point++){ 01422 for(int row = 0; row < matDim; row++){ 01423 outputData(cell, point, row) = 0.0; 01424 for(int col = 0; col < matDim; col++){ 01425 outputData(cell, point, row) += \ 01426 inputDataLeft(cell, 0, col, row)*inputDataRight(point, col); 01427 }// col 01428 } //row 01429 }// point 01430 }// cell 01431 } //transpose 01432 else { 01433 TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument, 01434 ">>> ERROR (ArrayTools::matvecProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'."); 01435 } 01436 break; 01437 01438 default: 01439 TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument, 01440 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataLeft rank 2, 3 or 4 required.") 01441 } // switch inputDataLeft rank 01442 } // end constant inputDataLeft case 01443 } // inputDataRight rank 2 01444 else { 01445 TEST_FOR_EXCEPTION( true, std::invalid_argument, 01446 ">>> ERROR (ArrayTools::matvecProductDataData): inputDataRight rank 2 or 3 required.") 01447 }// rank error 01448 } 01449 01450 01451 01452 template<class Scalar, 01453 class ArrayOutFields, 01454 class ArrayInData, 01455 class ArrayInFields> 01456 void ArrayTools::matmatProductDataField(ArrayOutFields & outputFields, 01457 const ArrayInData & inputData, 01458 const ArrayInFields & inputFields, 01459 const char transpose){ 01460 #ifdef HAVE_INTREPID_DEBUG 01461 std::string errmsg = ">>> ERROR (ArrayTools::matmatProductDataField):"; 01462 /* 01463 * Check array rank and spatial dimension range (if applicable) 01464 * (1) inputData(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. data 01465 * (2) inputFields(C,F,P,D,D) or (F,P,D,D); 01466 * (3) outputFields(C,F,P,D,D) 01467 */ 01468 // (1) inputData is (C,P), (C, P, D) or (C, P, D, D) and 1 <= D <= 3 is required 01469 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputData, 2,4), 01470 std::invalid_argument, errmsg); 01471 if(inputData.rank() > 2) { 01472 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputData, 2, 1,3), 01473 std::invalid_argument, errmsg); 01474 } 01475 if(inputData.rank() == 4) { 01476 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputData, 3, 1,3), 01477 std::invalid_argument, errmsg); 01478 } 01479 // (2) inputFields is (C,F,P,D,D) or (F,P,D,D) and 1 <= (dimension(rank-1), (rank-2)) <= 3 is required. 01480 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputFields, 4,5), std::invalid_argument, errmsg); 01481 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputFields, inputFields.rank()-1, 1,3), 01482 std::invalid_argument, errmsg); 01483 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputFields, inputFields.rank()-2, 1,3), 01484 std::invalid_argument, errmsg); 01485 // (3) outputFields is (C,F,P,D,D) 01486 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputFields, 5,5), std::invalid_argument, errmsg); 01487 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputFields, 3, 1,3), 01488 std::invalid_argument, errmsg); 01489 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputFields, 4, 1,3), 01490 std::invalid_argument, errmsg); 01491 /* 01492 * Dimension cross-checks: 01493 * (1) inputData vs. inputFields 01494 * (2) outputFields vs. inputData 01495 * (3) outputFields vs. inputFields 01496 * 01497 * Cross-check (2): outputFields(C,F,P,D,D) vs. inputData(C,P), (C,P,D) or (C,P,D,D): 01498 * dimensions C, and D must match in all cases, dimension P must match only when non-constant 01499 * data is specified (P>1). Do not check P dimensions with constant data, i.e., when P=1 in 01500 * inputData(C,1,...) 01501 */ 01502 if(inputData.dimension(1) > 1){ // check P dimension if P>1 in inputData 01503 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01504 outputFields, 2, 01505 inputData, 1), 01506 std::invalid_argument, errmsg); 01507 } 01508 if(inputData.rank() == 2) { // inputData(C,P) -> C match 01509 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01510 outputFields, 0, 01511 inputData, 0), 01512 std::invalid_argument, errmsg); 01513 } 01514 if(inputData.rank() == 3){ // inputData(C,P,D) -> C, D match 01515 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01516 outputFields, 0,3,4, 01517 inputData, 0,2,2), 01518 std::invalid_argument, errmsg); 01519 01520 } 01521 if(inputData.rank() == 4){ // inputData(C,P,D,D) -> C, D, D match 01522 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01523 outputFields, 0,3,4, 01524 inputData, 0,2,3), 01525 std::invalid_argument, errmsg); 01526 } 01527 /* 01528 * Cross-checks (1,3): 01529 */ 01530 if( inputFields.rank() == 5) { 01531 // 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 01532 if(inputData.dimension(1) > 1){ // check P dimension if P>1 in inputData 01533 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01534 inputData, 1, 01535 inputFields, 2), 01536 std::invalid_argument, errmsg); 01537 } 01538 if(inputData.rank() == 2){ // inputData(C,P) -> C match 01539 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01540 inputData, 0, 01541 inputFields, 0), 01542 std::invalid_argument, errmsg); 01543 } 01544 if(inputData.rank() == 3){ // inputData(C,P,D) -> C, D match 01545 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01546 inputData, 0,2,2, 01547 inputFields, 0,3,4), 01548 std::invalid_argument, errmsg); 01549 } 01550 if(inputData.rank() == 4){ // inputData(C,P,D,D) -> C, D, D match 01551 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01552 inputData, 0,2,3, 01553 inputFields, 0,3,4), 01554 std::invalid_argument, errmsg); 01555 } 01556 // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(C,F,P,D,D): all dimensions C, F, P, D must match 01557 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputFields, inputFields), 01558 std::invalid_argument, errmsg); 01559 } 01560 else{ 01561 // 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 01562 if(inputData.dimension(1) > 1){ // check P if P>1 in inputData 01563 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01564 inputData, 1, 01565 inputFields, 1), 01566 std::invalid_argument, errmsg); 01567 } 01568 if(inputData.rank() == 3){ 01569 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01570 inputData, 2,2, 01571 inputFields, 2,3), 01572 std::invalid_argument, errmsg); 01573 } 01574 if(inputData.rank() == 4){ 01575 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01576 inputData, 2,3, 01577 inputFields, 2,3), 01578 std::invalid_argument, errmsg); 01579 } 01580 // Cross-check (3): outputFields(C,F,P,D,D) vs. inputFields(F,P,D,D): dimensions F, P, D must match 01581 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01582 outputFields, 1,2,3,4, 01583 inputFields, 0,1,2,3), 01584 std::invalid_argument, errmsg); 01585 } 01586 #endif 01587 int dataRank = inputData.rank(); 01588 int numDataPts = inputData.dimension(1); 01589 int inRank = inputFields.rank(); 01590 int numCells = outputFields.dimension(0); 01591 int numFields = outputFields.dimension(1); 01592 int numPoints = outputFields.dimension(2); 01593 int matDim = outputFields.dimension(3); 01594 01595 /********************************************************************************************* 01596 * inputFields is (C,F,P,D,D) * 01597 *********************************************************************************************/ 01598 if(inRank == 5){ 01599 if(numDataPts != 1){ // non-constant data 01600 01601 switch(dataRank){ 01602 case 2: 01603 for(int cell = 0; cell < numCells; cell++) { 01604 for(int field = 0; field < numFields; field++) { 01605 for(int point = 0; point < numPoints; point++) { 01606 for( int row = 0; row < matDim; row++) { 01607 for( int col = 0; col < matDim; col++) { 01608 outputFields(cell, field, point, row, col) = \ 01609 inputData(cell, point)*inputFields(cell, field, point, row, col); 01610 }// Col-loop 01611 } // Row-loop 01612 } // P-loop 01613 } // F-loop 01614 }// C-loop 01615 break; 01616 01617 case 3: 01618 for(int cell = 0; cell < numCells; cell++) { 01619 for(int field = 0; field < numFields; field++) { 01620 for(int point = 0; point < numPoints; point++) { 01621 for( int row = 0; row < matDim; row++) { 01622 for( int col = 0; col < matDim; col++) { 01623 outputFields(cell, field, point, row, col) = \ 01624 inputData(cell, point, row)*inputFields(cell, field, point, row, col); 01625 }// Col-loop 01626 } // Row-loop 01627 } // P-loop 01628 } // F-loop 01629 }// C-loop 01630 break; 01631 01632 case 4: 01633 if ((transpose == 'n') || (transpose == 'N')) { 01634 for(int cell = 0; cell < numCells; cell++){ 01635 for(int field = 0; field < numFields; field++){ 01636 for(int point = 0; point < numPoints; point++){ 01637 for(int row = 0; row < matDim; row++){ 01638 for(int col = 0; col < matDim; col++){ 01639 outputFields(cell, field, point, row, col) = 0.0; 01640 for(int i = 0; i < matDim; i++){ 01641 outputFields(cell, field, point, row, col) += \ 01642 inputData(cell, point, row, i)*inputFields(cell, field, point, i, col); 01643 }// i 01644 } // col 01645 } //row 01646 }// point 01647 }// field 01648 }// cell 01649 } // no transpose 01650 else if ((transpose == 't') || (transpose == 'T')) { 01651 for(int cell = 0; cell < numCells; cell++){ 01652 for(int field = 0; field < numFields; field++){ 01653 for(int point = 0; point < numPoints; point++){ 01654 for(int row = 0; row < matDim; row++){ 01655 for(int col = 0; col < matDim; col++){ 01656 outputFields(cell, field, point, row, col) = 0.0; 01657 for(int i = 0; i < matDim; i++){ 01658 outputFields(cell, field, point, row, col) += \ 01659 inputData(cell, point, i, row)*inputFields(cell, field, point, i, col); 01660 }// i 01661 } // col 01662 } //row 01663 }// point 01664 }// field 01665 }// cell 01666 } //transpose 01667 else { 01668 TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument, 01669 ">>> ERROR (ArrayTools::matmatProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'."); 01670 } 01671 break; 01672 01673 default: 01674 TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument, 01675 ">>> ERROR (ArrayTools::matmatProductDataField): inputData rank 2, 3 or 4 required.") 01676 } // switch inputData rank 01677 } 01678 else{ // constant data case 01679 switch(dataRank){ 01680 case 2: 01681 for(int cell = 0; cell < numCells; cell++) { 01682 for(int field = 0; field < numFields; field++) { 01683 for(int point = 0; point < numPoints; point++) { 01684 for( int row = 0; row < matDim; row++) { 01685 for( int col = 0; col < matDim; col++) { 01686 outputFields(cell, field, point, row, col) = \ 01687 inputData(cell, 0)*inputFields(cell, field, point, row, col); 01688 }// Col-loop 01689 } // Row-loop 01690 } // P-loop 01691 } // F-loop 01692 }// C-loop 01693 break; 01694 01695 case 3: 01696 for(int cell = 0; cell < numCells; cell++) { 01697 for(int field = 0; field < numFields; field++) { 01698 for(int point = 0; point < numPoints; point++) { 01699 for( int row = 0; row < matDim; row++) { 01700 for( int col = 0; col < matDim; col++) { 01701 outputFields(cell, field, point, row, col) = \ 01702 inputData(cell, 0, row)*inputFields(cell, field, point, row, col); 01703 }// Col-loop 01704 } // Row-loop 01705 } // P-loop 01706 } // F-loop 01707 }// C-loop 01708 break; 01709 01710 case 4: 01711 if ((transpose == 'n') || (transpose == 'N')) { 01712 for(int cell = 0; cell < numCells; cell++){ 01713 for(int field = 0; field < numFields; field++){ 01714 for(int point = 0; point < numPoints; point++){ 01715 for(int row = 0; row < matDim; row++){ 01716 for(int col = 0; col < matDim; col++){ 01717 outputFields(cell, field, point, row, col) = 0.0; 01718 for(int i = 0; i < matDim; i++){ 01719 outputFields(cell, field, point, row, col) += \ 01720 inputData(cell, 0, row, i)*inputFields(cell, field, point, i, col); 01721 }// i 01722 } // col 01723 } //row 01724 }// point 01725 }// field 01726 }// cell 01727 } // no transpose 01728 else if ((transpose == 't') || (transpose == 'T')) { 01729 for(int cell = 0; cell < numCells; cell++){ 01730 for(int field = 0; field < numFields; field++){ 01731 for(int point = 0; point < numPoints; point++){ 01732 for(int row = 0; row < matDim; row++){ 01733 for(int col = 0; col < matDim; col++){ 01734 outputFields(cell, field, point, row, col) = 0.0; 01735 for(int i = 0; i < matDim; i++){ 01736 outputFields(cell, field, point, row, col) += \ 01737 inputData(cell, 0, i, row)*inputFields(cell, field, point, i, col); 01738 }// i 01739 } // col 01740 } //row 01741 }// point 01742 }// field 01743 }// cell 01744 } //transpose 01745 else { 01746 TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument, 01747 ">>> ERROR (ArrayTools::matmatProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'."); 01748 } 01749 break; 01750 01751 default: 01752 TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument, 01753 ">>> ERROR (ArrayTools::matmatProductDataField): inputData rank 2, 3 or 4 required.") 01754 } // switch inputData rank 01755 } // end constant data case 01756 } // inputFields rank 5 01757 /********************************************************************************************** 01758 * inputFields is (F,P,D,D) * 01759 *********************************************************************************************/ 01760 else if(inRank == 4) { 01761 if(numDataPts != 1){ // non-constant data 01762 01763 switch(dataRank){ 01764 case 2: 01765 for(int cell = 0; cell < numCells; cell++) { 01766 for(int field = 0; field < numFields; field++) { 01767 for(int point = 0; point < numPoints; point++) { 01768 for( int row = 0; row < matDim; row++) { 01769 for( int col = 0; col < matDim; col++) { 01770 outputFields(cell, field, point, row, col) = \ 01771 inputData(cell, point)*inputFields(field, point, row, col); 01772 }// Col-loop 01773 } // Row-loop 01774 } // P-loop 01775 } // F-loop 01776 }// C-loop 01777 break; 01778 01779 case 3: 01780 for(int cell = 0; cell < numCells; cell++) { 01781 for(int field = 0; field < numFields; field++) { 01782 for(int point = 0; point < numPoints; point++) { 01783 for( int row = 0; row < matDim; row++) { 01784 for( int col = 0; col < matDim; col++) { 01785 outputFields(cell, field, point, row, col) = \ 01786 inputData(cell, point, row)*inputFields(field, point, row, col); 01787 }// Col-loop 01788 } // Row-loop 01789 } // P-loop 01790 } // F-loop 01791 }// C-loop 01792 break; 01793 01794 case 4: 01795 if ((transpose == 'n') || (transpose == 'N')) { 01796 for(int cell = 0; cell < numCells; cell++){ 01797 for(int field = 0; field < numFields; field++){ 01798 for(int point = 0; point < numPoints; point++){ 01799 for(int row = 0; row < matDim; row++){ 01800 for(int col = 0; col < matDim; col++){ 01801 outputFields(cell, field, point, row, col) = 0.0; 01802 for(int i = 0; i < matDim; i++){ 01803 outputFields(cell, field, point, row, col) += \ 01804 inputData(cell, point, row, i)*inputFields(field, point, i, col); 01805 }// i 01806 } // col 01807 } //row 01808 }// point 01809 }// field 01810 }// cell 01811 } // no transpose 01812 else if ((transpose == 't') || (transpose == 'T')) { 01813 for(int cell = 0; cell < numCells; cell++){ 01814 for(int field = 0; field < numFields; field++){ 01815 for(int point = 0; point < numPoints; point++){ 01816 for(int row = 0; row < matDim; row++){ 01817 for(int col = 0; col < matDim; col++){ 01818 outputFields(cell, field, point, row, col) = 0.0; 01819 for(int i = 0; i < matDim; i++){ 01820 outputFields(cell, field, point, row, col) += \ 01821 inputData(cell, point, i, row)*inputFields(field, point, i, col); 01822 }// i 01823 } // col 01824 } //row 01825 }// point 01826 }// field 01827 }// cell 01828 } //transpose 01829 else { 01830 TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument, 01831 ">>> ERROR (ArrayTools::matmatProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'."); 01832 } 01833 break; 01834 01835 default: 01836 TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument, 01837 ">>> ERROR (ArrayTools::matmatProductDataField): inputData rank 2, 3 or 4 required.") 01838 } // switch inputData rank 01839 } 01840 else{ // constant data case 01841 switch(dataRank){ 01842 case 2: 01843 for(int cell = 0; cell < numCells; cell++) { 01844 for(int field = 0; field < numFields; field++) { 01845 for(int point = 0; point < numPoints; point++) { 01846 for( int row = 0; row < matDim; row++) { 01847 for( int col = 0; col < matDim; col++) { 01848 outputFields(cell, field, point, row, col) = \ 01849 inputData(cell, 0)*inputFields(field, point, row, col); 01850 }// Col-loop 01851 } // Row-loop 01852 } // P-loop 01853 } // F-loop 01854 }// C-loop 01855 break; 01856 01857 case 3: 01858 for(int cell = 0; cell < numCells; cell++) { 01859 for(int field = 0; field < numFields; field++) { 01860 for(int point = 0; point < numPoints; point++) { 01861 for( int row = 0; row < matDim; row++) { 01862 for( int col = 0; col < matDim; col++) { 01863 outputFields(cell, field, point, row, col) = \ 01864 inputData(cell, 0, row)*inputFields(field, point, row, col); 01865 }// Col-loop 01866 } // Row-loop 01867 } // P-loop 01868 } // F-loop 01869 }// C-loop 01870 break; 01871 01872 case 4: 01873 if ((transpose == 'n') || (transpose == 'N')) { 01874 for(int cell = 0; cell < numCells; cell++){ 01875 for(int field = 0; field < numFields; field++){ 01876 for(int point = 0; point < numPoints; point++){ 01877 for(int row = 0; row < matDim; row++){ 01878 for(int col = 0; col < matDim; col++){ 01879 outputFields(cell, field, point, row, col) = 0.0; 01880 for(int i = 0; i < matDim; i++){ 01881 outputFields(cell, field, point, row, col) += \ 01882 inputData(cell, 0, row, i)*inputFields(field, point, i, col); 01883 }// i 01884 } // col 01885 } //row 01886 }// point 01887 }// field 01888 }// cell 01889 } // no transpose 01890 else if ((transpose == 't') || (transpose == 'T')) { 01891 for(int cell = 0; cell < numCells; cell++){ 01892 for(int field = 0; field < numFields; field++){ 01893 for(int point = 0; point < numPoints; point++){ 01894 for(int row = 0; row < matDim; row++){ 01895 for(int col = 0; col < matDim; col++){ 01896 outputFields(cell, field, point, row, col) = 0.0; 01897 for(int i = 0; i < matDim; i++){ 01898 outputFields(cell, field, point, row, col) += \ 01899 inputData(cell, 0, i, row)*inputFields(field, point, i, col); 01900 }// i 01901 } // col 01902 } //row 01903 }// point 01904 }// field 01905 }// cell 01906 } //transpose 01907 else { 01908 TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument, 01909 ">>> ERROR (ArrayTools::matmatProductDataField): The transpose flag must be 'n', 'N', 't' or 'T'."); 01910 } 01911 break; 01912 01913 default: 01914 TEST_FOR_EXCEPTION( !( (dataRank == 2) || (dataRank == 3) || (dataRank == 4) ), std::invalid_argument, 01915 ">>> ERROR (ArrayTools::matmatProductDataField): inputData rank 2, 3 or 4 required.") 01916 } // switch inputData rank 01917 } // end constant data case 01918 } // inputFields rank 4 01919 else { 01920 TEST_FOR_EXCEPTION( true, std::invalid_argument, 01921 ">>> ERROR (ArrayTools::matmatProductDataField): inputFields rank 4 or 5 required.") 01922 }// rank error 01923 } 01924 01925 01926 01927 template<class Scalar, 01928 class ArrayOutData, 01929 class ArrayInDataLeft, 01930 class ArrayInDataRight> 01931 void ArrayTools::matmatProductDataData(ArrayOutData & outputData, 01932 const ArrayInDataLeft & inputDataLeft, 01933 const ArrayInDataRight & inputDataRight, 01934 const char transpose){ 01935 #ifdef HAVE_INTREPID_DEBUG 01936 std::string errmsg = ">>> ERROR (ArrayTools::matmatProductDataData):"; 01937 /* 01938 * Check array rank and spatial dimension range (if applicable) 01939 * (1) inputDataLeft(C,P), (C,P,D) or (C,P,D,D); P=1 is admissible to allow multiply by const. left data 01940 * (2) inputDataRight(C,P,D,D) or (P,D,D); 01941 * (3) outputData(C,P,D,D) 01942 */ 01943 // (1) inputDataLeft is (C,P), (C,P,D) or (C,P,D,D) and 1 <= D <= 3 is required 01944 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataLeft, 2,4), 01945 std::invalid_argument, errmsg); 01946 if(inputDataLeft.rank() > 2) { 01947 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataLeft, 2, 1,3), 01948 std::invalid_argument, errmsg); 01949 } 01950 if(inputDataLeft.rank() == 4) { 01951 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataLeft, 3, 1,3), 01952 std::invalid_argument, errmsg); 01953 } 01954 // (2) inputDataRight is (C,P,D,D) or (P,D,D) and 1 <= (D=dimension(rank-1),(rank-2)) <= 3 is required. 01955 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, inputDataRight, 3,4), 01956 std::invalid_argument, errmsg); 01957 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataRight, inputDataRight.rank()-1, 1,3), 01958 std::invalid_argument, errmsg); 01959 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, inputDataRight, inputDataRight.rank()-2, 1,3), 01960 std::invalid_argument, errmsg); 01961 // (3) outputData is (C,P,D,D) 01962 TEST_FOR_EXCEPTION( !requireRankRange(errmsg, outputData, 4, 4), std::invalid_argument, errmsg); 01963 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputData, 2, 1,3), 01964 std::invalid_argument, errmsg); 01965 TEST_FOR_EXCEPTION( !requireDimensionRange(errmsg, outputData, 3, 1,3), 01966 std::invalid_argument, errmsg); 01967 /* 01968 * Dimension cross-checks: 01969 * (1) inputDataLeft vs. inputDataRight 01970 * (2) outputData vs. inputDataLeft 01971 * (3) outputData vs. inputDataRight 01972 * 01973 * Cross-check (2): outputData(C,P,D,D) vs. inputDataLeft(C,P), (C,P,D) or (C,P,D,D): 01974 * dimensions C, and D must match in all cases, dimension P must match only when non-constant 01975 * data is specified (P>1). Do not check P dimensions with constant left data, i.e., when P=1 in 01976 * inputDataLeft(C,1,...) 01977 */ 01978 if(inputDataLeft.dimension(1) > 1){ // check P dimension if P>1 in inputDataLeft 01979 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01980 outputData, 1, 01981 inputDataLeft, 1), 01982 std::invalid_argument, errmsg); 01983 } 01984 if(inputDataLeft.rank() == 2){ // inputDataLeft(C,P): check C 01985 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01986 outputData, 0, 01987 inputDataLeft, 0), 01988 std::invalid_argument, errmsg); 01989 } 01990 if(inputDataLeft.rank() == 3){ // inputDataLeft(C,P,D): check C and D 01991 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01992 outputData, 0,2,3, 01993 inputDataLeft, 0,2,2), 01994 std::invalid_argument, errmsg); 01995 } 01996 if(inputDataLeft.rank() == 4){ // inputDataLeft(C,P,D,D): check C and D 01997 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 01998 outputData, 0,2,3, 01999 inputDataLeft, 0,2,3), 02000 std::invalid_argument, errmsg); 02001 } 02002 /* 02003 * Cross-checks (1,3): 02004 */ 02005 if( inputDataRight.rank() == 4) { 02006 // 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 02007 if(inputDataLeft.dimension(1) > 1){ // check P dimension if P>1 in inputDataLeft 02008 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 02009 inputDataLeft, 1, 02010 inputDataRight, 1), 02011 std::invalid_argument, errmsg); 02012 } 02013 if(inputDataLeft.rank() == 2){ // inputDataLeft(C,P): check C 02014 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 02015 inputDataLeft, 0, 02016 inputDataRight, 0), 02017 std::invalid_argument, errmsg); 02018 } 02019 if(inputDataLeft.rank() == 3){ // inputDataLeft(C,P,D): check C and D 02020 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 02021 inputDataLeft, 0,2,2, 02022 inputDataRight, 0,2,3), 02023 std::invalid_argument, errmsg); 02024 } 02025 if(inputDataLeft.rank() == 4){ // inputDataLeft(C,P,D,D): check C and D 02026 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 02027 inputDataLeft, 0,2,3, 02028 inputDataRight, 0,2,3), 02029 std::invalid_argument, errmsg); 02030 } 02031 // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(C,P,D,D): all dimensions C, P, D must match 02032 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, outputData, inputDataRight), 02033 std::invalid_argument, errmsg); 02034 } 02035 else{ 02036 // Cross-check (1): inputDataLeft(C,P), (C,P,D), or (C,P,D,D) vs. inputDataRight(P,D,D): dimensions P, D must match 02037 if(inputDataLeft.dimension(1) > 1){ // check P if P>1 in inputData 02038 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 02039 inputDataLeft, 1, 02040 inputDataRight, 0), 02041 std::invalid_argument, errmsg); 02042 } 02043 if(inputDataLeft.rank() == 3){ // inputDataLeft(C,P,D): check D 02044 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 02045 inputDataLeft, 2,2, 02046 inputDataRight, 1,2), 02047 std::invalid_argument, errmsg); 02048 } 02049 if(inputDataLeft.rank() == 4){ // inputDataLeft(C,P,D,D): check D 02050 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 02051 inputDataLeft, 2,3, 02052 inputDataRight, 1,2), 02053 std::invalid_argument, errmsg); 02054 } 02055 // Cross-check (3): outputData(C,P,D,D) vs. inputDataRight(P,D,D): dimensions P, D must match 02056 TEST_FOR_EXCEPTION( !requireDimensionMatch(errmsg, 02057 outputData, 1,2,3, 02058 inputDataRight, 0,1,2), 02059 std::invalid_argument, errmsg); 02060 } 02061 #endif 02062 int dataLeftRank = inputDataLeft.rank(); 02063 int numDataLeftPts = inputDataLeft.dimension(1); 02064 int dataRightRank = inputDataRight.rank(); 02065 int numCells = outputData.dimension(0); 02066 int numPoints = outputData.dimension(1); 02067 int matDim = outputData.dimension(2); 02068 02069 /********************************************************************************************* 02070 * inputDataRight is (C,P,D,D) * 02071 *********************************************************************************************/ 02072 if(dataRightRank == 4){ 02073 if(numDataLeftPts != 1){ // non-constant data 02074 02075 switch(dataLeftRank){ 02076 case 2: 02077 for(int cell = 0; cell < numCells; cell++) { 02078 for(int point = 0; point < numPoints; point++) { 02079 for( int row = 0; row < matDim; row++) { 02080 for( int col = 0; col < matDim; col++) { 02081 outputData(cell, point, row, col) = \ 02082 inputDataLeft(cell, point)*inputDataRight(cell, point, row, col); 02083 }// Col-loop 02084 } // Row-loop 02085 } // P-loop 02086 }// C-loop 02087 break; 02088 02089 case 3: 02090 for(int cell = 0; cell < numCells; cell++) { 02091 for(int point = 0; point < numPoints; point++) { 02092 for( int row = 0; row < matDim; row++) { 02093 for( int col = 0; col < matDim; col++) { 02094 outputData(cell, point, row, col) = \ 02095 inputDataLeft(cell, point, row)*inputDataRight(cell, point, row, col); 02096 }// Col-loop 02097 } // Row-loop 02098 } // P-loop 02099 }// C-loop 02100 break; 02101 02102 case 4: 02103 if ((transpose == 'n') || (transpose == 'N')) { 02104 for(int cell = 0; cell < numCells; cell++){ 02105 for(int point = 0; point < numPoints; point++){ 02106 for(int row = 0; row < matDim; row++){ 02107 for(int col = 0; col < matDim; col++){ 02108 outputData(cell, point, row, col) = 0.0; 02109 for(int i = 0; i < matDim; i++){ 02110 outputData(cell, point, row, col) += \ 02111 inputDataLeft(cell, point, row, i)*inputDataRight(cell, point, i, col); 02112 }// i 02113 } // col 02114 } //row 02115 }// point 02116 }// cell 02117 } // no transpose 02118 else if ((transpose == 't') || (transpose == 'T')) { 02119 for(int cell = 0; cell < numCells; cell++){ 02120 for(int point = 0; point < numPoints; point++){ 02121 for(int row = 0; row < matDim; row++){ 02122 for(int col = 0; col < matDim; col++){ 02123 outputData(cell, point, row, col) = 0.0; 02124 for(int i = 0; i < matDim; i++){ 02125 outputData(cell, point, row, col) += \ 02126 inputDataLeft(cell, point, i, row)*inputDataRight(cell, point, i, col); 02127 }// i 02128 } // col 02129 } //row 02130 }// point 02131 }// cell 02132 } //transpose 02133 else { 02134 TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument, 02135 ">>> ERROR (ArrayTools::matmatProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'."); 02136 } 02137 break; 02138 02139 default: 02140 TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument, 02141 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft rank 2, 3 or 4 required.") 02142 } // switch inputData rank 02143 } 02144 else{ // constant data case 02145 switch(dataLeftRank){ 02146 case 2: 02147 for(int cell = 0; cell < numCells; cell++) { 02148 for(int point = 0; point < numPoints; point++) { 02149 for( int row = 0; row < matDim; row++) { 02150 for( int col = 0; col < matDim; col++) { 02151 outputData(cell, point, row, col) = \ 02152 inputDataLeft(cell, 0)*inputDataRight(cell, point, row, col); 02153 }// Col-loop 02154 } // Row-loop 02155 } // P-loop 02156 }// C-loop 02157 break; 02158 02159 case 3: 02160 for(int cell = 0; cell < numCells; cell++) { 02161 for(int point = 0; point < numPoints; point++) { 02162 for( int row = 0; row < matDim; row++) { 02163 for( int col = 0; col < matDim; col++) { 02164 outputData(cell, point, row, col) = \ 02165 inputDataLeft(cell, 0, row)*inputDataRight(cell, point, row, col); 02166 }// Col-loop 02167 } // Row-loop 02168 } // P-loop 02169 }// C-loop 02170 break; 02171 02172 case 4: 02173 if ((transpose == 'n') || (transpose == 'N')) { 02174 for(int cell = 0; cell < numCells; cell++){ 02175 for(int point = 0; point < numPoints; point++){ 02176 for(int row = 0; row < matDim; row++){ 02177 for(int col = 0; col < matDim; col++){ 02178 outputData(cell, point, row, col) = 0.0; 02179 for(int i = 0; i < matDim; i++){ 02180 outputData(cell, point, row, col) += \ 02181 inputDataLeft(cell, 0, row, i)*inputDataRight(cell, point, i, col); 02182 }// i 02183 } // col 02184 } //row 02185 }// point 02186 }// cell 02187 } // no transpose 02188 else if ((transpose == 't') || (transpose == 'T')) { 02189 for(int cell = 0; cell < numCells; cell++){ 02190 for(int point = 0; point < numPoints; point++){ 02191 for(int row = 0; row < matDim; row++){ 02192 for(int col = 0; col < matDim; col++){ 02193 outputData(cell, point, row, col) = 0.0; 02194 for(int i = 0; i < matDim; i++){ 02195 outputData(cell, point, row, col) += \ 02196 inputDataLeft(cell, 0, i, row)*inputDataRight(cell, point, i, col); 02197 }// i 02198 } // col 02199 } //row 02200 }// point 02201 }// cell 02202 } //transpose 02203 else { 02204 TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument, 02205 ">>> ERROR (ArrayTools::matmatProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'."); 02206 } 02207 break; 02208 02209 default: 02210 TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument, 02211 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft rank 2, 3 or 4 required.") 02212 } // switch inputDataLeft rank 02213 } // end constant data case 02214 } // inputDataRight rank 4 02215 /********************************************************************************************** 02216 * inputDataRight is (P,D,D) * 02217 *********************************************************************************************/ 02218 else if(dataRightRank == 3) { 02219 if(numDataLeftPts != 1){ // non-constant data 02220 02221 switch(dataLeftRank){ 02222 case 2: 02223 for(int cell = 0; cell < numCells; cell++) { 02224 for(int point = 0; point < numPoints; point++) { 02225 for( int row = 0; row < matDim; row++) { 02226 for( int col = 0; col < matDim; col++) { 02227 outputData(cell, point, row, col) = \ 02228 inputDataLeft(cell, point)*inputDataRight(point, row, col); 02229 }// Col-loop 02230 } // Row-loop 02231 } // P-loop 02232 }// C-loop 02233 break; 02234 02235 case 3: 02236 for(int cell = 0; cell < numCells; cell++) { 02237 for(int point = 0; point < numPoints; point++) { 02238 for( int row = 0; row < matDim; row++) { 02239 for( int col = 0; col < matDim; col++) { 02240 outputData(cell, point, row, col) = \ 02241 inputDataLeft(cell, point, row)*inputDataRight(point, row, col); 02242 }// Col-loop 02243 } // Row-loop 02244 } // P-loop 02245 }// C-loop 02246 break; 02247 02248 case 4: 02249 if ((transpose == 'n') || (transpose == 'N')) { 02250 for(int cell = 0; cell < numCells; cell++){ 02251 for(int point = 0; point < numPoints; point++){ 02252 for(int row = 0; row < matDim; row++){ 02253 for(int col = 0; col < matDim; col++){ 02254 outputData(cell, point, row, col) = 0.0; 02255 for(int i = 0; i < matDim; i++){ 02256 outputData(cell, point, row, col) += \ 02257 inputDataLeft(cell, point, row, i)*inputDataRight(point, i, col); 02258 }// i 02259 } // col 02260 } //row 02261 }// point 02262 }// cell 02263 } // no transpose 02264 else if ((transpose == 't') || (transpose == 'T')) { 02265 for(int cell = 0; cell < numCells; cell++){ 02266 for(int point = 0; point < numPoints; point++){ 02267 for(int row = 0; row < matDim; row++){ 02268 for(int col = 0; col < matDim; col++){ 02269 outputData(cell, point, row, col) = 0.0; 02270 for(int i = 0; i < matDim; i++){ 02271 outputData(cell, point, row, col) += \ 02272 inputDataLeft(cell, point, i, row)*inputDataRight(point, i, col); 02273 }// i 02274 } // col 02275 } //row 02276 }// point 02277 }// cell 02278 } //transpose 02279 else { 02280 TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument, 02281 ">>> ERROR (ArrayTools::matmatProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'."); 02282 } 02283 break; 02284 02285 default: 02286 TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument, 02287 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft rank 2, 3 or 4 required.") 02288 } // switch inputDataLeft rank 02289 } 02290 else{ // constant data case 02291 switch(dataLeftRank){ 02292 case 2: 02293 for(int cell = 0; cell < numCells; cell++) { 02294 for(int point = 0; point < numPoints; point++) { 02295 for( int row = 0; row < matDim; row++) { 02296 for( int col = 0; col < matDim; col++) { 02297 outputData(cell, point, row, col) = \ 02298 inputDataLeft(cell, 0)*inputDataRight(point, row, col); 02299 }// Col-loop 02300 } // Row-loop 02301 } // P-loop 02302 }// C-loop 02303 break; 02304 02305 case 3: 02306 for(int cell = 0; cell < numCells; cell++) { 02307 for(int point = 0; point < numPoints; point++) { 02308 for( int row = 0; row < matDim; row++) { 02309 for( int col = 0; col < matDim; col++) { 02310 outputData(cell, point, row, col) = \ 02311 inputDataLeft(cell, 0, row)*inputDataRight(point, row, col); 02312 }// Col-loop 02313 } // Row-loop 02314 } // P-loop 02315 }// C-loop 02316 break; 02317 02318 case 4: 02319 if ((transpose == 'n') || (transpose == 'N')) { 02320 for(int cell = 0; cell < numCells; cell++){ 02321 for(int point = 0; point < numPoints; point++){ 02322 for(int row = 0; row < matDim; row++){ 02323 for(int col = 0; col < matDim; col++){ 02324 outputData(cell, point, row, col) = 0.0; 02325 for(int i = 0; i < matDim; i++){ 02326 outputData(cell, point, row, col) += \ 02327 inputDataLeft(cell, 0, row, i)*inputDataRight(point, i, col); 02328 }// i 02329 } // col 02330 } //row 02331 }// point 02332 }// cell 02333 } // no transpose 02334 else if ((transpose == 't') || (transpose == 'T')) { 02335 for(int cell = 0; cell < numCells; cell++){ 02336 for(int point = 0; point < numPoints; point++){ 02337 for(int row = 0; row < matDim; row++){ 02338 for(int col = 0; col < matDim; col++){ 02339 outputData(cell, point, row, col) = 0.0; 02340 for(int i = 0; i < matDim; i++){ 02341 outputData(cell, point, row, col) += \ 02342 inputDataLeft(cell, 0, i, row)*inputDataRight(point, i, col); 02343 }// i 02344 } // col 02345 } //row 02346 }// point 02347 }// cell 02348 } //transpose 02349 else { 02350 TEST_FOR_EXCEPTION( !( (transpose == 'n') || (transpose == 'N') || (transpose == 't') || (transpose == 'T') ), std::invalid_argument, 02351 ">>> ERROR (ArrayTools::matmatProductDataData): The transpose flag must be 'n', 'N', 't' or 'T'."); 02352 } 02353 break; 02354 02355 default: 02356 TEST_FOR_EXCEPTION( !( (dataLeftRank == 2) || (dataLeftRank == 3) || (dataLeftRank == 4) ), std::invalid_argument, 02357 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataLeft rank 2, 3 or 4 required.") 02358 } // switch inputDataLeft rank 02359 } // end constant data case 02360 } // inputDataRight rank 3 02361 else { 02362 TEST_FOR_EXCEPTION( true, std::invalid_argument, 02363 ">>> ERROR (ArrayTools::matmatProductDataData): inputDataRight rank 3 or 4 required.") 02364 }// rank error 02365 } 02366 02367 02368 02369 02370 02371 } // end namespace Intrepid 02372 02373 02374 02375 02376 02377 02378 02379 02380 02381 02382 02383 02384 02385 02386 02387 02388 02389 02390 02391 02392
1.7.4